text:
module HessianTest

import Calculus

using Base.Test
using ForwardDiff
using ForwardDiff: Dual, Tag
using StaticArrays

include(joinpath(dirname(@__FILE__), "utils.jl"))

#############################
# rosenbrock hardcoded test #
#############################

f = DiffBase.rosenbrock_1
x = [0.1, 0.2, 0.3]
v = f(x)
g = [-9.4, 15.6, 52.0]
h = [-66.0  -40.0    0.0;
     -40.0  130.0  -80.0;
       0.0  -80.0  200.0]

for c in (1, 2, 3), tag in (nothing, f)
    println("  ...running hardcoded test with chunk size = $c and tag = $tag")
    cfg = ForwardDiff.HessianConfig(tag, x, ForwardDiff.Chunk{c}())
    resultcfg = ForwardDiff.HessianConfig(tag, DiffBase.HessianResult(x), x, ForwardDiff.Chunk{c}())

    D = Dual{typeof(Tag(Void, eltype(x))), eltype(x), c}
    @test eltype(cfg) == Dual{typeof(Tag(typeof(tag), Dual{Void,eltype(x),0})), D, c}
    @test eltype(resultcfg) == eltype(cfg)

    @test isapprox(h, ForwardDiff.hessian(f, x))
    @test isapprox(h, ForwardDiff.hessian(f, x, cfg))

    out = similar(x, 3, 3)
    ForwardDiff.hessian!(out, f, x)
    @test isapprox(out, h)

    out = similar(x, 3, 3)
    ForwardDiff.hessian!(out, f, x, cfg)
    @test isapprox(out, h)

    out = DiffBase.HessianResult(x)
    ForwardDiff.hessian!(out, f, x)
    @test isapprox(DiffBase.value(out), v)
    @test isapprox(DiffBase.gradient(out), g)
    @test isapprox(DiffBase.hessian(out), h)

    out = DiffBase.HessianResult(x)
    ForwardDiff.hessian!(out, f, x, resultcfg)
    @test isapprox(DiffBase.value(out), v)
    @test isapprox(DiffBase.gradient(out), g)
    @test isapprox(DiffBase.hessian(out), h)
end

########################
# test vs. Calculus.jl #
########################

for f in DiffBase.VECTOR_TO_NUMBER_FUNCS
    v = f(X)
    g = ForwardDiff.gradient(f, X)
    h = ForwardDiff.hessian(f, X)
    # finite difference approximation error is really bad for Hessians...
    @test isapprox(h, Calculus.hessian(f, X), atol=0.02)
    for c in CHUNK_SIZES, tag in (nothing, f)
        println("  ...testing $f with chunk size = $c and tag = $tag")
        cfg = ForwardDiff.HessianConfig(tag, X, ForwardDiff.Chunk{c}())
        resultcfg = ForwardDiff.HessianConfig(tag, DiffBase.HessianResult(X), X, ForwardDiff.Chunk{c}())

        out = ForwardDiff.hessian(f, X, cfg)
        @test isapprox(out, h)

        out = similar(X, length(X), length(X))
        ForwardDiff.hessian!(out, f, X, cfg)
        @test isapprox(out, h)

        out = DiffBase.HessianResult(X)
        ForwardDiff.hessian!(out, f, X, resultcfg)
        @test isapprox(DiffBase.value(out), v)
        @test isapprox(DiffBase.gradient(out), g)
        @test isapprox(DiffBase.hessian(out), h)
    end
end

##########################################
# test specialized StaticArray codepaths #
##########################################

x = rand(3, 3)
sx = StaticArrays.SArray{Tuple{3,3}}(x)
out = similar(x, 9, 9)
actual = ForwardDiff.hessian(prod, x)

@test ForwardDiff.hessian(prod, sx) == actual

ForwardDiff.hessian!(out, prod, sx)

@test out == actual

result = DiffBase.HessianResult(x)
sresult = DiffBase.HessianResult(sx)

ForwardDiff.hessian!(result, prod, x)
ForwardDiff.hessian!(sresult, prod, sx)

@test DiffBase.value(sresult) == DiffBase.value(result)
@test DiffBase.gradient(sresult) == DiffBase.gradient(result)
@test DiffBase.hessian(sresult) == DiffBase.hessian(result)

end # module


text:
module UProbes

using LLVM
import Libdl

export @probe, @query

# Links:
# * https://lwn.net/Articles/753601/
# * systemtap-sdt-dev for dtrace util on Linux
# * GDB support https://sourceware.org/gdb/onlinedocs/gdb/Static-Probe-Points.html
# * Interesting bcc utils https://github.com/iovisor/bcc/pull/774
# * https://github.com/opendtrace/toolkit
# * bpftrace: https://github.com/iovisor/bpftrace/blob/master/docs/reference_guide.md
# * Both these tools create elf files dynamically
#   * https://github.com/chrisa/libusdt
#   * https://github.com/sthima/libstapsdt
# * https://cuviper.github.io/rust-libprobe/doc/src/probe/src/platform/systemtap.rs.html
# * https://sourceware.org/systemtap/man/stapprobes.3stap.html#lbAO (see `process.mark`)
# * https://sourceware.org/systemtap/wiki/AddingUserSpaceProbingToApps
# * https://sourceware.org/systemtap/wiki/UserSpaceProbeImplementation
# * https://sourceware.org/gdb/onlinedocs/gdb/Static-Probe-Points.html

if Sys.iswindows()
macro probe(args...)
    :()
end
macro quert(args...)
    :(false)
end

else

macro probe(provider, name, args...)
    quote
        $__probe(Val($provider), Val($name), $(map(esc, args)...))
    end
end

macro query(provider, name, types...)
    quote
        $__query(Val($provider), Val($name), Tuple{$(map(esc, types)...)})
    end
end

end # iswindows

@generated function __probe(::Val{provider}, ::Val{name}, args...) where {provider, name}
    dlptr = cache_dl(provider, name, args)
    dlsym = Libdl.dlsym(dlptr, Symbol(join(("__uprobe", provider, name), "_")))
    quote
        ccall($dlsym, Nothing, ($(args...),), $((:(args[$i]) for i in 1:length(args))...))
    end
end

# If a semaphore is associated with a probe, it will be of type unsigned short.
# A semaphore may gate invocations of a probe; it must be set to a non-zero
# value to guarantee that the probe will be hit. Semaphores are treated
# as a counter; your tool should increment the semaphore to enable it,
# and decrement the semaphore when finished.
@generated function __query(::Val{provider}, ::Val{name}, ::Type{args}) where {provider, name, args}
    args = (args.parameters..., )
    dlptr = cache_dl(provider, name, args)
    dlsym = Libdl.dlsym(dlptr, Symbol(join((provider, name, "semaphore"), "_")))
    quote
        unsafe_load(convert(Ptr{UInt16}, $dlsym)) !== UInt16(0)
    end
end

const __probes = Dict{Tuple{Symbol, Symbol, Tuple}, Ptr{Nothing}}()
function cache_dl(provider, name, args)
    key = (provider, name, args)
    if !haskey(__probes, key)
        dlptr = emit_probe(provider, name, args)
        __probes[key] = dlptr
        return dlptr
    end
    return __probes[key]
end

argstring(::Type{T}, i) where T = error("UProbes.jl doesn't know how to handle $T")
argstring(::Type{T}, i) where T <: Signed   = string("-", sizeof(T), raw"@$", i)
argstring(::Type{T}, i) where T <: Unsigned = string(     sizeof(T), raw"@$", i)

if Sys.ARCH === :powerpc64le || Sys.ARCH === :ppc64le
    constraint(::Type{T}) where T <: Integer = "nZr"
else
    constraint(::Type{T}) where T <: Integer = "nor"
end

function emit_probe(provider::Symbol, name::Symbol, args::Tuple; file = tempname())
    @assert length(args) <= 12

    argstr = join((argstring(args[i], i-1) for i in 1:length(args)), ' ')
    constr = join((constraint(args[i])     for i in 1:length(args)), ',')
    addr = string(".", sizeof(Int), "byte")

    asm = """
990:    nop
        .pushsection .note.stapsdt,"?","note"
        .balign 4
        .4byte 992f-991f, 994f-993f, 3
991:    .asciz "stapsdt"
992:    .balign 4
993:    $addr 990b
        $addr _.stapsdt.base
        $addr $(provider)_$(name)_semaphore
        .asciz "$provider"
        .asciz "$name"
        .asciz "$argstr"
994:    .balign 4
        .popsection
.ifndef _.stapsdt.base
        .pushsection .stapsdt.base,"aG","progbits",.stapsdt.base,comdat
        .weak _.stapsdt.base
        .hidden _.stapsdt.base
_.stapsdt.base: .space 1
        .size _.stapsdt.base, 1
        .popsection
.endif
"""

    ctx = LLVM.Interop.JuliaContext()
    mod = LLVM.Module("uprobe_$(provider)_$(name)", ctx)

    # Create semaphore variable
    int16_t = LLVM.Int16Type(ctx)
    semaphore = LLVM.GlobalVariable(mod, int16_t, "$(provider)_$(name)_semaphore")
    section!(semaphore, ".probes")
    linkage!(semaphore, LLVM.API.LLVMExternalLinkage)
    initializer!(semaphore, ConstantInt(int16_t, 0))

    # GDB is rather unhappy if there is no `.data` section
    gdb_unhappy = LLVM.GlobalVariable(mod, int16_t, "are_you_happy_now")
    section!(gdb_unhappy, ".data")
    linkage!(gdb_unhappy, LLVM.API.LLVMExternalLinkage)
    initializer!(gdb_unhappy, ConstantInt(int16_t, 0))

    # create function that will do a call to nop assembly
    rettyp = convert(LLVMType, Nothing)
    argtyp = LLVMType[convert.(LLVMType, args)...]

    ft = LLVM.FunctionType(rettyp, argtyp)
    f = LLVM.Function(mod, string("__uprobe_", provider, "_", name), ft)
    linkage!(f, LLVM.API.LLVMExternalLinkage)

    inline_asm = InlineAsm(ft, asm, constr, true)

    # generate IR
    Builder(ctx) do builder
        entry = BasicBlock(f, "entry", ctx)
        position!(builder, entry)

        val = call!(builder, inline_asm, collect(parameters(f)))
        ret!(builder)
    end

    triple = LLVM.triple()
    target = LLVM.Target(triple)
    objfile = tempname()
    TargetMachine(target, triple, "", "", LLVM.API.LLVMCodeGenLevelDefault, LLVM.API.LLVMRelocPIC) do tm
        LLVM.emit(tm, mod, LLVM.API.LLVMObjectFile, objfile)
    end

    run(`ld -shared $objfile -o $file`)
    return Libdl.dlopen(file, Libdl.RTLD_LOCAL)
end

end #module

ext: jl


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(csch(a+b*x), x) == :(-1*b^-1*arctanh(cosh(a+b*x)))
@test integrate(csch(a+b*x)^2, x) == :(-1*b^-1*coth(a+b*x))
@test integrate(csch(a+b*x)^3, x) == :((1/2)*b^-1*arctanh(cosh(a+b*x))+-1//2*b^-1*coth(a+b*x)*csch(a+b*x))
@test integrate(csch(a+b*x)^4, x) == :(b^-1*coth(a+b*x)+-1//3*b^-1*coth(a+b*x)^3)
@test integrate(csch(a+b*x)^5, x) == :(-3//8*b^-1*arctanh(cosh(a+b*x))+-1//4*b^-1*csch(a+b*x)^3*coth(a+b*x)+3//8*b^-1*coth(a+b*x)*csch(a+b*x))
@test integrate(csch(a+b*x)^6, x) == :(-1*b^-1*coth(a+b*x)+-1//5*b^-1*coth(a+b*x)^5+2//3*b^-1*coth(a+b*x)^3)
@test integrate(csch(a+b*x)^5//2, x) == :(-2//3*b^-1*csch(a+b*x)^3//2*cosh(a+b*x)+2//3*im*b^-1*(im*sinh(a+b*x))^(1/2)*csch(a+b*x)^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*a+(1/2)*im*b*x,2))
@test integrate(csch(a+b*x)^3//2, x) == :(-2*b^-1*csch(a+b*x)^(1/2)*cosh(a+b*x)+-2*im*b^-1*(im*sinh(a+b*x))^-1//2*csch(a+b*x)^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*a+(1/2)*im*b*x,2))
@test integrate(csch(a+b*x)^(1/2), x) == :(-2*im*b^-1*(im*sinh(a+b*x))^(1/2)*csch(a+b*x)^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*a+(1/2)*im*b*x,2))
@test integrate(csch(a+b*x)^-1//2, x) == :(-2*im*b^-1*(im*sinh(a+b*x))^-1//2*csch(a+b*x)^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*a+(1/2)*im*b*x,2))
@test integrate(csch(a+b*x)^-3//2, x) == :(2//3*b^-1*csch(a+b*x)^-1//2*cosh(a+b*x)+2//3*im*b^-1*(im*sinh(a+b*x))^(1/2)*csch(a+b*x)^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*a+(1/2)*im*b*x,2))
@test integrate(csch(a+b*x)^-5//2, x) == :(2//5*b^-1*csch(a+b*x)^-3//2*cosh(a+b*x)+6//5*im*b^-1*(im*sinh(a+b*x))^-1//2*csch(a+b*x)^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*a+(1/2)*im*b*x,2))
@test integrate((b*csch(c+d*x))^7//2, x) == :(-2//5*b*d^-1*(b*csch(c+d*x))^5//2*cosh(c+d*x)+6//5*b^3*d^-1*(b*csch(c+d*x))^(1/2)*cosh(c+d*x)+6//5*im*b^4*d^-1*(im*sinh(c+d*x))^-1//2*(b*csch(c+d*x))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^5//2, x) == :(-2//3*b*d^-1*(b*csch(c+d*x))^3//2*cosh(c+d*x)+2//3*im*b^2*d^-1*(im*sinh(c+d*x))^(1/2)*(b*csch(c+d*x))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^3//2, x) == :(-2*b*d^-1*(b*csch(c+d*x))^(1/2)*cosh(c+d*x)+-2*im*b^2*d^-1*(im*sinh(c+d*x))^-1//2*(b*csch(c+d*x))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^(1/2), x) == :(-2*im*d^-1*(im*sinh(c+d*x))^(1/2)*(b*csch(c+d*x))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^-1//2, x) == :(-2*im*d^-1*(im*sinh(c+d*x))^-1//2*(b*csch(c+d*x))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^-3//2, x) == :(2//3*b^-1*d^-1*(b*csch(c+d*x))^-1//2*cosh(c+d*x)+2//3*im*b^-2*d^-1*(im*sinh(c+d*x))^(1/2)*(b*csch(c+d*x))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^-5//2, x) == :(2//5*b^-1*d^-1*(b*csch(c+d*x))^-3//2*cosh(c+d*x)+6//5*im*b^-2*d^-1*(im*sinh(c+d*x))^-1//2*(b*csch(c+d*x))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((b*csch(c+d*x))^-7//2, x) == :(-10//21*b^-3*d^-1*(b*csch(c+d*x))^-1//2*cosh(c+d*x)+2//7*b^-1*d^-1*(b*csch(c+d*x))^-5//2*cosh(c+d*x)+-10//21*im*b^-4*d^-1*(im*sinh(c+d*x))^(1/2)*(b*csch(c+d*x))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*c+(1/2)*im*d*x,2))
@test integrate((-1*csch(x)^2)^5//2, x) == :(3//8*arcsin(coth(x))+1//4*(-1*csch(x)^2)^3//2*coth(x)+3//8*(-1*csch(x)^2)^(1/2)*coth(x))
@test integrate((-1*csch(x)^2)^3//2, x) == :((1/2)*arcsin(coth(x))+(1/2)*(-1*csch(x)^2)^(1/2)*coth(x))
@test integrate((-1*csch(x)^2)^(1/2), x) == :(arcsin(coth(x)))
@test integrate((-1*csch(x)^2)^-1//2, x) == :((-1*csch(x)^2)^-1//2*coth(x))
@test integrate((-1*csch(x)^2)^-3//2, x) == :(1//3*(-1*csch(x)^2)^-3//2*coth(x)+2//3*(-1*csch(x)^2)^-1//2*coth(x))
@test integrate((-1*csch(x)^2)^-5//2, x) == :(1//5*(-1*csch(x)^2)^-5//2*coth(x)+4//15*(-1*csch(x)^2)^-3//2*coth(x)+8//15*(-1*csch(x)^2)^-1//2*coth(x))
@test integrate((-1*csch(x)^2)^-7//2, x) == :(1//7*(-1*csch(x)^2)^-7//2*coth(x)+6//35*(-1*csch(x)^2)^-5//2*coth(x)+8//35*(-1*csch(x)^2)^-3//2*coth(x)+16//35*(-1*csch(x)^2)^-1//2*coth(x))
@test integrate((a*csch(x)^2)^5//2, x) == :(-3//8*a^5//2*arctanh(a^(1/2)*(a*csch(x)^2)^-1//2*coth(x))+-1//4*a*(a*csch(x)^2)^3//2*coth(x)+3//8*a^2*(a*csch(x)^2)^(1/2)*coth(x))
@test integrate((a*csch(x)^2)^3//2, x) == :((1/2)*a^3//2*arctanh(a^(1/2)*(a*csch(x)^2)^-1//2*coth(x))+-1//2*a*(a*csch(x)^2)^(1/2)*coth(x))
@test integrate((a*csch(x)^2)^(1/2), x) == :(-1*a^(1/2)*arctanh(a^(1/2)*(a*csch(x)^2)^-1//2*coth(x)))
@test integrate((a*csch(x)^2)^-1//2, x) == :((a*csch(x)^2)^-1//2*coth(x))
@test integrate((a*csch(x)^2)^-3//2, x) == :(1//3*(a*csch(x)^2)^-3//2*coth(x)+-2//3*a^-1*(a*csch(x)^2)^-1//2*coth(x))
@test integrate((a*csch(x)^2)^-5//2, x) == :(1//5*(a*csch(x)^2)^-5//2*coth(x)+-4//15*a^-1*(a*csch(x)^2)^-3//2*coth(x)+8//15*a^-2*(a*csch(x)^2)^-1//2*coth(x))
@test integrate((a*csch(x)^2)^-7//2, x) == :(1//7*(a*csch(x)^2)^-7//2*coth(x)+-16//35*a^-3*(a*csch(x)^2)^-1//2*coth(x)+-6//35*a^-1*(a*csch(x)^2)^-5//2*coth(x)+8//35*a^-2*(a*csch(x)^2)^-3//2*coth(x))
@test integrate((a*csch(x)^3)^5//2, x) == :(-154//585*a^2*(a*csch(x)^3)^(1/2)*coth(x)+-2//13*a^2*(a*csch(x)^3)^(1/2)*csch(x)^4*coth(x)+22//117*a^2*(a*csch(x)^3)^(1/2)*csch(x)^2*coth(x)+154//195*a^2*(a*csch(x)^3)^(1/2)*cosh(x)*sinh(x)+-154//195*im*a^2*(im*sinh(x))^-1//2*(a*csch(x)^3)^(1/2)*sinh(x)^2*Elliptic.E(1//4*pi+-1//2*im*x,2))
@test integrate((a*csch(x)^3)^3//2, x) == :(10//21*a*(a*csch(x)^3)^(1/2)*cosh(x)+-2//7*a*(a*csch(x)^3)^(1/2)*coth(x)*csch(x)+10//21*im*a*(im*sinh(x))^(1/2)*(a*csch(x)^3)^(1/2)*Elliptic.F(1//4*pi+-1//2*im*x,2)*sinh(x))
@test integrate((a*csch(x)^3)^-1//2, x) == :(2//3*(a*csch(x)^3)^-1//2*coth(x)+-2//3*im*(im*sinh(x))^(1/2)*(a*csch(x)^3)^-1//2*csch(x)^2*Elliptic.F(1//4*pi+-1//2*im*x,2))
@test integrate((a*csch(x)^3)^-3//2, x) == :(-14//45*a^-1*(a*csch(x)^3)^-1//2*cosh(x)+2//9*a^-1*(a*csch(x)^3)^-1//2*sinh(x)^2*cosh(x)+14//15*im*a^-1*(im*sinh(x))^-1//2*(a*csch(x)^3)^-1//2*Elliptic.E(1//4*pi+-1//2*im*x,2)*csch(x))
@test integrate((a*csch(x)^3)^-5//2, x) == :(-26//77*a^-2*(a*csch(x)^3)^-1//2*coth(x)+-26//165*a^-2*(a*csch(x)^3)^-1//2*sinh(x)^3*cosh(x)+2//15*a^-2*(a*csch(x)^3)^-1//2*sinh(x)^5*cosh(x)+78//385*a^-2*(a*csch(x)^3)^-1//2*cosh(x)*sinh(x)+26//77*im*a^-2*(im*sinh(x))^(1/2)*(a*csch(x)^3)^-1//2*csch(x)^2*Elliptic.F(1//4*pi+-1//2*im*x,2))
@test integrate((a*csch(x)^4)^7//2, x) == :(-1*a^3*(a*csch(x)^4)^(1/2)*cosh(x)*sinh(x)+-3*a^3*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^3+2*a^3*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)+-5//3*a^3*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^7+-1//13*a^3*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^11+6//11*a^3*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^9+20//7*a^3*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^5)
@test integrate((a*csch(x)^4)^5//2, x) == :(-1*a^2*(a*csch(x)^4)^(1/2)*cosh(x)*sinh(x)+-6//5*a^2*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^3+-1//9*a^2*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^7+4//3*a^2*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)+4//7*a^2*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^5)
@test integrate((a*csch(x)^4)^3//2, x) == :(-1*a*(a*csch(x)^4)^(1/2)*cosh(x)*sinh(x)+-1//5*a*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x)^3+2//3*a*(a*csch(x)^4)^(1/2)*cosh(x)^2*coth(x))
@test integrate((a*csch(x)^4)^(1/2), x) == :(-1*(a*csch(x)^4)^(1/2)*cosh(x)*sinh(x))
@test integrate((a*csch(x)^4)^-1//2, x) == :((1/2)*(a*csch(x)^4)^-1//2*coth(x)+-1//2*x*(a*csch(x)^4)^-1//2*csch(x)^2)
@test integrate((a*csch(x)^4)^-3//2, x) == :(5//16*a^-1*(a*csch(x)^4)^-1//2*coth(x)+-5//16*x*a^-1*(a*csch(x)^4)^-1//2*csch(x)^2+-5//24*a^-1*(a*csch(x)^4)^-1//2*cosh(x)*sinh(x)+1//6*a^-1*(a*csch(x)^4)^-1//2*sinh(x)^3*cosh(x))
@test integrate((a*csch(x)^4)^-5//2, x) == :(63//256*a^-2*(a*csch(x)^4)^-1//2*coth(x)+-63//256*x*a^-2*(a*csch(x)^4)^-1//2*csch(x)^2+-21//128*a^-2*(a*csch(x)^4)^-1//2*cosh(x)*sinh(x)+-9//80*a^-2*(a*csch(x)^4)^-1//2*sinh(x)^5*cosh(x)+1//10*a^-2*(a*csch(x)^4)^-1//2*sinh(x)^7*cosh(x)+21//160*a^-2*(a*csch(x)^4)^-1//2*sinh(x)^3*cosh(x))
@test integrate((a+im*a*csch(a+b*x))^-1, x) == :(x*a^-1+-1*b^-1*(a+im*a*csch(a+b*x))^-1*coth(a+b*x))
@test integrate((a+-1*im*a*csch(a+b*x))^-1, x) == :(x*a^-1+-1*b^-1*(a+-1*im*a*csch(a+b*x))^-1*coth(a+b*x))
@test integrate((a+im*a*csch(c+d*x))^5//2, x) == :(2*a^5//2*d^-1*arctanh(a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x))+2//3*a^2*d^-1*(a+im*a*csch(c+d*x))^(1/2)*coth(c+d*x)+14//3*a^3*d^-1*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x))
@test integrate((a+im*a*csch(c+d*x))^3//2, x) == :(2*a^3//2*d^-1*arctanh(a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x))+2*a^2*d^-1*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x))
@test integrate((a+im*a*csch(c+d*x))^(1/2), x) == :(2*a^(1/2)*d^-1*arctanh(a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x)))
@test integrate((a+im*a*csch(c+d*x))^-1//2, x) == :(2*a^-1//2*d^-1*arctanh(a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x))+-1*2^(1/2)*a^-1//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x)))
@test integrate((a+im*a*csch(c+d*x))^-3//2, x) == :(2*a^-3//2*d^-1*arctanh(a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x))+-1//2*d^-1*(a+im*a*csch(c+d*x))^-3//2*coth(c+d*x)+-5//4*2^(1/2)*a^-3//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+im*a*csch(c+d*x))^-1//2*coth(c+d*x)))
@test integrate((a+-1*im*a*csch(c+d*x))^(1/2), x) == :(2*a^(1/2)*d^-1*arctanh(a^(1/2)*(a+-1*im*a*csch(c+d*x))^-1//2*coth(c+d*x)))
@test integrate((a+-1*im*a*csch(c+d*x))^-1//2, x) == :(2*a^-1//2*d^-1*arctanh(a^(1/2)*(a+-1*im*a*csch(c+d*x))^-1//2*coth(c+d*x))+-1*2^(1/2)*a^-1//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+-1*im*a*csch(c+d*x))^-1//2*coth(c+d*x)))
@test integrate((3+3*im*csch(x))^(1/2), x) == :(2*3^(1/2)*arctanh((1+im*csch(x))^-1//2*coth(x)))
@test integrate((3+-3*im*csch(x))^(1/2), x) == :(2*3^(1/2)*arctanh((1+-1*im*csch(x))^-1//2*coth(x)))
@test integrate((-3+3*im*csch(x))^(1/2), x) == :(-2*3^(1/2)*arctan((-1+im*csch(x))^-1//2*coth(x)))
@test integrate((-3+-3*im*csch(x))^(1/2), x) == :(-2*3^(1/2)*arctan((-1+-1*im*csch(x))^-1//2*coth(x)))
@test integrate((im+csch(x))^-1*sinh(x)^4, x) == :(-4*cosh(x)+4//3*cosh(x)^3+-15//8*im*x+-1*(im+csch(x))^-1*sinh(x)^3*cosh(x)+-5//4*im*sinh(x)^3*cosh(x)+15//8*im*cosh(x)*sinh(x))
@test integrate((im+csch(x))^-1*sinh(x)^3, x) == :(-3//2*x+4*im*cosh(x)+-4//3*im*cosh(x)^3+3//2*cosh(x)*sinh(x)+-1*(im+csch(x))^-1*sinh(x)^2*cosh(x))
@test integrate((im+csch(x))^-1*sinh(x)^2, x) == :(2*cosh(x)+3//2*im*x+-1*(im+csch(x))^-1*cosh(x)*sinh(x)+-3//2*im*cosh(x)*sinh(x))
@test integrate((im+csch(x))^-1*sinh(x), x) == :(x+-1*(im+csch(x))^-1*cosh(x)+-2*im*cosh(x))
@test integrate((im+csch(x))^-1*csch(x), x) == :(im*(im+csch(x))^-1*coth(x))
@test integrate((im+csch(x))^-1*csch(x)^2, x) == :(-1*arctanh(cosh(x))+(im+csch(x))^-1*coth(x))
@test integrate((im+csch(x))^-1*csch(x)^3, x) == :(-1*coth(x)+im*arctanh(cosh(x))+-1*im*(im+csch(x))^-1*coth(x))
@test integrate((im+csch(x))^-1*csch(x)^4, x) == :(3//2*arctanh(cosh(x))+2*im*coth(x)+-3//2*coth(x)*csch(x)+(im+csch(x))^-1*csch(x)^2*coth(x))
@test integrate((a+b*csch(c+d*x))^4, x) == :(x*a^4+-1//3*b^2*d^-1*(a+b*csch(c+d*x))^2*coth(c+d*x)+-1//3*b^2*d^-1*(-2*b^2+17*a^2)*coth(c+d*x)+-2*a*b*d^-1*(-1*b^2+2*a^2)*arctanh(cosh(c+d*x))+-4//3*a*b^3*d^-1*coth(c+d*x)*csch(c+d*x))
@test integrate((a+b*csch(c+d*x))^3, x) == :(x*a^3+-5//2*a*b^2*d^-1*coth(c+d*x)+-1//2*b*d^-1*(-1*b^2+6*a^2)*arctanh(cosh(c+d*x))+-1//2*b^2*d^-1*(a+b*csch(c+d*x))*coth(c+d*x))
@test integrate((a+b*csch(c+d*x))^2, x) == :(x*a^2+-1*b^2*d^-1*coth(c+d*x)+-2*a*b*d^-1*arctanh(cosh(c+d*x)))
@test integrate(a+b*csch(c+d*x), x) == :(a*x+-1*b*d^-1*arctanh(cosh(c+d*x)))
@test integrate((a+b*csch(c+d*x))^-1, x) == :(x*a^-1+2*b*a^-1*d^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*csch(c+d*x))^-2, x) == :(x*a^-2+-1*a^-1*b^2*d^-1*(a+b*csch(c+d*x))^-1*(a^2+b^2)^-1*coth(c+d*x)+2*b*a^-2*d^-1*(a^2+b^2)^-3//2*(b^2+2*a^2)*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*csch(c+d*x))^-3, x) == :(x*a^-3+b*a^-3*d^-1*(a^2+b^2)^-5//2*(2*b^4+6*a^4+5*a^2*b^2)*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*c+(1/2)*d*x)))+-1//2*a^-1*b^2*d^-1*(a+b*csch(c+d*x))^-2*(a^2+b^2)^-1*coth(c+d*x)+-1//2*a^-2*b^2*d^-1*(a+b*csch(c+d*x))^-1*(a^2+b^2)^-2*(2*b^2+5*a^2)*coth(c+d*x))
@test integrate((a+b*csch(x))^-1*sinh(x)^3, x) == :(-1//3*a^-3*(-3*b^2+2*a^2)*cosh(x)+1//3*a^-1*sinh(x)^2*cosh(x)+(1/2)*b*x*a^-4*(a^2+-2*b^2)+-2*a^-4*b^4*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x)))+-1//2*b*a^-2*cosh(x)*sinh(x))
@test integrate((a+b*csch(x))^-1*sinh(x)^2, x) == :((1/2)*a^-1*cosh(x)*sinh(x)+-1*b*a^-2*cosh(x)+-1//2*x*a^-3*(a^2+-2*b^2)+2*a^-3*b^3*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*sinh(x), x) == :(a^-1*cosh(x)+-1*b*x*a^-2+-2*a^-2*b^2*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*csch(x), x) == :(-2*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*csch(x)^2, x) == :(-1*b^-1*arctanh(cosh(x))+2*a*b^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*csch(x)^3, x) == :(-1*b^-1*coth(x)+a*b^-2*arctanh(cosh(x))+-2*a^2*b^-2*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*csch(x)^4, x) == :(a*b^-2*coth(x)+-1//2*b^-1*coth(x)*csch(x)+-1//2*b^-3*(-1*b^2+2*a^2)*arctanh(cosh(x))+2*a^3*b^-3*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((im+csch(x))^-1*cosh(x)^4, x) == :(1//3*cosh(x)^3+1//8*im*x+-1//4*im*cosh(x)^3*sinh(x)+1//8*im*cosh(x)*sinh(x))
@test integrate((im+csch(x))^-1*cosh(x)^3, x) == :((1/2)*sinh(x)^2+-1//3*im*sinh(x)^3)
@test integrate((im+csch(x))^-1*cosh(x)^2, x) == :((1/2)*im*x+-1//2*im*cosh(x)*sinh(x)+cosh(x))
@test integrate((im+csch(x))^-1*cosh(x), x) == :(-1*im*sinh(x)+log(im+-1*sinh(x)))
@test integrate((im+csch(x))^-1*sech(x), x) == :(-1//2*sech(x)^2+-1//2*im*arctan(sinh(x))+(1/2)*im*sech(x)*tanh(x))
@test integrate((im+csch(x))^-1*sech(x)^2, x) == :(-1//3*sech(x)^3+-1//3*im*tanh(x)^3)
@test integrate((im+csch(x))^-1*sech(x)^3, x) == :(-1//4*sech(x)^4+-1//8*im*arctan(sinh(x))+-1//8*im*sech(x)*tanh(x)+1//4*im*sech(x)^3*tanh(x))
@test integrate((im+csch(x))^-1*sech(x)^4, x) == :(-1//5*sech(x)^5+-1//3*im*tanh(x)^3+1//5*im*tanh(x)^5)
@test integrate((a+b*csch(x))^-1*cosh(x)^5, x) == :(1//5*a^-1*sinh(x)^5+a^-5*(a^2+b^2)^2*sinh(x)+-1//4*b*a^-2*sinh(x)^4+1//3*a^-3*sinh(x)^3*(b^2+2*a^2)+-1*b*a^-6*(a^2+b^2)^2*log(b+a*sinh(x))+-1//2*b*a^-4*sinh(x)^2*(b^2+2*a^2))
@test integrate((a+b*csch(x))^-1*cosh(x)^4, x) == :(-1//8*a^-4*(8*b*(a^2+b^2)+-1*a*(3*a^2+4*b^2)*sinh(x))*cosh(x)+-1//12*a^-2*cosh(x)^3*(4b+-3*a*sinh(x))+1//8*x*a^-5*(3*a^4+8*b^4+12*a^2*b^2)+2*b*a^-5*(a^2+b^2)^3//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*cosh(x)^3, x) == :(1//3*a^-1*sinh(x)^3+a^-3*(a^2+b^2)*sinh(x)+-1//2*b*a^-2*sinh(x)^2+-1*b*a^-4*(a^2+b^2)*log(b+a*sinh(x)))
@test integrate((a+b*csch(x))^-1*cosh(x)^2, x) == :((1/2)*x*a^-3*(a^2+2*b^2)+-1//2*a^-2*(2b+-1*a*sinh(x))*cosh(x)+2*b*a^-3*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*cosh(x), x) == :(a^-1*sinh(x)+-1*b*a^-2*log(b+a*sinh(x)))
@test integrate((a+b*csch(x))^-1*sech(x), x) == :((2b+2*im*a)^-1*log(im+-1*sinh(x))+-1*(-2b+2*im*a)^-1*log(im+sinh(x))+-1*b*(a^2+b^2)^-1*log(b+a*sinh(x)))
@test integrate((a+b*csch(x))^-1*sech(x)^2, x) == :(-1*(a^2+b^2)^-1*(b+-1*a*sinh(x))*sech(x)+2*a*b*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*sech(x)^3, x) == :(-1*(2*a^2+2*b^2)^-1*sech(x)^2*(b+-1*a*sinh(x))+-1*b*a^2*(a^2+b^2)^-2*log(b+a*sinh(x))+-1//4*im*a*(a+-1*im*b)^-2*log(im+-1*sinh(x))+1//4*im*a*(a+im*b)^-2*log(im+sinh(x)))
@test integrate((a+b*csch(x))^-1*sech(x)^4, x) == :(-1*(3*a^2+3*b^2)^-1*sech(x)^3*(b+-1*a*sinh(x))+-1//3*(a^2+b^2)^-2*(3*b*a^2+-1*a*(-1*b^2+2*a^2)*sinh(x))*sech(x)+2*b*a^3*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*sech(x)^5, x) == :(-1*(4*a^2+4*b^2)^-1*sech(x)^4*(b+-1*a*sinh(x))+-1//8*(a^2+b^2)^-2*sech(x)^2*(4*b*a^2+-1*a*(-1*b^2+3*a^2)*sinh(x))+-1*b*a^4*(a^2+b^2)^-3*log(b+a*sinh(x))+-1//16*a*(a+-1*im*b)^-3*(b+3*im*a)*log(im+-1*sinh(x))+1//16*a*(-1b+im*a)^-3*(3a+im*b)*log(im+sinh(x)))
@test integrate((im+csch(x))^-1*tanh(x)^5, x) == :(-1*im*(4+-4*im*sinh(x))^-1+-15*im*(16+16*im*sinh(x))^-1+-21//32*im*log(im+-1*sinh(x))+-11//32*im*log(im+sinh(x))+-1//24*im*(1+im*sinh(x))^-3+1//32*im*(1+-1*im*sinh(x))^-2+9//32*im*(1+im*sinh(x))^-2)
@test integrate((im+csch(x))^-1*tanh(x)^4, x) == :(tanh(x)^3*(-4//15*csch(x)+1//3*im)+tanh(x)^5*(-1//5*csch(x)+1//5*im)+(im+-8//15*csch(x))*tanh(x)+-1*im*x)
@test integrate((im+csch(x))^-1*tanh(x)^3, x) == :(-1*im*(8+-8*im*sinh(x))^-1+-3*im*(4+4*im*sinh(x))^-1+-11//16*im*log(im+-1*sinh(x))+-5//16*im*log(im+sinh(x))+1//8*im*(1+im*sinh(x))^-2)
@test integrate((im+csch(x))^-1*tanh(x)^2, x) == :(tanh(x)^3*(-1//3*csch(x)+1//3*im)+(im+-2//3*csch(x))*tanh(x)+-1*im*x)
@test integrate((im+csch(x))^-1*tanh(x), x) == :(-1*im*(2+2*im*sinh(x))^-1+-3//4*im*log(im+-1*sinh(x))+-1//4*im*log(im+sinh(x)))
@test integrate((im+csch(x))^-1*coth(x), x) == :(-1*im*log(im+-1*sinh(x)))
@test integrate((im+csch(x))^-1*coth(x)^2, x) == :(-1*arctanh(cosh(x))+-1*im*x)
@test integrate((im+csch(x))^-1*coth(x)^3, x) == :(-1*csch(x)+-1*im*log(sinh(x)))
@test integrate((im+csch(x))^-1*coth(x)^4, x) == :(-1//2*arctanh(cosh(x))+(1/2)*(-1*csch(x)+2im)*coth(x)+-1*im*x)
@test integrate((im+csch(x))^-1*coth(x)^5, x) == :(-1*csch(x)+-1//3*csch(x)^3+(1/2)*im*csch(x)^2+-1*im*log(sinh(x)))
@test integrate((im+csch(x))^-1*coth(x)^6, x) == :(-3//8*arctanh(cosh(x))+-1*im*x+1//8*(-3*csch(x)+8im)*coth(x)+1//12*coth(x)^3*(-3*csch(x)+4im))
@test integrate((a+b*csch(x))^-1*tanh(x)^5, x) == :(a^-1*log(sinh(x))+-1*b^5*(a^2+b^2)^-3*arctan(sinh(x))+-1*(4*a^2+4*b^2)^-1*tanh(x)^4*(a+-1*b*csch(x))+-3*b*(8*a^2+8*b^2)^-1*arctan(sinh(x))+-1//2*b^3*(a^2+b^2)^-2*arctan(sinh(x))+-1//2*(a^2+b^2)^-2*tanh(x)^2*(a*(a^2+2*b^2)+-1*b^3*csch(x))+a^-1*b^6*(a^2+b^2)^-3*log(a+b*csch(x))+-1*a*(a^2+b^2)^-3*(a^4+3*b^4+3*a^2*b^2)*log(tanh(x))+3*b*(8*a^2+8*b^2)^-1*sech(x)*tanh(x))
@test integrate((a+b*csch(x))^-1*tanh(x)^4, x) == :(a*x*(a^2+b^2)^-1+b*(a^2+b^2)^-1*sech(x)+b^3*(a^2+b^2)^-2*sech(x)+-1*a*(a^2+b^2)^-1*tanh(x)+-1*a*(3*a^2+3*b^2)^-1*tanh(x)^3+-1*b*(3*a^2+3*b^2)^-1*sech(x)^3+a*x*b^2*(a^2+b^2)^-2+x*a^-1*b^4*(a^2+b^2)^-2+-1*a*b^2*(a^2+b^2)^-2*tanh(x)+2*a^-1*b^5*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*tanh(x)^3, x) == :(a^-1*log(sinh(x))+-1*b*(2*a^2+2*b^2)^-1*arctan(sinh(x))+-1*b^3*(a^2+b^2)^-2*arctan(sinh(x))+-1*(2*a^2+2*b^2)^-1*tanh(x)^2*(a+-1*b*csch(x))+a^-1*b^4*(a^2+b^2)^-2*log(a+b*csch(x))+-1*a*(a^2+b^2)^-2*(a^2+2*b^2)*log(tanh(x)))
@test integrate((a+b*csch(x))^-1*tanh(x)^2, x) == :(a*x*(a^2+b^2)^-1+b*(a^2+b^2)^-1*sech(x)+-1*a*(a^2+b^2)^-1*tanh(x)+x*a^-1*b^2*(a^2+b^2)^-1+2*a^-1*b^3*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*tanh(x), x) == :(a^-1*log(sinh(x))+-1*a*(a^2+b^2)^-1*log(tanh(x))+-1*b*(a^2+b^2)^-1*arctan(sinh(x))+a^-1*b^2*(a^2+b^2)^-1*log(a+b*csch(x)))
@test integrate((a+b*csch(x))^-1*coth(x), x) == :(a^-1*log(a+b*csch(x))+a^-1*log(sinh(x)))
@test integrate((a+b*csch(x))^-1*coth(x)^2, x) == :(x*a^-1+-1*b^-1*arctanh(cosh(x))+2*a^-1*b^-1*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*coth(x)^3, x) == :(a^-1*log(sinh(x))+(a^-1+a*b^-2)*log(a+b*csch(x))+-1*b^-1*csch(x))
@test integrate((a+b*csch(x))^-1*coth(x)^4, x) == :(x*a^-1+a*b^-2*coth(x)+-1//2*b^-1*coth(x)*csch(x)+-1//2*b^-3*(2*a^2+3*b^2)*arctanh(cosh(x))+2*a^-1*b^-3*(a^2+b^2)^3//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x))))
@test integrate((a+b*csch(x))^-1*coth(x)^5, x) == :(a^-1*log(sinh(x))+-1//3*b^-1*csch(x)^3+(1/2)*a*b^-2*csch(x)^2+-1*b^-3*(a^2+2*b^2)*csch(x)+a^-1*b^-4*(a^2+b^2)^2*log(a+b*csch(x)))
@test integrate((a+b*csch(x))^-1*coth(x)^6, x) == :(x*a^-1+-3//8*b^-1*arctanh(cosh(x))+(1/2)*b^-3*(a^2+3*b^2)*arctanh(cosh(x))+-1*a*b^-2*coth(x)+-1*b^-5*(a^4+3*b^4+3*a^2*b^2)*arctanh(cosh(x))+-1//4*b^-1*csch(x)^3*coth(x)+1//3*a*b^-2*coth(x)^3+3//8*b^-1*coth(x)*csch(x)+a*b^-4*(a^2+3*b^2)*coth(x)+2*a^-1*b^-5*(a^2+b^2)^5//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*tanh((1/2)*x)))+-1//2*b^-3*(a^2+3*b^2)*coth(x)*csch(x))
@test integrate((a+b*csch(x))^-1*coth(x)^7, x) == :(a^-1*log(sinh(x))+-1//5*b^-1*csch(x)^5+-1*b^-5*(a^4+3*b^4+3*a^2*b^2)*csch(x)+-1//3*b^-3*csch(x)^3*(a^2+3*b^2)+1//4*a*b^-2*csch(x)^4+a^-1*b^-6*(a^2+b^2)^3*log(a+b*csch(x))+(1/2)*a*b^-4*csch(x)^2*(a^2+3*b^2))
@test integrate((csch(a*c+b*c*x)^2)^7//2*exp(c*(a+b*x)), x) == :(-48*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-4*sinh(a*c+b*c*x)+-32//3*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-6*sinh(a*c+b*c*x)+64//3*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-3*sinh(a*c+b*c*x)+192//5*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-5*sinh(a*c+b*c*x))
@test integrate((csch(a*c+b*c*x)^2)^5//2*exp(c*(a+b*x)), x) == :(-8*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-2*sinh(a*c+b*c*x)+-4*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-4*sinh(a*c+b*c*x)+32//3*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-3*sinh(a*c+b*c*x))
@test integrate((csch(a*c+b*c*x)^2)^3//2*exp(c*(a+b*x)), x) == :(-2*b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-2*exp(4*c*(a+b*x))*sinh(a*c+b*c*x))
@test integrate((csch(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(csch(a*c+b*c*x)^2)^(1/2)*log(1+-1*exp(2*c*(a+b*x)))*sinh(a*c+b*c*x))
@test integrate((csch(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x)), x) == :(-1//2*x*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)+1//4*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(2*c*(a+b*x)))
@test integrate((csch(a*c+b*c*x)^2)^-3//2*exp(c*(a+b*x)), x) == :(3//8*x*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)+-3//16*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(2*c*(a+b*x))+1//16*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(-2*c*(a+b*x))+1//32*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(4*c*(a+b*x)))
@test integrate((csch(a*c+b*c*x)^2)^-5//2*exp(c*(a+b*x)), x) == :(-5//16*x*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)+-5//64*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(-2*c*(a+b*x))+-5//128*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(4*c*(a+b*x))+1//128*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(-4*c*(a+b*x))+1//192*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(6*c*(a+b*x))+5//32*b^-1*c^-1*(csch(a*c+b*c*x)^2)^-1//2*csch(a*c+b*c*x)*exp(2*c*(a+b*x)))
@test integrate(x^5*csch(2*log(c*x))^-1//2, x) == :(1//7*x^6*csch(2*log(c*x))^-1//2+-2//21*c^-4*x^2*csch(2*log(c*x))^-1//2+2//21*c^-7*x^-1*(1+-1*c^-4*x^-4)^-1//2*csch(2*log(c*x))^-1//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(x^4*csch(2*log(c*x))^-1//2, x) == :(1//6*c^-4*x^5*csch(2*log(c*x))^-1//2*(c^4+-1*x^-4))
@test integrate(x^3*csch(2*log(c*x))^-1//2, x) == :(-2//5*c^-4*csch(2*log(c*x))^-1//2+1//5*x^4*csch(2*log(c*x))^-1//2+-2//5*c^-5*x^-1*(1+-1*c^-4*x^-4)^-1//2*csch(2*log(c*x))^-1//2*Elliptic.E(arccsc(c*x),-1)+2//5*c^-5*x^-1*(1+-1*c^-4*x^-4)^-1//2*csch(2*log(c*x))^-1//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(x^2*csch(2*log(c*x))^-1//2, x) == :(1//4*x^3*csch(2*log(c*x))^-1//2+-1//4*c^-4*x^-1*(1+-1*c^-4*x^-4)^-1//2*csch(2*log(c*x))^-1//2*arctanh((1+-1*c^-4*x^-4)^(1/2)))
@test integrate(x*csch(2*log(c*x))^-1//2, x) == :(1//3*x^2*csch(2*log(c*x))^-1//2+2//3*c^-3*x^-1*(1+-1*c^-4*x^-4)^-1//2*csch(2*log(c*x))^-1//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(csch(2*log(c*x))^-1//2, x) == :((1/2)*x*csch(2*log(c*x))^-1//2+(1/2)*c^-2*x^-1*(1+-1*c^-4*x^-4)^-1//2*csch(2*log(c*x))^-1//2*arccsc(c^2*x^2))
@test integrate(x^-1*csch(2*log(c*x))^(1/2), x) == :(im*(im*sinh(2*log(c*x)))^(1/2)*csch(2*log(c*x))^(1/2)*Elliptic.F(1//4*pi+-1*im*log(c*x),2))
@test integrate(x^-2*csch(2*log(c*x))^(1/2), x) == :(-1//2*x*c^2*(1+-1*c^-4*x^-4)^(1/2)*csch(2*log(c*x))^(1/2)*arccsc(c^2*x^2))
@test integrate(x^-3*csch(2*log(c*x))^(1/2), x) == :(x*c^3*(1+-1*c^-4*x^-4)^(1/2)*csch(2*log(c*x))^(1/2)*Elliptic.F(arccsc(c*x),-1)+-1*x*c^3*(1+-1*c^-4*x^-4)^(1/2)*csch(2*log(c*x))^(1/2)*Elliptic.E(arccsc(c*x),-1))
@test integrate(x^-4*csch(2*log(c*x))^(1/2), x) == :(x*csch(2*log(c*x))^(1/2)*((1/2)*c^4+-1//2*x^-4))
@test integrate(x^-5*csch(2*log(c*x))^(1/2), x) == :(csch(2*log(c*x))^(1/2)*(-1//3*x^-4+1//3*c^4)+-1//3*x*c^5*(1+-1*c^-4*x^-4)^(1/2)*csch(2*log(c*x))^(1/2)*Elliptic.F(arccsc(c*x),-1))
@test integrate(x^8*csch(2*log(c*x))^-3//2, x) == :(1//12*x^9*csch(2*log(c*x))^-3//2+-1*x^5*(-16*x^-4+16*c^4)^-1*csch(2*log(c*x))^-3//2+1//32*x*c^-4*(c^4+-1*x^-4)^-1*csch(2*log(c*x))^-3//2+1//32*c^-12*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*arctanh((1+-1*c^-4*x^-4)^(1/2)))
@test integrate(x^7*csch(2*log(c*x))^-3//2, x) == :(1//11*x^8*csch(2*log(c*x))^-3//2+-6*x^4*(-77*x^-4+77*c^4)^-1*csch(2*log(c*x))^-3//2+4//77*c^-4*(c^4+-1*x^-4)^-1*csch(2*log(c*x))^-3//2+-4//77*c^-11*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(x^6*csch(2*log(c*x))^-3//2, x) == :(1//10*c^-4*x^7*csch(2*log(c*x))^-3//2*(c^4+-1*x^-4))
@test integrate(x^5*csch(2*log(c*x))^-3//2, x) == :(1//9*x^6*csch(2*log(c*x))^-3//2+-2*x^2*(-15*x^-4+15*c^4)^-1*csch(2*log(c*x))^-3//2+4//15*c^-4*x^-2*(c^4+-1*x^-4)^-1*csch(2*log(c*x))^-3//2+-4//15*c^-9*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*Elliptic.F(arccsc(c*x),-1)+4//15*c^-9*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*Elliptic.E(arccsc(c*x),-1))
@test integrate(x^4*csch(2*log(c*x))^-3//2, x) == :(1//8*x^5*csch(2*log(c*x))^-3//2+-3*x*(-16*x^-4+16*c^4)^-1*csch(2*log(c*x))^-3//2+3//16*c^-8*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*arctanh((1+-1*c^-4*x^-4)^(1/2)))
@test integrate(x^3*csch(2*log(c*x))^-3//2, x) == :(-2*(-7*x^-4+7*c^4)^-1*csch(2*log(c*x))^-3//2+1//7*x^4*csch(2*log(c*x))^-3//2+-4//7*c^-7*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(x^2*csch(2*log(c*x))^-3//2, x) == :(1//6*x^3*csch(2*log(c*x))^-3//2+-1*x^-1*(-2*x^-4+2*c^4)^-1*csch(2*log(c*x))^-3//2+-1//2*c^-6*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*arccsc(c^2*x^2))
@test integrate(x*csch(2*log(c*x))^-3//2, x) == :(1//5*x^2*csch(2*log(c*x))^-3//2+-6*x^-2*(-5*x^-4+5*c^4)^-1*csch(2*log(c*x))^-3//2+-12//5*c^-5*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*Elliptic.E(arccsc(c*x),-1)+12//5*c^-5*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(csch(2*log(c*x))^-3//2, x) == :(1//4*x*csch(2*log(c*x))^-3//2+3*x^-3*(-4*x^-4+4*c^4)^-1*csch(2*log(c*x))^-3//2+-3//4*c^-4*x^-3*(1+-1*c^-4*x^-4)^-3//2*csch(2*log(c*x))^-3//2*arctanh((1+-1*c^-4*x^-4)^(1/2)))
@test integrate(x^-1*csch(2*log(c*x))^3//2, x) == :(-1*csch(2*log(c*x))^(1/2)*cosh(2*log(c*x))+im*(im*sinh(2*log(c*x)))^-1//2*csch(2*log(c*x))^-1//2*Elliptic.E(1//4*pi+-1*im*log(c*x),2))
@test integrate(x^-2*csch(2*log(c*x))^3//2, x) == :(x^3*csch(2*log(c*x))^3//2*((1/2)*x^-4+-1//2*c^4))
@test integrate(x^-3*csch(2*log(c*x))^3//2, x) == :(x^2*csch(2*log(c*x))^3//2*((1/2)*x^-4+-1//2*c^4)+(1/2)*c^5*x^3*(1+-1*c^-4*x^-4)^3//2*csch(2*log(c*x))^3//2*Elliptic.F(arccsc(c*x),-1))
@test integrate(x^-4*csch(2*log(c*x))^3//2, x) == :(x*csch(2*log(c*x))^3//2*((1/2)*x^-4+-1//2*c^4)+(1/2)*c^6*x^3*(1+-1*c^-4*x^-4)^3//2*csch(2*log(c*x))^3//2*arccsc(c^2*x^2))
@test integrate(-1*(1+-1*b^2*n^2)*csch(a+b*log(c*x^n))+2*b^2*n^2*csch(a+b*log(c*x^n))^3, x) == :(-1*x*csch(a+b*log(c*x^n))+-1*b*n*x*coth(a+b*log(c*x^n))*csch(a+b*log(c*x^n)))
@test integrate(csch(a+2*log(c*x^(1/2)))^3, x) == :(-2*c^6*(c^4+-1*x^-2*exp(-2a))^-2*exp(-1a))
@test integrate(csch(a+2*log(c*x^-1//2))^3, x) == :(2*c^2*(-1*c^4*x^-2+exp(-2a))^-2*exp(-3a))
@test integrate(csch(a+n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(-1*x*(c*x^n)^(-2*n^-1*(2+-1p)^-1)*(2+-2p)^-1*csch(a+-1*n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+-1*(c*x^n)^(2*n^-1*(2+-1p)^-1)*exp(-2a))*(2+-1p)*exp(2a))
@test integrate(csch(a+-1*n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(x*(2+-2p)^-1*csch(a+n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+-1*(c*x^n)^(-2*n^-1*(2+-1p)^-1)*exp(-2a))*(2+-1p))
@test integrate(x^-1*csch(a+b*log(c*x^n)), x) == :(-1*b^-1*n^-1*arctanh(cosh(a+b*log(c*x^n))))
@test integrate(x^-1*csch(a+b*log(c*x^n))^2, x) == :(-1*b^-1*n^-1*coth(a+b*log(c*x^n)))
@test integrate(x^-1*csch(a+b*log(c*x^n))^3, x) == :((1/2)*b^-1*n^-1*arctanh(cosh(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*coth(a+b*log(c*x^n))*csch(a+b*log(c*x^n)))
@test integrate(x^-1*csch(a+b*log(c*x^n))^4, x) == :(b^-1*n^-1*coth(a+b*log(c*x^n))+-1//3*b^-1*n^-1*coth(a+b*log(c*x^n))^3)
@test integrate(x^-1*csch(a+b*log(c*x^n))^5, x) == :(-3//8*b^-1*n^-1*arctanh(cosh(a+b*log(c*x^n)))+-1//4*b^-1*n^-1*csch(a+b*log(c*x^n))^3*coth(a+b*log(c*x^n))+3//8*b^-1*n^-1*coth(a+b*log(c*x^n))*csch(a+b*log(c*x^n)))
@test integrate(x^-1*csch(a+b*log(c*x^n))^5//2, x) == :(-2//3*b^-1*n^-1*csch(a+b*log(c*x^n))^3//2*cosh(a+b*log(c*x^n))+2//3*im*b^-1*n^-1*(im*sinh(a+b*log(c*x^n)))^(1/2)*csch(a+b*log(c*x^n))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*a+(1/2)*im*b*log(c*x^n),2))
@test integrate(x^-1*csch(a+b*log(c*x^n))^3//2, x) == :(-2*b^-1*n^-1*csch(a+b*log(c*x^n))^(1/2)*cosh(a+b*log(c*x^n))+-2*im*b^-1*n^-1*(im*sinh(a+b*log(c*x^n)))^-1//2*csch(a+b*log(c*x^n))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*a+(1/2)*im*b*log(c*x^n),2))
@test integrate(x^-1*csch(a+b*log(c*x^n))^(1/2), x) == :(-2*im*b^-1*n^-1*(im*sinh(a+b*log(c*x^n)))^(1/2)*csch(a+b*log(c*x^n))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*a+(1/2)*im*b*log(c*x^n),2))
@test integrate(x^-1*csch(a+b*log(c*x^n))^-1//2, x) == :(-2*im*b^-1*n^-1*(im*sinh(a+b*log(c*x^n)))^-1//2*csch(a+b*log(c*x^n))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*a+(1/2)*im*b*log(c*x^n),2))
@test integrate(x^-1*csch(a+b*log(c*x^n))^-3//2, x) == :(2//3*b^-1*n^-1*csch(a+b*log(c*x^n))^-1//2*cosh(a+b*log(c*x^n))+2//3*im*b^-1*n^-1*(im*sinh(a+b*log(c*x^n)))^(1/2)*csch(a+b*log(c*x^n))^(1/2)*Elliptic.F(-1//4*pi+(1/2)*im*a+(1/2)*im*b*log(c*x^n),2))
@test integrate(x^-1*csch(a+b*log(c*x^n))^-5//2, x) == :(2//5*b^-1*n^-1*csch(a+b*log(c*x^n))^-3//2*cosh(a+b*log(c*x^n))+6//5*im*b^-1*n^-1*(im*sinh(a+b*log(c*x^n)))^-1//2*csch(a+b*log(c*x^n))^-1//2*Elliptic.E(-1//4*pi+(1/2)*im*a+(1/2)*im*b*log(c*x^n),2))

max_stars_repo_path: test/6 Hyperbolic functions/6.6 Hyperbolic cosecant/6.6.3 Hyperbolic cosecant functions.jl


text:
using BenchmarkTools, Test, Knet

m=rand(1000,1000)
c=rand(1000)
r=c'

1                               # μs
@show @benchmark hcat($m,$m)    # 1221
@show @benchmark hcat($m,$c)    # 609
@show @benchmark hcat($c,$m)    # 609
@show @benchmark hcat($c,$c)    # 8
@show @benchmark vcat($m,$m)    # 2317
@show @benchmark vcat($m,$r)    # 671
@show @benchmark vcat($r,$m)    # 672
@show @benchmark vcat($c,$c)    # 1

km = KnetArray(m)
kc = KnetArray(c)
kr = KnetArray(r)
gs() = Knet.cudaDeviceSynchronize()
macro benchmarx(x); :(@benchmark ($x;gs())); end

@show @benchmarx hcat($km,$km) # 225
@show @benchmarx hcat($km,$kc) # 122
@show @benchmarx hcat($kc,$km) # 126
@show @benchmarx hcat($kc,$kc) # 22
@show @benchmarx vcat($km,$km) # 602
@show @benchmarx vcat($km,$kr) # 316
@show @benchmarx vcat($kr,$km) # 316
@show @benchmarx vcat($kc,$kc) # 18

nothing

#=
# commit 359d3646 2018-09-22, julia 1.0.0 vs commit 4aa5f92f 2018-08-14, julia 0.6.4
# GPU:V100, CPU:Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz

@benchmark(hcat($(Expr(:$, :m)), $(Expr(:$, :m)))) = Trial(1.013 ms)	= Trial(938.682 μs)
@benchmark(hcat($(Expr(:$, :m)), $(Expr(:$, :c)))) = Trial(512.029 μs)  = Trial(516.471 μs)
@benchmark(hcat($(Expr(:$, :c)), $(Expr(:$, :m)))) = Trial(493.630 μs)  = Trial(469.592 μs)
@benchmark(hcat($(Expr(:$, :c)), $(Expr(:$, :c)))) = Trial(2.864 μs)    = Trial(7.658 μs)   ++++
@benchmark(vcat($(Expr(:$, :m)), $(Expr(:$, :m)))) = Trial(2.264 ms)    = Trial(2.259 ms)
@benchmark(vcat($(Expr(:$, :m)), $(Expr(:$, :r)))) = Trial(993.510 μs)  = Trial(761.153 μs) ----
@benchmark(vcat($(Expr(:$, :r)), $(Expr(:$, :m)))) = Trial(990.959 μs)  = Trial(760.438 μs) ----
@benchmark(vcat($(Expr(:$, :c)), $(Expr(:$, :c)))) = Trial(887.162 ns)  = Trial(976.700 ns) ++++
@benchmarx(hcat($(Expr(:$, :km)), $(Expr(:$, :km)))) = Trial(60.120 μs)   = Trial(60.645 μs)
@benchmarx(hcat($(Expr(:$, :km)), $(Expr(:$, :kc)))) = Trial(130.411 μs)  = Trial(197.104 μs) ++++
@benchmarx(hcat($(Expr(:$, :kc)), $(Expr(:$, :km)))) = Trial(135.312 μs)  = Trial(101.766 μs) ----
@benchmarx(hcat($(Expr(:$, :kc)), $(Expr(:$, :kc)))) = Trial(15.883 μs)   = Trial(17.072 μs)
@benchmarx(vcat($(Expr(:$, :km)), $(Expr(:$, :km)))) = Trial(144.315 μs)  = Trial(5.020 s)  ????
@benchmarx(vcat($(Expr(:$, :km)), $(Expr(:$, :kr)))) = Trial(170.106 μs)  = Trial(5.589 s)  ????
@benchmarx(vcat($(Expr(:$, :kr)), $(Expr(:$, :km)))) = Trial(179.613 μs)  = Trial(5.044 s)  ????
@benchmarx(vcat($(Expr(:$, :kc)), $(Expr(:$, :kc)))) = Trial(17.047 μs)   = Trial(17.001 μs)
=#


text:
function initialize_converter!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{PSY.AverageConverter, O, IC, DC, P, F},
    },
    inner_vars::AbstractVector,
) where {
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
} end

function initialize_converter!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{PSY.RenewableEnergyConverterTypeA, O, IC, DC, P, F},
    },
    inner_vars::AbstractVector,
) where {
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}
    #Get inner vars
    V_R = inner_vars[Vr_inv_var]
    V_I = inner_vars[Vi_inv_var]
    θ = atan(V_I / V_R)
    V_t = sqrt(V_R^2 + V_I^2)
    Iq_external = inner_vars[Ii_cnv_var]
    Ip_external = inner_vars[Ir_cnv_var]
    #Reference Transformation
    Ip = Ip_external * cos(-θ) - Iq_external * sin(-θ)
    Iq = Ip_external * sin(-θ) + Iq_external * cos(-θ)
    converter = PSY.get_converter(dynamic_device)
    Io_lim = PSY.get_Io_lim(converter)
    Vo_lim = PSY.get_Vo_lim(converter)

    # Lv_pnt0 is unused in the initialization
    _, Lv_pnt1 = PSY.get_Lv_pnts(converter)

    if (Iq < Io_lim) || (V_t > Vo_lim) || (V_t < Lv_pnt1)
        error("Power flow solution outside of inverter limits. Update parameters.")
    end

    #Update converter states
    converter_ix = get_local_state_ix(dynamic_device, PSY.RenewableEnergyConverterTypeA)
    converter_states = @view device_states[converter_ix]
    converter_states[1] = Ip #Ip_cnv
    converter_states[2] = -Iq #Iq_cnv
    converter_states[3] = V_t #Vmeas

    #Update entry to converter currents (inner currents)
    inner_vars[Id_ic_var] = Ip
    # Multiplied by -1 because PSS/e doesn't do a conjugate.
    inner_vars[Iq_ic_var] = -Iq
end


text:
"""
It's like `Random.shuffle`. When a `Matrix{Residue}` is used, you can indicate if the gaps
should remain their positions using the last boolean argument. The previous argument should
be the dimension to shuffle, 1 for shuffling residues in a sequence (row) or 2 for shuffling
residues in a column.

```jldoctest
julia> using MIToS.MSA

julia> using Random

julia> msa = hcat(res"RRE",res"DDK", res"G--")
3×3 Array{Residue,2}:
 R  D  G
 R  D  -
 E  K  -

julia> Random.seed!(42);

julia> shuffle(msa, 1, true)
3×3 Array{Residue,2}:
 G  D  R
 D  R  -
 E  K  -

julia> Random.seed!(42);

julia> shuffle(msa, 1, false)
3×3 Array{Residue,2}:
 G  D  R
 D  -  R
 -  E  K

```
"""
function Random.shuffle!(r::AbstractRNG, msa::Matrix{Residue},
                       dim::Int, fixedgaps::Bool=true)
    nseq, ncol = size(msa)
    @assert dim == 1 || dim == 2 "The dimension must be 1 (sequences) or 2 (columns)"
    if fixedgaps
        mask = msa .!= GAP
        if dim == 2
            @inbounds for i in 1:ncol
                shuffle!(view(msa,mask[:,i],i))
            end
        elseif dim == 1
            @inbounds for i in 1:nseq
                shuffle!(view(msa,i,mask[i,:]))
            end
        end
    else
        if dim == 2
            @inbounds for i in 1:ncol
                shuffle!(view(msa,:,i))
            end
        elseif dim == 1
            @inbounds for i in 1:nseq
                shuffle!(view(msa,i,:))
            end
        end
    end
    msa
end

function Random.shuffle!(msa::Matrix{Residue}, args...)
    shuffle!(Random.GLOBAL_RNG, msa, args...)
end

"""
It's like `shuffle` but in-place. When a `Matrix{Residue}` or a `AbstractAlignedObject`
(sequence or MSA) is used, you can indicate if the gaps should remain their positions
using the last boolean argument.
"""
function Random.shuffle(r::AbstractRNG, msa::Matrix{Residue}, args...)
    shuffle!(r, copy(msa), args...)
end

function Random.shuffle(msa::Matrix{Residue}, args...)
    shuffle!(Random.GLOBAL_RNG, copy(msa), args...)
end

function Random.shuffle(r::AbstractRNG,
                      msa::Union{AbstractAlignedObject,
                                 NamedResidueMatrix{Array{Residue,2}}}, args...)
    shuffle(r, copy(getresidues(msa)), args...)
end

function Random.shuffle(msa::Union{AbstractAlignedObject,
                                   NamedResidueMatrix{Array{Residue,2}}}, args...)
    shuffle(Random.GLOBAL_RNG, copy(getresidues(msa)), args...)
end


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ aa854f02-46bf-4a2b-b1e9-5fb52b99925c
begin
    import Pkg
    Pkg.activate(Base.current_project())
    Pkg.instantiate()
end

# ╔═╡ 8b4426b8-588b-49ad-83c5-a79ed698d704
using PlutoUI

# ╔═╡ d47c1cf4-4844-4c66-8a21-b62ad7740dcf
using JLD2

# ╔═╡ 3c07c33b-7f91-4e2d-a27a-09a8924f328c
using DataFrames

# ╔═╡ b59dd2c0-4c3e-4c3d-a4d4-f8da9dcbe286
using CSV

# ╔═╡ 76e2e659-c708-4a08-8263-def6c3dfa930
using Dates

# ╔═╡ 41b0f95e-d806-4167-8199-7392717782f9
using OffsetArrays

# ╔═╡ 342cbcff-5a05-4199-b5b5-6cc0fa70c202
using TemplateMatching

# ╔═╡ 4249c937-1b33-4aca-8351-d98a19c584fa
using Tables

# ╔═╡ 87df9037-d497-44e8-836d-3c34a58c5d07
using StatsBase

# ╔═╡ c2f22858-6d3c-4f92-a05b-659deae566f6
md"## Reading Continuous Data"

# ╔═╡ a86b6c8a-9b90-4388-8208-d257ec6951d2
md"""Path of the JLD2 file containing continuous data: $(@bind datapath TextField(default="../data/2021-01-12_20-25-30/data.jld2"))"""

# ╔═╡ f3f17c65-dab6-46da-83a9-87aad1181524
data = load(datapath, "data")

# ╔═╡ f7e0f70d-4069-4039-90b6-d3da559bdcb1
samplefreq = 2 # MHz

# ╔═╡ c7f0582a-8997-4641-b607-0b95a7f73a6c
v_p = 0.67345 # in cm/us

# ╔═╡ 57b35604-727b-4eac-bcbd-ea302e4c79a2
md"## Reading events catalogue"

# ╔═╡ b62d7318-3a6d-4442-b6cc-3166b87dff8a
md"""Path of the CSV catalogue $(@bind cataloguepath TextField(default="../data/2021-01-12_20-25-30/cat.csv"))"""

# ╔═╡ d5beeb75-7ca6-4f28-a63b-6e3070fb4c73
begin
	import Base.+, Base.show

	+(x::Dates.UTInstant{Microsecond}, d::Microsecond) = Dates.UTInstant(x.periods + d)

	function show(io::IO, ::MIME"text/plain", x::Dates.UTInstant{Microsecond})
		date = DateTime(Dates.UTM(round(Int, 1e-3 * x.periods.value)))
		show(io, MIME"text/plain"(), date)
	end
end

# ╔═╡ 1bf27aad-b754-4acc-b08d-fe1a86791de8
function DateTimeMicrosecond(y, m, d, h, mi, s, us)
	rata = us + 1_000_000 * (s + 60mi + 3600h + 86400 * Dates.totaldays(y, m, d))
	Dates.UTInstant{Microsecond}(Microsecond(rata))
end

# ╔═╡ 96cf6e3b-1c2b-49f6-979d-0fce28a55a06
starttime = DateTimeMicrosecond(2021, 01, 12, 20, 25, 30, 0)

# ╔═╡ 5d5f3ca7-260d-4dbb-bde1-59f2bd58c0c9
catalogue = let
	columns = [:Year, :Month, :Day, :Hour, :Minute, :Second, :North, :East, :Up]
	df = CSV.read(cataloguepath, DataFrame, select=columns)
	sec = floor.(Int, df.Second)
	usec = round.(Int, 1e6 * (df.Second - sec))
	catalogue = DataFrame()
	catalogue.datetime = DateTimeMicrosecond.(df.Year, df.Month, df.Day, df.Hour, df.Minute, sec, usec)
	catalogue.sample = map(x -> round(Int, x.value * samplefreq), catalogue.datetime .- starttime)
	catalogue.north = df.North
	catalogue.east = df.East
	catalogue.up = df.Up
	catalogue
end

# ╔═╡ 6ac6cce3-09ae-421a-b502-b5c8807bf555
md"## Reading Sensors Positions"

# ╔═╡ cfceae8c-0967-4a47-93f5-2b6e3073624c
md"""Path of the CSV of sensors positions $(@bind sensorspospath TextField(default="../data/2021-01-12_20-25-30/passive-xyz.csv"))"""

# ╔═╡ be090485-7917-464a-8e8c-aec4a2c009d3
sensors_positions = let
	sensors_positions = CSV.read(sensorspospath, DataFrame, header=[:north, :east, :up])
	for s = [:north, :east, :up]
		sensors_positions[!, s] .*= 100
	end
	sensors_positions.sensor = axes(sensors_positions, 1) .- 1
	filter!(r -> r.sensor in keys(data), sensors_positions)
	sensors_positions
end

# ╔═╡ 368f7fc1-0afc-439f-905a-76af04e8ca8d
md"## Processing the catalogue"

# ╔═╡ 033473d3-6958-4e7b-8c95-a9e15d882118
md"""
Template pre:$(@bind t_pre Slider(0:1000, default=100, show_value=true))
"""

# ╔═╡ 1ba147d6-2609-4d83-a734-6bf91445a23c
md"""Template post:$(@bind t_post Slider(0:1000, default=500, show_value=true))
"""

# ╔═╡ 7d1f1903-8342-4715-8e88-3e3d063ea5db
md"""
Signal height threshold $(@bind height_threshold Slider(0.0:0.05:1, default=0.3, show_value=true))
"""

# ╔═╡ 8790874d-322e-4aea-a60b-4717e653afc9
md"""Distance (in unit of template length) $(@bind reldistance Slider(0:0.5:5, default=2, show_value=true))"""

# ╔═╡ 0d122aef-1767-4389-b043-97d07ecef1fe
md"""
Tolerance $(@bind tolerance Slider(0:10, default=5, show_value=true))
"""

# ╔═╡ ef346f24-ab46-407b-95a0-2bce10b41ac8
md"""Cross-correlation threshold $(@bind cc_threshold Slider(0.0:0.05:1, default=0.5, show_value=true))"""

# ╔═╡ c4b39f91-8f0d-4f54-90b4-3c035a5b125f
md"""Minimum number of channels above threshold $(@bind nch_min Slider(0:length(data), default=4, show_value=true))"""

# ╔═╡ afe1c12f-23d3-47af-bfe8-b59ca7018593
function process_match(data_stream, template, match_offsets, guess, correlation_threshold, nch_min)
	toas = estimatetoa.(data_stream, template, match_offsets, tolerance)
	readings = copy(sensors_positions)
	readings.index = 1:nrow(readings)
	readings.toas = [(sample + t_pre) / samplefreq for (sample, _) in toas]
	readings.cc = [cc for (_, cc) in toas]
	filter!(row -> row.cc > correlation_threshold, readings)
	crosscorrelation = mean(readings.cc)
	relative_magnitude = magnitude(data_stream[readings.index], template[readings.index], match_offsets[readings.index])
	nch = nrow(readings)
	if nch >= nch_min
		readings_matrix = eachrow(Matrix(readings[!, [:north, :east, :up, :toas]]))
		candidate = locate(readings_matrix, v_p, guess)
		north, east, up, origin_time = candidate.minimizer
		multilateration_residual = candidate.minimum
	else
		north, east, up, origin_time = guess
		multilateration_residual = missing
	end
	(; north, east, up, origin_time, relative_magnitude, crosscorrelation, nch, multilateration_residual)
end

# ╔═╡ 666c8d5f-c1c5-4cb4-9364-bf997b18a683
augmented_catalogue = let
	data_stream = collect(values(data))
	matches_vec = Vector{DataFrame}(undef, nrow(catalogue))
	distance_threshold = reldistance * (t_pre + t_post - 1)
	Threads.@threads for n = eachindex(matches_vec)
		# Preparing template
		displacement = DataFrame()
		for s in [:north, :east, :up]
			displacement[!, s] = sensors_positions[!, s] .- catalogue[n, s]
		end
		distances = map(Base.splat(hypot), eachrow(displacement[!, [:north, :east, :up]]))
		offsets = @. round(Int, samplefreq * (distances / v_p))
		shifts = catalogue[n, :sample] .+ offsets
		cuts = [shift - t_pre: shift + t_post for shift in shifts]
		template_data = [view(series, cut) for (series, cut) in zip(data_stream, cuts)]
		# Computing cross-correlation and finding peaks
		signal = correlatetemplate(data_stream, template_data, offsets, tolerance)
		peaks, heights = findpeaks(signal, height_threshold, distance_threshold)
		# Computing magnitudes and relocating events
		matches = DataFrame()
		matches.peak_sample = peaks .+ t_pre
		matches.peak_height = heights
		matches.template .= n
		x0 = Vector(catalogue[n, [:north, :east, :up]])
		matches_data = [process_match(data_stream,
			                          template_data,
			                          peak .+ offsets,
			                          [x0; (peak + t_pre) / samplefreq],
			                          cc_threshold,
			                          nch_min)
			            for peak in peaks]
		matches = hcat(matches, DataFrame(matches_data))
		matches_vec[n] = matches
	end
	reduce(vcat, matches_vec)
end

# ╔═╡ b50df30e-f7c3-44f2-b759-3e4ca31003d9
md"""Path where to output the augmented catalogue $(@bind outputpath TextField(default="../data/2021-01-12_20-25-30/augmented_catalogue.csv"))"""

# ╔═╡ 3784a5b2-8621-4bac-87a1-4decc25cf4f5
CSV.write(outputpath, augmented_catalogue)

# ╔═╡ Cell order:
# ╠═aa854f02-46bf-4a2b-b1e9-5fb52b99925c
# ╠═8b4426b8-588b-49ad-83c5-a79ed698d704
# ╟─c2f22858-6d3c-4f92-a05b-659deae566f6
# ╟─a86b6c8a-9b90-4388-8208-d257ec6951d2
# ╠═d47c1cf4-4844-4c66-8a21-b62ad7740dcf
# ╠═f3f17c65-dab6-46da-83a9-87aad1181524
# ╠═96cf6e3b-1c2b-49f6-979d-0fce28a55a06
# ╠═f7e0f70d-4069-4039-90b6-d3da559bdcb1
# ╠═c7f0582a-8997-4641-b607-0b95a7f73a6c
# ╟─57b35604-727b-4eac-bcbd-ea302e4c79a2
# ╠═3c07c33b-7f91-4e2d-a27a-09a8924f328c
# ╠═b59dd2c0-4c3e-4c3d-a4d4-f8da9dcbe286
# ╠═76e2e659-c708-4a08-8263-def6c3dfa930
# ╟─b62d7318-3a6d-4442-b6cc-3166b87dff8a
# ╠═1bf27aad-b754-4acc-b08d-fe1a86791de8
# ╠═d5beeb75-7ca6-4f28-a63b-6e3070fb4c73
# ╠═5d5f3ca7-260d-4dbb-bde1-59f2bd58c0c9
# ╟─6ac6cce3-09ae-421a-b502-b5c8807bf555
# ╟─cfceae8c-0967-4a47-93f5-2b6e3073624c
# ╠═be090485-7917-464a-8e8c-aec4a2c009d3
# ╠═368f7fc1-0afc-439f-905a-76af04e8ca8d
# ╠═41b0f95e-d806-4167-8199-7392717782f9
# ╠═342cbcff-5a05-4199-b5b5-6cc0fa70c202
# ╠═4249c937-1b33-4aca-8351-d98a19c584fa
# ╠═87df9037-d497-44e8-836d-3c34a58c5d07
# ╟─033473d3-6958-4e7b-8c95-a9e15d882118
# ╟─1ba147d6-2609-4d83-a734-6bf91445a23c
# ╟─7d1f1903-8342-4715-8e88-3e3d063ea5db
# ╟─8790874d-322e-4aea-a60b-4717e653afc9
# ╟─0d122aef-1767-4389-b043-97d07ecef1fe
# ╟─ef346f24-ab46-407b-95a0-2bce10b41ac8
# ╟─c4b39f91-8f0d-4f54-90b4-3c035a5b125f
# ╠═afe1c12f-23d3-47af-bfe8-b59ca7018593
# ╠═666c8d5f-c1c5-4cb4-9364-bf997b18a683
# ╟─b50df30e-f7c3-44f2-b759-3e4ca31003d9
# ╠═3784a5b2-8621-4bac-87a1-4decc25cf4f5


text:
"""
    gavrikov(δ,θ,Tvn,T₀)

Correlation function for detonation cell width proposed by Gavrikov et al
COMBUSTION AND FLAME 120:19�33 (2000) based on using a reaction zone length based on time
to limiting reactant consumption in constant volume explosion approximation using vn
postshock velocity to convert time to distance.
Tested against a range of fuel-oxidizer diluent mixtures

INPUT:
    δ = reaction zone length based on time to 50consumption of limiting
    reactant from CV computation and delta = time * w_VN
    θ = Eₐ/RT_VN,  effective reduced activation energy based on CV
    computation
    Tvn = von Neumann (postshock temperature behind CJ shock wave)
    T₀ = initial temperature
"""
function gavrikov(δ,θ,Tvn,T₀)

    # Constants
    a = -0.007843787493
    b = 0.1777662961
    c = 0.02371845901
    d = 1.477047968
    e = 0.1545112957
    f = 0.01547021569
    g = -1.446582357
    h = 8.730494354
    i = 4.599907939
    j = 7.443410379
    k = 0.4058325462
    m = 1.453392165
    # Define nondimensional parameters
    X = θ
    Y = Tvn/T₀
    z = Y*(a*Y-b) + X*(c*X-d + (e-f*Y)*Y) + g*log(Y) + h*log(X) + Y*(i/X - k*Y/X^m) - j
    λ = δ*10^z
end

"""
     ng(Δᵢ,χ)

Correlation function for detonation cell size from
Ng, Hoi Dick, Yiguang Ju, and John H. S. Lee. 2007. Assessment of Detonation Hazards in
High-Pressure Hydrogen Storage from Chemical Sensitivity Analysis.
INTERNATIONAL JOURNAL OF HYDROGEN ENERGY 32 (1):93-99.
Tested only against low pressure H2-air data.

INPUT:
    Δᵢ = induction zone length based on peak thermicity in ZND simulation
    χ = 𝜺ᵢ*Δᵢ/Δᵣ where
          𝜺ᵢ = reduced effective activation energy from CV computation
          Δᵢ = distance to peak thermicity from ZND computation
          Δᵣ = w_vN/σ̇_max from ZND computation

See Ng et al. Combustion Theory and Modeling 2005 for a discussion of
the χ parameter.
"""
function ng(Δᵢ,χ)

    # Constants
    A₀ = 30.465860763763
    a₁ = 89.55438805808153
    a₂ = -130.792822369483
    a₃ = 42.02450507117405
    b₁ = -0.02929128383850
    b₂ = 1.0263250730647101E-5
    b₃ = -1.031921244571857E-9
    λ = Δᵢ*(A₀ + ((a₃/χ + a₂)/χ + a₁)/χ + ((b₃*χ + b₂)*χ + b₁)*χ)
end


text:
@testset ExtendedTestSet "Integrated" begin
    @testset ExtendedTestSet "number" begin
        v1 = [
            -1024 -1024 -1024 -1024
            -1024 800 800 -1024
            -1024 800 800 -1024
            -1024 -1024 -1024 -1024
        ]
        vol = cat(v1, v1; dims=3)
        alg = Integrated(vol)
        S_Bkg = -1024
        S_Obj = 800
        answer = 8
        @test score(S_Bkg, S_Obj, alg) == answer
    end

    @testset ExtendedTestSet "volume" begin
        v1 = [
            -1024 -1024 -1024 -1024
            -1024 800 800 -1024
            -1024 800 800 -1024
            -1024 -1024 -1024 -1024
        ]
        vol = cat(v1, v1; dims=3)
        alg = Integrated(vol)
        S_Bkg = -1024
        S_Obj = 800
        size = [0.5, 0.5, 0.5]
        answer = 8 * 0.5 * 0.5 * 0.5
        @test score(S_Bkg, S_Obj, size, alg) == answer
    end

    @testset ExtendedTestSet "mass" begin
        v1 = [
            -1024 -1024 -1024 -1024
            -1024 800 800 -1024
            -1024 800 800 -1024
            -1024 -1024 -1024 -1024
        ]
        vol = cat(v1, v1; dims=3)
        alg = Integrated(vol)
        S_Bkg = -1024
        S_Obj = 800
        size = [0.5, 0.5, 0.5]
        ρ = 10.0
        answer = 8 * 0.5 * 0.5 * 0.5 * 10
        @test score(S_Bkg, S_Obj, size, ρ, alg) == answer
    end
end


text:
@testset "NegBinomialLikelihood{<:LogisticLink}" begin
    test_interface(NegBinomialLikelihood(10.0), NegativeBinomial)
    test_auglik(NegBinomialLikelihood(LogisticLink(), 10))
    test_auglik(NegBinomialLikelihood(LogisticLink(), 5.5))
end


text:
using Test
using Suppressor

@testset "Images" begin
include("arrays.jl")
include("algorithms.jl")
@suppress_err include("exposure.jl") # deprecated
include("edge.jl")
include("corner.jl")
include("writemime.jl")
# include("deprecated.jl")

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -25.420040452488465 <= q <= 21.65452636072683)

                     add_NL_constraint(m, :(gelu(-0.9974376164994152 + 0.7314791518637902*gelu(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + 0.4470938857218645*gelu(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + 0.10010754470216687*gelu(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.8594790184444285*gelu(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + -0.4658863904521331*gelu(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + gelu(0.15396139316976853 + -0.4228510094486002*gelu(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.8035341866400141*gelu(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.07617073018159282*gelu(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + 0.3940887684836425*gelu(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.6667343490308526*gelu(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + gelu(0.9665409577623398 + -0.449175249339572*gelu(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.2640107596849468*gelu(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + 0.7284419447850534*gelu(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + 0.5662248826250997*gelu(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + -0.21454884110410477*gelu(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + gelu(-0.1582635851091898 + 0.5103556986884796*gelu(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.6017569427714062*gelu(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.23749573380464106*gelu(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.2344836196764124*gelu(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.07066758651932892*gelu(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + gelu(-0.47638200661148833 + 0.9126631115783219*gelu(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + 0.35414400688245484*gelu(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.36037763135774403*gelu(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.18657915153735738*gelu(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.7872605277677782*gelu(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Distributions, NetCDF

"""
    FuzzyClassification(M,Sinv,Rrs)

Apply fuzzy membership classifier (`M` + `Sinv`) to a vector of remotely sensed
reflectcances (`Rrs`) and returns a membership vector (values between 0 and 1).
"""
function FuzzyClassification(M,Sinv,Rrs)
        f=Array{Any,1}(undef,length(M))
        for ii=1:length(M)
                X=vec(Rrs)-M[ii]
                Z=transpose(X)*Sinv[ii]*X
                f[ii]=ccdf(Chisq(6),Z)
        end

        return f
end

"""
    Jackson2017()

`Fuzzy logic` classifiers defined in [Moore et al 2009](https://doi.org/10.1016/j.rse.2009.07.016) and
[Jackson et al 2017](http://dx.doi.org/10.1016/j.rse.2017.03.036) can be used to
assign optical class memberships from an `Rrs` vector. This function provides
vector `M` and inverse matrix `Sinv` for the Jackson et al 2017 classifier.

Credits: T. Jackson kindly provided the `J17.nc` classifier file
"""
function Jackson2017()
        fil=joinpath(dirname(pathof(OceanColorData)),"../examples/J17.nc")
        tmpM = ncread(fil, "cluster_means")
        tmpSinv = ncread(fil, "inverse_covariance")

        M=Array{Any,1}(undef,14)
        Sinv=Array{Any,1}(undef,14)
        for ii=1:length(M)
                M[ii]=vec(tmpM[ii,:])
                Sinv[ii]=tmpSinv[1:6,1:6,ii]
        end

        return M,Sinv
end

"""
    Moore2009()

`Fuzzy logic` classifiers defined in [Moore et al 2009](https://doi.org/10.1016/j.rse.2009.07.016) and
[Jackson et al 2017](http://dx.doi.org/10.1016/j.rse.2017.03.036) can be used to
assign optical class memberships from an `Rrs` vector. This function provides
vector `M` and inverse matrix `Sinv` for the Moore et al 2009 classifier.
"""
function Moore2009()

        #Moore et al 2009 mean vectors:
        M=Array{Any,1}(undef,8)
        M[1]=vec([0.0234 0.0192 0.0129 0.0075 0.0031 0.0002])
        M[2]=vec([0.0162 0.0141 0.0112 0.0073 0.0034 0.0002])
        M[3]=vec([0.0107 0.0098 0.0092 0.0070 0.0039 0.0003])
        M[4]=vec([0.0065 0.0064 0.0070 0.0064 0.0048 0.0006])
        M[5]=vec([0.0033 0.0034 0.0042 0.0042 0.0043 0.0009])
        M[6]=vec([0.0064 0.0074 0.0105 0.0116 0.0140 0.0041])
        M[7]=vec([0.0121 0.0140 0.0192 0.0204 0.0231 0.0084])
        M[8]=vec([0.0184 0.0230 0.0333 0.0359 0.0409 0.0137])

        #Moore et al 2009 covariance matrices:
        S=Array{Any,1}(undef,8)
        S[1]=[0.00000959 0.00000556 0.00000138 -0.00000034 -0.00000024 -0.00000003;
        0.00000556 0.00000493 0.00000193 0.00000060 0.00000023 0.00000001;
        0.00000138 0.00000193 0.00000282 0.00000223 0.00000119 0.00000007;
        -0.00000034 0.00000060 0.00000223 0.00000232 0.00000119 0.00000007;
        -0.00000024 0.00000023 0.00000119 0.00000119 0.00000071 0.00000005;
        -0.00000003 0.00000001 0.00000007 0.00000007 0.00000005 0.00000001]
        S[2]=[0.00000346 0.00000186 -0.00000011 -0.00000060 -0.00000062 -0.00000005;
        0.00000186 0.00000228 0.00000086 0.00000033 -0.00000007 0.00000003;
        -0.00000011 0.00000086 0.00000231 0.00000221 0.00000145 0.00000023;
        -0.00000060 0.00000033 0.00000221 0.00000266 0.00000191 0.00000034;
        -0.00000062 -0.00000007 0.00000145 0.00000191 0.00000175 0.00000041;
        -0.00000005 0.00000003 0.00000023 0.00000034 0.00000041 0.00000021]
        S[3]=[0.00000241 0.00000144 0.00000035 -0.00000031 -0.00000063 -0.00000006;
        0.00000144 0.00000138 0.00000076 0.00000015 -0.00000021 -0.00000001;
        0.00000035 0.00000076 0.00000161 0.00000156 0.00000121 0.00000016;
        -0.00000031 0.00000015 0.00000156 0.00000227 0.00000209 0.00000031;
        -0.00000063 -0.00000021 0.00000121 0.00000209 0.00000225 0.00000037;
        -0.00000006 -0.00000001 0.00000016 0.00000031 0.00000037 0.00000013]
        S[4]=[0.00000166 0.00000091 0.00000034 -0.00000009 -0.00000080 -0.00000025;
        0.00000091 0.00000097 0.00000071 0.00000025 -0.00000041 -0.00000015;
        0.00000034 0.00000071 0.00000118 0.00000103 0.00000072 0.00000003;
        -0.00000009 0.00000025 0.00000103 0.00000137 0.00000162 0.00000025;
        -0.00000080 -0.00000041 0.00000072 0.00000162 0.00000290 0.00000065;
        -0.00000025 -0.00000015 0.00000003 0.00000025 0.00000065 0.00000050]
        S[5]=[0.00000178 0.00000132 0.00000104 0.00000081 0.00000018 -0.00000014;
        0.00000132 0.00000127 0.00000121 0.00000099 0.00000034 -0.00000010;
        0.00000104 0.00000121 0.00000150 0.00000142 0.00000110 0.00000013;
        0.00000081 0.00000099 0.00000142 0.00000158 0.00000177 0.00000042;
        0.00000018 0.00000034 0.00000110 0.00000177 0.00000351 0.00000131;
        -0.00000014 -0.00000010 0.00000013 0.00000042 0.00000131 0.00000081]
        S[6]=[0.00000715 0.00000586 0.00000409 0.00000292 0.00000005 -0.00000075;
        0.00000586 0.00000589 0.00000520 0.00000398 0.00000027 -0.00000114;
        0.00000409 0.00000520 0.00000634 0.00000541 0.00000188 -0.00000097;
        0.00000292 0.00000398 0.00000541 0.00000528 0.00000392 0.00000070;
        0.00000005 0.00000027 0.00000188 0.00000392 0.00000995 0.00000657;
        -0.00000075 -0.00000114 -0.00000097 0.00000070 0.00000657 0.00000819]
        S[7]=[0.00002625 0.00001981 0.00001058 0.00000544 -0.00000654 -0.00001122;
        0.00001981 0.00001745 0.00001314 0.00000822 -0.00000431 -0.00001228;
        0.00001058 0.00001314 0.00001629 0.00001226 0.00000035 -0.00001311;
        0.00000544 0.00000822 0.00001226 0.00001170 0.00000742 -0.00000500;
        -0.00000654 -0.00000431 0.00000035 0.00000742 0.00002241 0.00001782;
        -0.00001122 -0.00001228 -0.00001311 -0.00000500 0.00001782 0.00003987]
        S[8]=[0.00001186 0.00001134 0.00001139 0.00000919 0.00000395 -0.00000186;
        0.00001134 0.00001484 0.00002034 0.00001907 0.00001531 0.00000087;
        0.00001139 0.00002034 0.00003467 0.00003546 0.00003555 0.00000708;
        0.00000919 0.00001907 0.00003546 0.00003907 0.00004604 0.00001733;
        0.00000395 0.00001531 0.00003555 0.00004604 0.00007306 0.00004953;
        -0.00000186 0.00000087 0.00000708 0.00001733 0.00004953 0.00006542];

        Sinv=inv.(S)

        return M,Sinv
end


text:
__precompile__()

module LightGBM

using Libdl

const LGBM_library = find_library(["lib_lightgbm.so", "lib_lightgbm.dll", "lib_lightgbm.dylib"], [ENV["LIGHTGBM_PATH"]])

function __init__()
   global LGBM_library
   if !haskey(ENV, "LIGHTGBM_PATH")
      error("Environment variable LIGHTGBM_PATH not found. ",
         "Set this variable to point to the LightGBM directory prior to loading LightGBM.jl ",
         "(e.g. `ENV[\"LIGHTGBM_PATH\"] = \"../LightGBM\"`).")
   else
      slib = find_library(["lib_lightgbm.so", "lib_lightgbm.dll", "lib_lightgbm.dylib"], [ENV["LIGHTGBM_PATH"]])
      if slib == ""
         error("Could not open the LightGBM library at $(ENV["LIGHTGBM_PATH"]). ",
            "Set this variable to point to the LightGBM directory prior to loading LightGBM.jl ",
            "(e.g. `ENV[\"LIGHTGBM_PATH\"] = \"../LightGBM\"`).")
      end
   end
end

using Compat

include("wrapper.jl")
include("estimators.jl")
include("utils.jl")
include("fit.jl")
include("predict.jl")
include("cv.jl")
include("search_cv.jl")

export fit, predict, cv, search_cv, savemodel, loadmodel
export LGBMEstimator, LGBMRegression, LGBMBinary, LGBMMulticlass

end # module LightGBM


text:
# ------------------------------------------------------------------------------
# multigrid
# ------------------------------------------------------------------------------

"""
```julia
Multigrid(fineoperator, coarseoperator, smoother, prolongation, multigridtype)
```

Multigrid preconditioner for finite element operators

# Arguments:
- `fineoperator`:       finite element operator to precondition
- `coarseoperator`:     coarse grid representation of finite element operator to
                            precondition
- `smoother`:           error relaxation operator, such as Jacobi
- `prolongationbases`:  element prolongation bases from coarse to fine grid

# Returns:
- Multigrid preconditioner object
"""
mutable struct Multigrid <: AbstractPreconditioner
    # data never changed
    fineoperator::Operator
    coarseoperator::Any
    smoother::AbstractPreconditioner
    prolongationbases::AbstractArray{AbstractBasis}
    multigridtype::MultigridType.MgridType

    # data empty until assembled
    prolongationmatrix::AbstractArray{Float64}
    nodecoordinatedifferences::AbstractArray{Float64}

    # inner constructor
    Multigrid(
        fineoperator::Operator,
        coarseoperator::Any,
        smoother::AbstractPreconditioner,
        prolongationbases::AbstractArray,
        multigridtype::MultigridType.MgridType,
    ) = (
        # check smoother for fine grid
        if fineoperator != smoother.operator
            error("smoother must be for fine grid operator") # COV_EXCL_LINE
        end;

        # check coarse operator is operator or Multigrid
        if !isa(coarseoperator, Operator) && !isa(coarseoperator, Multigrid)
            error("coarse operator must be an operator or multigrid") # COV_EXCL_LINE
        end;

        # check agreement in number of fields
        if length(prolongationbases) != length(fineoperator.outputs) ||
           length(prolongationbases) != length(coarseoperator.outputs)
            error("operators and prolongation bases must have same number of fields") # COV_EXCL_LINE
        end;

        # check dimensions
        for basis in prolongationbases
            if fineoperator.inputs[1].basis.dimension != basis.dimension
                error("fine grid and prolongation space dimensions must agree") #COV_EXCL_LINE
            end
        end;

        # constructor
        new(fineoperator, coarseoperator, smoother, prolongationbases, multigridtype)
    )
end

# printing
# COV_EXCL_START
function Base.show(io::IO, preconditioner::Multigrid)
    if preconditioner.multigridtype == MultigridType.hmultigrid
        print(io, "h")
    elseif preconditioner.multigridtype == MultigridType.pmultigrid
        print(io, "p")
    end
    print(io, "-multigrid preconditioner")
end
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# data for computing symbols
# ------------------------------------------------------------------------------

"""
```julia
getnodecoordinateddifferences(multigrid)
```

Compute or retrieve the array of differences in coordinates between nodes

# Arguments:
- `preconditioner`:  preconditioner to compute node coordinate differences

# Returns:
- Array of differences in coordinates between nodes
"""
function getnodecoordinatedifferences(multigrid::Multigrid)
    # assemble if needed
    if !isdefined(multigrid, :nodecoordinatedifferences)
        # setup for computation
        dimension = multigrid.fineoperator.dimension
        inputcoordinates = multigrid.coarseoperator.inputcoordinates
        outputcoordinates = multigrid.fineoperator.outputcoordinates
        lengths = [
            maximum(inputcoordinates[:, d]) - minimum(inputcoordinates[:, d]) for
            d = 1:dimension
        ]

        # fill matrix
        numberrows = size(outputcoordinates)[1]
        numbercolumns = size(inputcoordinates)[1]
        nodecoordinatedifferences = zeros(numberrows, numbercolumns, dimension)
        for i = 1:numberrows, j = 1:numbercolumns, k = 1:dimension
            nodecoordinatedifferences[i, j, k] =
                (inputcoordinates[j, k] - outputcoordinates[i, k]) / lengths[k]
        end

        # store
        multigrid.nodecoordinatedifferences = nodecoordinatedifferences
    end

    # return
    return getfield(multigrid, :nodecoordinatedifferences)
end

"""
```julia
getprolongationmatrix(multigrid)
```

Compute or retrieve the prolongation matrix

# Arguments:
- `preconditioner`:  preconditioner to compute prolongation matrix

# Returns:
- Matrix prolonging from coarse grid to fine grid

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
ctofbasis = TensorH1LagrangeBasis(2, 3, 1, 2; collocatedquadrature=true);

# operators
finediffusion = GalleryOperator("diffusion", 3, 3, mesh);
coarsediffusion = GalleryOperator("diffusion", 2, 3, mesh);

# smoother
jacobi = Jacobi(finediffusion);

# preconditioner
multigrid = PMultigrid(finediffusion, coarsediffusion, jacobi, [ctofbasis]);

# verify
u = ones(ctofbasis.numbernodes);
v = multigrid.prolongationmatrix * u;
@assert v' ≈ [4. 2. 4. 2. 1. 2. 4. 2. 4.].^-1

# output

```
"""
function getprolongationmatrix(multigrid::Multigrid)
    # assemble if needed
    if !isdefined(multigrid, :prolongationmatrix)
        # setup
        numberfinenodes = 0
        numbercoarsenodes = 0
        Pblocks = []

        # field prolongation matrices
        for basis in multigrid.prolongationbases
            # number of nodes
            numberfinenodes += basis.numberquadraturepoints * basis.numbercomponents
            numbercoarsenodes += basis.numbernodes * basis.numbercomponents

            # prolongation matrices
            push!(Pblocks, basis.interpolation)
        end

        prolongationmatrix = spzeros(numberfinenodes, numbercoarsenodes)
        currentrow = 1
        currentcolumn = 1
        for Pblock in Pblocks
            prolongationmatrix[
                currentrow:currentrow+size(Pblock)[1]-1,
                currentcolumn:currentcolumn+size(Pblock)[2]-1,
            ] = Pblock
            currentrow += size(Pblock)[1]
            currentcolumn += size(Pblock)[2]
        end

        # store
        prolongationmatrix[abs.(prolongationmatrix).<10*eps()] .= 0
        dropzeros!(prolongationmatrix)
        multigrid.prolongationmatrix =
            Diagonal(multigrid.fineoperator.multiplicity)^-1 * prolongationmatrix
    end

    # return
    return getfield(multigrid, :prolongationmatrix)
end

"""
```julia
computesymbolsprolongation(multigrid, θ)
```

Compute the symbol matrix for a multigrid prolongation operator

# Arguments:
- `preconditioner`:  preconditioner to compute prolongation symbol

# Arguments:
- `multigrid`: Multigrid operator to compute prolongation symbol matrix for
- `θ`:         Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the multigrid prolongation operator
"""
function computesymbolsprolongation(multigrid::Multigrid, θ::Array)
    # setup
    dimension = multigrid.fineoperator.dimension
    rowmodemap = multigrid.fineoperator.rowmodemap
    columnmodemap = multigrid.coarseoperator.columnmodemap
    prolongationmatrix = multigrid.prolongationmatrix
    numberrows, numbercolumns = size(prolongationmatrix)
    nodecoordinatedifferences = multigrid.nodecoordinatedifferences
    symbolmatrixnodes = zeros(ComplexF64, numberrows, numbercolumns)

    # compute
    for i = 1:numberrows, j = 1:numbercolumns
        symbolmatrixnodes[i, j] =
            prolongationmatrix[i, j] *
            ℯ^(im * sum([θ[k] * nodecoordinatedifferences[i, j, k] for k = 1:dimension]))
    end
    symbolmatrixmodes = rowmodemap * symbolmatrixnodes * columnmodemap

    # return
    return symbolmatrixmodes
end

# ------------------------------------------------------------------------------
# get/set property
# ------------------------------------------------------------------------------

function Base.getproperty(multigrid::Multigrid, f::Symbol)
    if f == :prolongationmatrix
        return getprolongationmatrix(multigrid)
    elseif f == :nodecoordinatedifferences
        return getnodecoordinatedifferences(multigrid)
    elseif f == :outputs # Used if nesting multigrid levels
        return multigrid.fineoperator.outputs
    elseif f == :rowmodemap # Used if nesting multigrid levels
        return multigrid.fineoperator.rowmodemap
    elseif f == :columnmodemap # Used if nesting multigrid levels
        return multigrid.fineoperator.columnmodemap
    elseif f == :inputcoordinates # Used if nesting multigrid levels
        return multigrid.fineoperator.inputcoordinates
    else
        return getfield(multigrid, f)
    end
end

function Base.setproperty!(multigrid::Multigrid, f::Symbol, value)
    if f == :fineoperator ||
       f == :coarseoperator ||
       f == :smoother ||
       f == :prolongationbases ||
       f == :multigridtype
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(multigrid, f, value)
    end
end

# ------------------------------------------------------------------------------
# compute symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbols(multigrid, p, v, θ)
```

Compute or retrieve the symbol matrix for a Jacobi preconditioned operator

# Arguments:
- `multigrid`:  multigrid preconditioner to compute symbol matrix for
- `p`:          smoothing paramater array
- `v`:          pre and post smooths iteration count array, 0 indicates no pre or post smoothing
- `θ`:          Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the multigrid preconditioned operator

# Example:
```jldoctest
using LinearAlgebra

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    ctofbasis = TensorH1LagrangeBasis(3, 5, 1, dimension; collocatedquadrature=true);

    # operators
    finediffusion = GalleryOperator("diffusion", 5, 5, mesh);
    coarsediffusion = GalleryOperator("diffusion", 3, 5, mesh);

    # smoother
    jacobi = Jacobi(finediffusion);

    # preconditioner
    multigrid = PMultigrid(finediffusion, coarsediffusion, jacobi, [ctofbasis]);

    # compute symbols
    A = computesymbols(multigrid, [1.0], [1, 1], π*ones(dimension));

    # verify
    eigenvalues = real(eigvals(A));
    if dimension == 1
       @assert maximum(eigenvalues) ≈ 0.64
    elseif dimension == 2
       @assert maximum(eigenvalues) ≈ 0.9082562365654528
    elseif dimension == 3
       @assert maximum(eigenvalues) ≈ 1.4359882222222669
    end
end

# output

```
"""
function computesymbols(multigrid::Multigrid, p::Array, v::Array{Int}, θ::Array)
    # validate number of parameters
    if length(v) != 2
        Throw(error("must specify number of pre and post smooths")) # COV_EXCL_LINE
    end

    # compute component symbols
    S_f = computesymbols(multigrid.smoother, p, θ)

    P_ctof = computesymbolsprolongation(multigrid, θ)
    R_ftoc = P_ctof'

    A_f = computesymbols(multigrid.fineoperator, θ)
    A_c_inv = []
    if isa(multigrid.coarseoperator, Operator)
        A_c_inv = computesymbols(multigrid.coarseoperator, θ)^-1
    elseif isa(multigrid.coarseoperator, Multigrid)
        A_c = computesymbols(multigrid.coarseoperator.fineoperator, θ)
        A_c_inv = (I - computesymbols(multigrid.coarseoperator, p, v, θ)) * A_c^-1
    else
        Throw(error("coarse operator not supported")) # COV_EXCL_LINE
    end

    # return
    return S_f^v[2] * (I - P_ctof * A_c_inv * R_ftoc * A_f) * S_f^v[1]
end

# ------------------------------------------------------------------------------


text:
# differentiated version of IR_stab.jl

function getIRA0_diff(params::ParamType{2},
                      q::AbstractArray{Tsol,1},
                      q_dot::AbstractMatrix,
                      A0::AbstractMatrix,
                      A0_dot::Abstract3DArray) where Tsol

  nd = size(q_dot, 2)  # number of dot vectors
  numDofPerNode = size(q, 1)

  @debug1 begin
    @assert size(q, 1) == size(q_dot, 1)
    @assert size(A0_dot, 3) >= nd
    @assert size(A0, 1) == size(A0_dot, 1)
    @assert size(A0, 2) == size(A0_dot, 2)
  end

  p_dot = params.get_ira0data.p_dot
  gamma = params.gamma
  gamma_1 = params.gamma_1

  # p_dot is needed for i=1:nd, so compute the derivative with respect to
  # q[1:4] and then multiply with q_dot[:, i] later
  p = calcPressure_diff(params, q, p_dot)

  rho = q[1]
  rhou = q[2]
  rhov = q[3]
  rhoe = q[4]

  rhoinv = 1/rho

  h = (rhoe + p)*rhoinv  # this isn't really h, but including the factor of
                         # 1/rho is convenient
  a2 = gamma*p*rhoinv  # speed of sound

  A0[1,1] = rho
  A0[2,1] = rhou
  A0[3,1] = rhov
  A0[4,1] = rhoe

  A0[1,2] = rhou
  A0[2,2] = rhou*rhou*rhoinv + p
  A0[3,2] = rhou*rhov*rhoinv
  A0[4,2] = rhou*h

  A0[1,3] = rhov
  A0[2,3] = rhou*rhov*rhoinv
  A0[3,3] = rhov*rhov*rhoinv + p
  A0[4,3] = rhov*h

  A0[1,4] = rhoe
  A0[2,4] = h*rhou
  A0[3,4] = h*rhov
  A0[4,4] = rho*h*h - a2*p/gamma_1

  for i=1:nd

    # compute p_dot
    p_dot_i = zero(Tsol)
    for j=1:numDofPerNode
      p_dot_i += p_dot[j]*q_dot[j, i]
    end

    rho_dot  = q_dot[1, i]
    rhou_dot = q_dot[2, i]
    rhov_dot = q_dot[3, i]
    rhoe_dot = q_dot[4, i]

    rhoinv = 1/rho
    rhoinv_dot = -rhoinv*rhoinv*rho_dot

    h = (rhoe + p)*rhoinv
    h_dot = (rhoe + p)*rhoinv_dot + (rhoe_dot + p_dot_i)*rhoinv

    a2 = gamma*p*rhoinv
    a2_dot = gamma*(p_dot_i*rhoinv + p*rhoinv_dot)

    A0_dot[1,1,i] = rho_dot
    A0_dot[2,1,i] = rhou_dot
    A0_dot[3,1,i] = rhov_dot
    A0_dot[4,1,i] = rhoe_dot

    A0_dot[1,2,i] = rhou_dot
    A0_dot[2,2,i] = 2*rhou*rhoinv*rhou_dot + rhou*rhou*rhoinv_dot + p_dot_i
    A0_dot[3,2,i] = rhov*rhoinv*rhou_dot   + rhou*rhoinv*rhov_dot + rhou*rhov*rhoinv_dot
    A0_dot[4,2,i] = h*rhou_dot + rhou*h_dot

    A0_dot[1,3,i] = rhov_dot
    A0_dot[2,3,i] = rhov*rhoinv*rhou_dot   + rhou*rhoinv*rhov_dot + rhou*rhov*rhoinv_dot
    A0_dot[3,3,i] = 2*rhov*rhoinv*rhov_dot + rhov*rhov*rhoinv_dot + p_dot_i
    A0_dot[4,3,i] = h*rhov_dot + rhov*h_dot

    A0_dot[1,4,i] = rhoe_dot
    A0_dot[2,4,i] = rhou*h_dot  + h*rhou_dot
    A0_dot[3,4,i] = rhov*h_dot  + h*rhov_dot
    A0_dot[4,4,i] = h*h*rho_dot + rho*2*h*h_dot - (p*a2_dot + p_dot_i*a2)/gamma_1
  end

  return nothing
end

function getIRA0_diff(params::ParamType{3},
                      q::AbstractArray{Tsol,1},
                      q_dot::AbstractMatrix,
                      A0::AbstractMatrix,
                      A0_dot::Abstract3DArray) where Tsol

  nd = size(q_dot, 2)  # number of dot vectors
  numDofPerNode = size(q, 1)

  @debug1 begin
    @assert size(q, 1) == size(q_dot, 1)
    @assert size(A0_dot, 3) == nd
    @assert size(A0, 1) == size(A0_dot, 1)
    @assert size(A0, 2) == size(A0_dot, 2)
  end

  p_dot = zeros(Tsol, numDofPerNode)
  gamma = params.gamma
  gamma_1 = params.gamma_1

  # p_dot is needed for i=1:nd, so compute the derivative with respect to
  # q[1:4] and then multiply with q_dot[:, i] later
  p = calcPressure_diff(params, q, p_dot)

  rho = q[1]
  rhou = q[2]
  rhov = q[3]
  rhow = q[4]
  rhoe = q[5]

  rhoinv = 1/rho

  h = (rhoe + p)*rhoinv
  a2 = gamma*p*rhoinv  # speed of sound

  A0[1,1] = rho
  A0[2,1] = rhou
  A0[3,1] = rhov
  A0[4,1] = rhow
  A0[5,1] = rhoe

  A0[1,2] = rhou
  A0[2,2] = rhou*rhou*rhoinv + p
  A0[3,2] = rhou*rhov*rhoinv
  A0[4,2] = rhou*rhow*rhoinv
  A0[5,2] = rhou*h

  A0[1,3] = rhov
  A0[2,3] = rhou*rhov/rho
  A0[3,3] = rhov*rhov*rhoinv + p
  A0[4,3] = rhov*rhow*rhoinv
  A0[5,3] = rhov*h

  A0[1,4] = rhow
  A0[2,4] = rhow*rhou*rhoinv
  A0[3,4] = rhow*rhov*rhoinv
  A0[4,4] = rhow*rhow*rhoinv + p
  A0[5,4] = rhow*h

  A0[1,5] = rhoe
  A0[2,5] = h*rhou
  A0[3,5] = h*rhov
  A0[4,5] = h*rhow
  A0[5,5] = rho*h*h - a2*p/gamma_1

  for i=1:nd

    # compute p_dot
    p_dot_i = zero(Tsol)
    for j=1:numDofPerNode
      p_dot_i += p_dot[j]*q_dot[j, i]
    end

    rho_dot  = q_dot[1, i]
    rhou_dot = q_dot[2, i]
    rhov_dot = q_dot[3, i]
    rhow_dot = q_dot[4, i]
    rhoe_dot = q_dot[5, i]

    rhoinv = 1/rho
    rhoinv_dot = -1*rhoinv*rhoinv*rho_dot

    h = (rhoe + p)*rhoinv
    h_dot = (rhoe + p)*rhoinv_dot + (rhoe_dot + p_dot_i)*rhoinv

    a2 = gamma*p*rhoinv
    a2_dot = gamma*(p_dot_i*rhoinv + p*rhoinv_dot)

    A0_dot[1,1,i] = rho_dot
    A0_dot[2,1,i] = rhou_dot
    A0_dot[3,1,i] = rhov_dot
    A0_dot[4,1,i] = rhow_dot
    A0_dot[5,1,i] = rhoe_dot

    A0_dot[1,2,i] = rhou_dot
    A0_dot[2,2,i] = 2*rhou*rhoinv*rhou_dot + rhou*rhou*rhoinv_dot + p_dot_i
    A0_dot[3,2,i] = rhov*rhoinv*rhou_dot   + rhou*rhoinv*rhov_dot + rhou*rhov*rhoinv_dot
    A0_dot[4,2,i] = rhow*rhoinv*rhou_dot   + rhou*rhoinv*rhow_dot + rhou*rhow*rhoinv_dot
    A0_dot[5,2,i] = h*rhou_dot + rhou*h_dot

    A0_dot[1,3,i] = rhov_dot
    A0_dot[2,3,i] = rhov*rhoinv*rhou_dot   + rhou*rhoinv*rhov_dot + rhou*rhov*rhoinv_dot
    A0_dot[3,3,i] = 2*rhov*rhoinv*rhov_dot + rhov*rhov*rhoinv_dot + p_dot_i
    A0_dot[4,3,i] = rhov*rhoinv*rhow_dot   + rhow*rhoinv*rhov_dot + rhov*rhow*rhoinv_dot
    A0_dot[5,3,i] = h*rhov_dot + rhov*h_dot

    A0_dot[1,4,i] = rhow_dot
    A0_dot[2,4,i] = rhow*rhoinv*rhou_dot   + rhou*rhoinv*rhow_dot + rhou*rhow*rhoinv_dot
    A0_dot[3,4,i] = rhow*rhoinv*rhov_dot   + rhov*rhoinv*rhow_dot + rhow*rhov*rhoinv_dot
    A0_dot[4,4,i] = 2*rhow*rhoinv*rhow_dot + rhow*rhow*rhoinv_dot + p_dot_i
    A0_dot[5,4,i] = h*rhow_dot + rhow*h_dot

    A0_dot[1,5,i] = rhoe_dot
    A0_dot[2,5,i] = rhou*h_dot  + h*rhou_dot
    A0_dot[3,5,i] = rhov*h_dot  + h*rhov_dot
    A0_dot[4,5,i] = rhow*h_dot  + h*rhow_dot
    A0_dot[5,5,i] = h*h*rho_dot + rho*2*h*h_dot - (p*a2_dot + p_dot_i*a2)/gamma_1
  end

  return nothing
end

"""
  Reverse mode wrt q of [`getIRA0`](@ref).

  **Inputs**

   * params
   * q: solution vector
   * A0_bar: seed matrix for reverse mode

  **Inputs/Outputs**

   * q_bar: dual part of `q`, to be updated (not overwritten), length `numDofPerNode`
   * A0: matrix, `numDofPerNode` to be overwritten
"""
function getIRA0_revq(params::ParamType{2},
                     q::AbstractArray{Tsol,1},
                     q_bar::AbstractArray{Tsol, 1},
                     A0::AbstractMatrix{Tsol},
                     A0_bar::AbstractMatrix{Tsol}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1
  p = calcPressure(params, q)

  rho = q[1]
  rhou = q[2]
  rhov = q[3]
  rhoe = q[4]

  rhoinv = 1/rho

  h = (rhoe + p)*rhoinv  # this isn't really h, but including the factor of
                         # 1/rho is convenient
  a2 = gamma*p*rhoinv  # speed of sound

  A0[1,1] = rho
  A0[2,1] = rhou
  A0[3,1] = rhov
  A0[4,1] = rhoe

  A0[1,2] = rhou
  A0[2,2] = rhou*rhou*rhoinv + p
  A0[3,2] = rhou*rhov*rhoinv
  A0[4,2] = rhou*h

  A0[1,3] = rhov
  A0[2,3] = rhou*rhov*rhoinv
  A0[3,3] = rhov*rhov*rhoinv + p
  A0[4,3] = rhov*h

  A0[1,4] = rhoe
  A0[2,4] = h*rhou
  A0[3,4] = h*rhov
  A0[4,4] = rho*h*h - a2*p/gamma_1

  # reverse sweep
  rho_bar = zero(Tsol)
  rhou_bar = zero(Tsol)
  rhov_bar = zero(Tsol)
  rhoe_bar = zero(Tsol)
  h_bar = zero(Tsol)
  rhoinv_bar = zero(Tsol)
  p_bar = zero(Tsol)
  a2_bar = zero(Tsol)

  # first column
  rho_bar  += A0_bar[1,1]
  rhou_bar += A0_bar[2,1]
  rhov_bar += A0_bar[3,1]
  rhoe_bar += A0_bar[4,1]

  # second column
  rhou_bar += A0_bar[1,2]

  rhou_bar   += 2*rhou*rhoinv*A0_bar[2,2]
  rhoinv_bar += rhou*rhou*A0_bar[2,2]
  p_bar      += A0_bar[2,2]

  rhou_bar   += rhov*rhoinv*A0_bar[3,2]
  rhov_bar   += rhou*rhoinv*A0_bar[3,2]
  rhoinv_bar += rhou*rhov*A0_bar[3,2]

  rhou_bar += h*A0_bar[4,2]
  h_bar    += rhou*A0_bar[4,2]

  # 3rd column
  rhov_bar += A0_bar[1,3]

  rhou_bar   += rhov*rhoinv*A0_bar[2,3]
  rhov_bar   += rhou*rhoinv*A0_bar[2,3]
  rhoinv_bar += rhou*rhov*A0_bar[2,3]

  rhov_bar   += 2*rhov*rhoinv*A0_bar[3,3]
  rhoinv_bar += rhov*rhov*A0_bar[3,3]
  p_bar      += A0_bar[3,3]

  rhov_bar += h*A0_bar[4,3]
  h_bar    += rhov*A0_bar[4,3]

  # 4th column
  rhoe_bar += A0_bar[1,4]

  h_bar     += rhou*A0_bar[2,4]
  rhou_bar  +=  h*A0_bar[2,4]

  h_bar    += rhov*A0_bar[3,4]
  rhov_bar += h*A0_bar[3,4]

  rho_bar += h*h*A0_bar[4,4]
  h_bar   += 2*rho*h*A0_bar[4,4]
  a2_bar  += -p*A0_bar[4,4]/gamma_1
  p_bar   += -a2*A0_bar[4,4]/gamma_1

  # a2
  p_bar      += gamma*rhoinv*a2_bar
  rhoinv_bar += gamma*p*a2_bar

  # h
  rhoe_bar   += rhoinv*h_bar
  p_bar      += rhoinv*h_bar
  rhoinv_bar += (rhoe + p)*h_bar

  # rhoinv
  rho_bar += -rhoinv*rhoinv*rhoinv_bar

  calcPressure_revq(params, q, q_bar, p_bar)

  q_bar[1] += rho_bar
  q_bar[2] += rhou_bar
  q_bar[3] += rhov_bar
  q_bar[4] += rhoe_bar

  return nothing
end

function getIRA0_revq(params::ParamType{3},
                     q::AbstractArray{Tsol,1},
                     q_bar::AbstractArray{Tsol, 1},
                     A0::AbstractMatrix{Tsol},
                     A0_bar::AbstractMatrix{Tsol}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1
  p = calcPressure(params, q)

  rho = q[1]
  rhou = q[2]
  rhov = q[3]
  rhow = q[4]
  rhoe = q[5]

  rhoinv = 1/rho

  h = (rhoe + p)*rhoinv  # this isn't really h, but including the factor of
                         # 1/rho is convenient
  a2 = gamma*p*rhoinv  # speed of sound

  A0[1,1] = rho
  A0[2,1] = rhou
  A0[3,1] = rhov
  A0[4,1] = rhow
  A0[5,1] = rhoe

  A0[1,2] = rhou
  A0[2,2] = rhou*rhou*rhoinv + p
  A0[3,2] = rhou*rhov*rhoinv
  A0[4,2] = rhou*rhow*rhoinv
  A0[5,2] = rhou*h

  A0[1,3] = rhov
  A0[2,3] = rhou*rhov*rhoinv
  A0[3,3] = rhov*rhov*rhoinv + p
  A0[4,3] = rhov*rhow*rhoinv
  A0[5,3] = rhov*h

  A0[1,4] = rhow
  A0[2,4] = rhow*rhou*rhoinv
  A0[3,4] = rhow*rhov*rhoinv
  A0[4,4] = rhow*rhow*rhoinv + p
  A0[5,4] = rhow*h

  A0[1,5] = rhoe
  A0[2,5] = h*rhou
  A0[3,5] = h*rhov
  A0[4,5] = h*rhow
  A0[5,5] = rho*h*h - a2*p/gamma_1

  # reverse sweep
  rho_bar = zero(Tsol)
  rhou_bar = zero(Tsol)
  rhov_bar = zero(Tsol)
  rhow_bar = zero(Tsol)
  rhoe_bar = zero(Tsol)
  h_bar = zero(Tsol)
  rhoinv_bar = zero(Tsol)
  p_bar = zero(Tsol)
  a2_bar = zero(Tsol)

  # first column
  rho_bar  += A0_bar[1,1]
  rhou_bar += A0_bar[2,1]
  rhov_bar += A0_bar[3,1]
  rhow_bar += A0_bar[4,1]
  rhoe_bar += A0_bar[5,1]

  # second column
  rhou_bar += A0_bar[1,2]

  rhou_bar   += 2*rhou*rhoinv*A0_bar[2,2]
  rhoinv_bar += rhou*rhou*A0_bar[2,2]
  p_bar      += A0_bar[2,2]

  rhou_bar   += rhov*rhoinv*A0_bar[3,2]
  rhov_bar   += rhou*rhoinv*A0_bar[3,2]
  rhoinv_bar += rhou*rhov*A0_bar[3,2]

  rhou_bar   += rhow*rhoinv*A0_bar[4,2]
  rhow_bar   += rhou*rhoinv*A0_bar[4,2]
  rhoinv_bar += rhou*rhow*A0_bar[4,2]

  rhou_bar += h*A0_bar[5,2]
  h_bar    += rhou*A0_bar[5,2]

  # 3rd column
  rhov_bar += A0_bar[1,3]

  rhou_bar   += rhov*rhoinv*A0_bar[2,3]
  rhov_bar   += rhou*rhoinv*A0_bar[2,3]
  rhoinv_bar += rhou*rhov*A0_bar[2,3]

  rhov_bar   += 2*rhov*rhoinv*A0_bar[3,3]
  rhoinv_bar += rhov*rhov*A0_bar[3,3]
  p_bar      += A0_bar[3,3]

  rhov_bar   += rhow*rhoinv*A0_bar[4,3]
  rhow_bar   += rhov*rhoinv*A0_bar[4,3]
  rhoinv_bar += rhov*rhow*A0_bar[4,3]

  rhov_bar += h*A0_bar[5,3]
  h_bar    += rhov*A0_bar[5,3]

  # 4th column
  rhow_bar += A0_bar[1,4]

  rhow_bar   += rhou*rhoinv*A0_bar[2,4]
  rhou_bar   += rhow*rhoinv*A0_bar[2,4]
  rhoinv_bar += rhow*rhou*A0_bar[2,4]

  rhow_bar   += rhov*rhoinv*A0_bar[3,4]
  rhov_bar   += rhow*rhoinv*A0_bar[3,4]
  rhoinv_bar += rhow*rhov*A0_bar[3,4]

  rhow_bar   += 2*rhow*rhoinv*A0_bar[4,4]
  rhoinv_bar += rhow*rhow*A0_bar[4,4]
  p_bar      += A0_bar[4,4]

  rhow_bar += h*A0_bar[5,4]
  h_bar    += rhow*A0_bar[5,4]

  # 5th column
  rhoe_bar += A0_bar[1,5]

  h_bar     += rhou*A0_bar[2,5]
  rhou_bar  +=  h*A0_bar[2,5]

  h_bar    += rhov*A0_bar[3,5]
  rhov_bar += h*A0_bar[3,5]

  h_bar    += rhow*A0_bar[4,5]
  rhow_bar += h*A0_bar[4,5]

  rho_bar += h*h*A0_bar[5,5]
  h_bar   += 2*rho*h*A0_bar[5,5]
  a2_bar  += -p*A0_bar[5,5]/gamma_1
  p_bar   += -a2*A0_bar[5,5]/gamma_1

  # a2
  p_bar      += gamma*rhoinv*a2_bar
  rhoinv_bar += gamma*p*a2_bar

  # h
  rhoe_bar   += rhoinv*h_bar
  p_bar      += rhoinv*h_bar
  rhoinv_bar += (rhoe + p)*h_bar

  # rhoinv
  rho_bar += -rhoinv*rhoinv*rhoinv_bar

  calcPressure_revq(params, q, q_bar, p_bar)

  q_bar[1] += rho_bar
  q_bar[2] += rhou_bar
  q_bar[3] += rhov_bar
  q_bar[4] += rhow_bar
  q_bar[5] += rhoe_bar

  return nothing
end