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

include("../../fancy_toys.jl") # for get_addable_spec

# return the platforms supported by libjulia
function libjulia_platforms(julia_version)
    platforms = supported_platforms(; experimental=julia_version ≥ v"1.7")

    # skip 32bit musl builds; they fail with this error:
    #    libunwind.so.8: undefined reference to `setcontext'
    filter!(p -> !(Sys.islinux(p) && libc(p) == "musl" && arch(p) == "i686"), platforms)

    # in Julia <= 1.3 skip PowerPC builds (see https://github.com/JuliaPackaging/Yggdrasil/pull/1795)
    if julia_version < v"1.4"
        filter!(p -> !(Sys.islinux(p) && arch(p) == "powerpc64le"), platforms)
    end

    return platforms
end

# Collection of sources required to build Julia
function build_julia(ARGS, version)
    name = "libjulia"

    checksums = Dict(
        v"1.3.1" => "3d9037d281fb41ad67b443f42d8a8e400b016068d142d6fafce1952253ae93db",
        v"1.4.2" => "76a94e06e68fb99822e0876a37c2ed3873e9061e895ab826fd8c9fc7e2f52795",
        v"1.5.3" => "be19630383047783d6f314ebe0bf5e3f95f82b0c203606ec636dced405aab1fe",
        v"1.6.0" => "c253360d29abb9a3a9e6e01493e0f7bb537d88014fd58ac561b5ba30fcb44cad",
        v"1.7.0-beta2" => "7bd9ae5b20f9a03ef97f830285007eaa318a780675ab59b6528cd8ef242aa8cb",
    )

    sources = [
        ArchiveSource("https://github.com/JuliaLang/julia/releases/download/v$(version)/julia-$(version).tar.gz", checksums[version]),
        DirectorySource("./bundled"),
    ]

    # Bash recipe for building across all platforms
    script = raw"""
    apk add coreutils libuv-dev utf8proc

    cd $WORKSPACE/srcdir/julia*
    version=$(cat VERSION)

    # Apply patches
    if [ -d $WORKSPACE/srcdir/patches ]; then
    for f in $WORKSPACE/srcdir/patches/*.patch; do
        echo "Applying path ${f}"
        atomic_patch -p1 ${f}
    done
    fi

    case ${target} in
        *linux*)
            OS=Linux
        ;;
        *mingw*)
            OS=WINNT
        ;;
        *darwin*)
            OS=Darwin
        ;;
        *freebsd*)
            OS=FreeBSD
        ;;
    esac

    cat << EOM >Make.host.user
    override CC=${CC_BUILD}
    override CXX=${CXX_BUILD}
    override AR=${AR_BUILD}
    USE_SYSTEM_LIBUV=1
    USE_SYSTEM_UTF8PROC=1
    # julia want's libuv.a
    override LIBUV=/usr/lib/libuv.so
    override LIBUTF8PROC=/usr/lib/libutf8proc.so.2

    override OS=Linux
    EOM

    # TODO: eventually we should get LLVM_CXXFLAGS from llvm-config from
    # a HostDependency, now that we have those
    LLVM_CXXFLAGS="-I${prefix}/include -fno-exceptions -fno-rtti -D_GNU_SOURCE -D__STDC_CONSTANT_MACROS -D__STDC_FORMAT_MACROS -D__STDC_LIMIT_MACROS"
    if [[ "${version}" == 1.[0-5].* ]]; then
        LLVM_CXXFLAGS="${LLVM_CXXFLAGS} -std=c++11""
    else
        LLVM_CXXFLAGS="${LLVM_CXXFLAGS} -std=c++14""
    fi
    LLVM_LDFLAGS="-L${prefix}/lib"
    LDFLAGS="-L${prefix}/lib"
    CFLAGS="-I${prefix}/include"
    # -NDEBUG below fixes the FreeBSD build of Julia 1.4 and 1.5
    CXXFLAGS="-I${prefix}/include -DNDEBUG"
    if [[ "${target}" == *mingw* ]]; then
        LLVMLINK="-L${prefix}/bin -lLLVM"
        LLVM_LDFLAGS="-L${prefix}/bin"
        LDFLAGS="-L${prefix}/bin"
    elif [[ "${target}" == *apple* ]]; then
        LLVMLINK="-L${prefix}/lib -lLLVM"
    else
        if [[ "${version}" == 1.3.* ]]; then
            LLVMLINK="-L${prefix}/lib -lLLVM-6.0"
        elif [[ "${version}" == 1.4.* ]]; then
            LLVMLINK="-L${prefix}/lib -lLLVM-8jl"
        elif [[ "${version}" == 1.5.* ]]; then
            LLVMLINK="-L${prefix}/lib -lLLVM-9jl"
        elif [[ "${version}" == 1.6.* ]]; then
            LLVMLINK="-L${prefix}/lib -lLLVM-11jl"
        elif [[ "${version}" == 1.7.* ]]; then
            LLVMLINK="-L${prefix}/lib -lLLVM-12jl"
        else
            echo "Error, LLVM version not specified"
            exit 1
        fi
    fi

    # enable extglob for BB_TRIPLET_LIBGFORTRAN_CXXABI
    shopt -s extglob
    if [[ "${version}" == 1.[0-5].* ]]; then
        BB_TRIPLET_LIBGFORTRAN_CXXABI=${bb_full_target/-julia_version+([^-])}
    else
        # Strip the OS version from Darwin and FreeBSD
        BB_TRIPLET_LIBGFORTRAN_CXXABI=$(echo ${bb_full_target/-julia_version+([^-])} | sed 's/\(darwin\|freebsd\)[0-9.]*/\1/')
    fi

    cat << EOM >Make.user
    USE_SYSTEM_LLVM=1
    USE_SYSTEM_LIBUNWIND=1

    USE_SYSTEM_PCRE=1
    USE_SYSTEM_OPENLIBM=1
    USE_SYSTEM_DSFMT=1
    USE_SYSTEM_GMP=1
    USE_SYSTEM_MPFR=1
    USE_SYSTEM_SUITESPARSE=1
    USE_SYSTEM_LIBUV=1
    USE_SYSTEM_UTF8PROC=1
    USE_SYSTEM_MBEDTLS=1
    USE_SYSTEM_LIBSSH2=1
    USE_SYSTEM_CURL=1
    USE_SYSTEM_LIBGIT2=1
    USE_SYSTEM_PATCHELF=1
    USE_SYSTEM_ZLIB=1
    USE_SYSTEM_P7ZIP=1

    override XC_HOST=${target}
    override OS=${OS}
    override BUILD_OS=Linux

    #llvm-config-host is not available
    override LLVMLINK=${LLVMLINK}
    override LLVM_CXXFLAGS=${LLVM_CXXFLAGS}
    override LLVM_LDFLAGS=${LLVM_LDFLAGS}

    # just nop this
    override LLVM_CONFIG_HOST=true

    # we only run flisp and we built that for Linux
    override spawn = \$(1)
    override cygpath_w = \$(1)

    # julia expects libuv-julia.a
    override LIBUV=${prefix}/lib/libuv.a

    override BB_TRIPLET_LIBGFORTRAN_CXXABI=${BB_TRIPLET_LIBGFORTRAN_CXXABI}
    override USE_BINARYBUILDER=1

    prefix=${prefix}
    LOCALBASE=${prefix}
    EOM

    # setting USE_SYSTEM_BLAS to get Julia to use OpenBLAS_jll doesn't work on macOS,
    # where the Julia build system instead tries to use a native BLAS, which however
    # then requires compiling a custom lapack, which fails for Julia 1.4
    if [[ "${target}" != *apple* ]]; then
        cat << EOM >>Make.user
        USE_SYSTEM_BLAS=1
        LIBBLASNAME=libopenblas
        USE_SYSTEM_LAPACK=1
        LIBLAPACKNAME=libopenblas
    EOM
    elif [[ "${version}" == 1.[0-6].* ]]; then
        if [[ "${version}" == 1.[0-4].* ]]; then
            # remove broken dylib placeholder to force static linking
            rm -f ${prefix}/lib/libosxunwind.dylib
        fi
        cat << EOM >>Make.user
        USECLANG=1

        # link against libosxunwind, see https://github.com/JuliaPackaging/Yggdrasil/pull/2164
        # and https://github.com/JuliaPackaging/Yggdrasil/pull/2190
        LIBUNWIND:=-losxunwind
        JCPPFLAGS+=-DLIBOSXUNWIND
    EOM
    fi

    # avoid linker errors related to atomic support in 32bit ARM builds
    if [[ "${bb_full_target}" == armv7l-* ]]; then
        echo "MARCH=armv7-a" >>Make.user
    fi

    # Add file to one of the `STD_LIB_PATH`
    if [[ "${target}" == *mingw* ]]; then
        cp /opt/*-w64-mingw32/*-w64-mingw32/sys-root/bin/libwinpthread-1.dll /opt/*-w64-mingw32/*-mingw32/sys-root/lib/
    fi

    make BUILDING_HOST_TOOLS=1 NO_GIT=1 -j${nproc} VERBOSE=1 -C src/flisp host/flisp
    make clean -C src
    make clean -C src/support
    make clean -C src/flisp

    # We don't trust the system libm in places
    # So we include a private copy of libopenlibm
    mkdir -p usr/lib
    cp ${prefix}/*/libopenlibm.a usr/lib/

    # Mac build complains about checksum
    rm -rf /workspace/srcdir/julia-1.5.1/deps/checksums/lapack-3.9.0.tgz

    # compile libjulia but don't try to build a sysimage
    if [[ "${version}" == 1.[0-5].* ]]; then
        MAKE_TARGET=julia-ui-release
    else
        MAKE_TARGET="julia-src-release julia-cli-release"
    fi
    make USE_CROSS_FLISP=1 NO_GIT=1 LDFLAGS="${LDFLAGS}" CFLAGS="${CFLAGS}" CXXFLAGS="${CXXFLAGS}" -j${nproc} VERBOSE=1 ${MAKE_TARGET}

    # 'manually' install libraries and headers
    mkdir -p ${libdir}
    mkdir -p ${includedir}/julia
    if [[ "${target}" == *mingw* ]]; then
        cp -r usr/bin/libjulia* ${bindir}/
    else
        cp -r usr/lib/libjulia* ${libdir}/
    fi

    cp -R -L usr/include/julia/* ${includedir}/julia
    install_license LICENSE.md
    """

    # HACK: JLLs are not allowed to use prerelease versions, so strip that out
    # e.g. for 1.7.0-beta2
    version = VersionNumber(version.major, version.minor, version.patch)

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = libjulia_platforms(version)

    # While the "official" Julia kernel ABI itself does not involve any C++
    # symbols on the linker level, `libjulia` still exports "unofficial" symbols
    # dependent on the C++ strings ABI (coming from LLVM related code). This
    # doesn't matter if the client code is pure C, but as soon as there are
    # other (actual) C++ dependencies, we must make sure to use the matching C++
    # strings ABI. Hence we must use `expand_cxxstring_abis` below.
    platforms = expand_cxxstring_abis(platforms)

    for p in platforms
        p["julia_version"] = string(version)
    end

    # The products that we will ensure are always built
    products = [
        LibraryProduct("libjulia", :libjulia; dont_dlopen=true),
    ]

    # Dependencies that must be installed before this package can be built/used

    dependencies = BinaryBuilder.AbstractDependency[
        Dependency("LibUnwind_jll"),
        BuildDependency("OpenLibm_jll"),
        BuildDependency("dSFMT_jll"),
        BuildDependency("LibUV_jll"),
        BuildDependency("utf8proc_jll"),
        BuildDependency("MbedTLS_jll"),
        BuildDependency("LibSSH2_jll"),
        BuildDependency("LibCURL_jll"),
        BuildDependency("Zlib_jll"),
        BuildDependency("p7zip_jll"),
        BuildDependency("MPFR_jll"),
        BuildDependency("GMP_jll"),
        BuildDependency("Objconv_jll"),
    ]
    if version < v"1.5.1"
        push!(dependencies, Dependency("LibOSXUnwind_jll", compat="0.0.5"))
    elseif version < v"1.7"
        push!(dependencies, Dependency("LibOSXUnwind_jll", compat="0.0.6"))
    end

    if version < v"1.6"
        push!(dependencies, BuildDependency("SuiteSparse_jll", compat="5.4.0"))
    else
        push!(dependencies, BuildDependency("SuiteSparse_jll"))
    end

    if version < v"1.7"
        push!(dependencies, BuildDependency("PCRE2_jll", compat="10.31"))
    #else
    #    push!(dependencies, BuildDependency("PCRE2_jll", compat="10.36"))
    end

    if version.major == 1 && version.minor == 3
        push!(dependencies, BuildDependency("OpenBLAS_jll", compat="0.3.5"))
        # there is no libLLVM_jll 6.0.1, so we use LLVM_jll instead
        push!(dependencies, Dependency("LLVM_jll", compat="6.0.1"))
        push!(dependencies, BuildDependency("LibGit2_jll", compat="0.28.2"))
    elseif version.major == 1 && version.minor == 4
        push!(dependencies, BuildDependency("OpenBLAS_jll", compat="0.3.5"))
        push!(dependencies, Dependency("libLLVM_jll", compat="8.0.1"))
        push!(dependencies, BuildDependency("LibGit2_jll", compat="0.28.2"))
    elseif version.major == 1 && version.minor == 5
        push!(dependencies, BuildDependency("OpenBLAS_jll", compat="0.3.9"))
        push!(dependencies, Dependency("libLLVM_jll", compat="9.0.1"))
        push!(dependencies, BuildDependency("LibGit2_jll", compat="0.28.2"))
    elseif version.major == 1 && version.minor == 6
        push!(dependencies, BuildDependency("OpenBLAS_jll", compat="0.3.10"))
        push!(dependencies, Dependency("libLLVM_jll", compat="11.0.0"))
        push!(dependencies, BuildDependency("LibGit2_jll", compat="1.0.1"))
    elseif version.major == 1 && version.minor == 7
        #push!(dependencies, BuildDependency("OpenBLAS_jll", compat="0.3.13"))
        #push!(dependencies, Dependency("libLLVM_jll", compat="12.0.0"))
        #push!(dependencies, BuildDependency("LibGit2_jll", compat="1.0.1"))

        # HACK: we can't install LLVM 12 JLLs for Julia 1.7 from within Julia 1.6. Similar
        # for several other standard JLLs.
        # So we use get_addable_spec below to "fake it" for now.
        # This means the resulting package has fewer dependencies declared, but at least it
        # will work and allow people to build JLL binaries ready for Julia 1.7
        push!(dependencies, BuildDependency(get_addable_spec("LLVM_full_jll", v"12.0.0+0")))

        # starting with Julia 1.7, we need LLVMLibUnwind_jll
        push!(dependencies, BuildDependency(get_addable_spec("LLVMLibUnwind_jll", v"11.0.1+1")))
    else
        error("Unsupported Julia version")
    end

    julia_compat = version ≥ v"1.7" ? "1.6" : "1.0"

    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
                   preferred_gcc_version=v"7", lock_microarchitecture=false, julia_compat)
end

max_stars_repo_name: galabovaa/Yggdrasil

max_issues_repo_name: galabovaa/Yggdrasil

max_forks_repo_name: galabovaa/Yggdrasil


text:
struct Orbital
    name  ::Symbol
    size::Int
    l   ::Int
end
const orbitals = [
    Orbital(:s, 1, 0),
    Orbital(:p, 3, 1),
    Orbital(:d, 5, 2),
    Orbital(:f, 7, 3),
    Orbital(:dz2, 1, 2),
    Orbital(:dxz, 1, 2),
    Orbital(:dyz, 1, 2),
    Orbital(Symbol("dx2-y2"), 1, 2),
    Orbital(:dxy, 1, 2),
]
orbital(s::Symbol) = getfirst(x -> x.name == s, orbitals)

Base.convert(::Type{Symbol}, x::Orbital) = x.name
orbsize(orb::Orbital) = orb.size
orbsize(orb::Symbol)  = orbsize(orbital(orb))

@with_kw struct Projection
    orb   ::Orbital = orbitals[1]
    start ::Int = 0
    last  ::Int = 0
end

orbital(proj::Projection) = proj.orb
orbsize(proj::Projection) = proj.last - proj.start + 1

"""
Adds projections to atoms.
"""
function addprojections!(atoms, projections_)
    t_start = 1
    for (proj_at, projs) in projections_
        for at in atoms
            if length(projs) <= length(projections(at))
                continue
            end
            for proj in projs
                if name(at) == proj_at
                    size   = orbsize(proj)
                    t_proj = Projection(orbital(proj), t_start, t_start + size - 1)
                    push!(projections(at), t_proj)
                    t_start += size
                end
            end
        end
    end
end

Base.zero(::Type{Projection}) = Projection(Orbital(:zero, 0, 0), 0, 0)
Base.zero(::Projection) = Projection(Orbital(:zero, 0, 0), 0, 0)

Base.range(proj::Projection)  = proj.start:proj.last

ext: jl


text:
function LinearAlgebra.det(M::Matrix{<:AbstractPolynomialLike})
    m = size(M)[1]
    if m > 2
        return sum((-1)^(i-1) * M[i,1] * LinearAlgebra.det(M[1:end .!= i, 2:end]) for i in 1:m)
    else
        return M[1,1] * M[2,2] - M[2,1] * M[1,2]
    end
end


text:
module LifeContingencies

using ActuaryUtilities
using MortalityTables
using Transducers
using Dates
using Yields

const mt = MortalityTables

export LifeContingency,
    Insurance, AnnuityDue, AnnuityImmediate,
    APV,
    SingleLife, Frasier, JointLife,
    LastSurvivor,
    survival,
    reserve_premium_net,
    insurance,
    annuity_due,
    annuity_immediate,
    premium_net,
    omega,
    survival,
    discount,
    benefit,
    probability,
    cashflows,
    cashflows,
    timepoints,
    present_value

# 'actuarial objects' that combine multiple forms of decrements (lapse, interest, death, etc)
abstract type Life end

"""
    struct SingleLife
        mort
        issue_age::Int
        alive::Bool
        fractional_assump::MortalityTables.DeathDistribution
    end

A `Life` object containing the necessary assumptions for contingent maths related to a single life. Use with a `LifeContingency` to do many actuarial present value calculations.

Keyword arguments:
- `mort` pass a mortality vector, which is an array of applicable mortality rates indexed by attained age
- `issue_age` is the assumed issue age for the `SingleLife` and is the basis of many contingency calculations.
- `alive` Default value is `true`. Useful for joint insurances with different status on the lives insured.
- `fractional_assump`. Default value is `Uniform()`. This is a `DeathDistribution` from the `MortalityTables.jl` package and is the assumption to use for non-integer ages/times.

# Examples
    using MortalityTables
    mort = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")

    SingleLife(
        mort       = mort.select[30],
        issue_age  = 30
    )
"""
struct SingleLife <: Life
    mort
    issue_age
    alive
    fractional_assump
end

function SingleLife(;mort,issue_age=nothing,alive=true,fractional_assump = mt.Uniform())
    return SingleLife(mort;issue_age,alive,fractional_assump)
end

function SingleLife(mort;issue_age=nothing,alive=true,fractional_assump = mt.Uniform())
    if isnothing(issue_age)
        issue_age = firstindex(mort)
    end

    if !(eltype(mort) <: Real)
        # most likely case is that mort is an array of vectors
        # use issue age to select the right one (assuming indexed with issue age
        return SingleLife(mort[issue_age],issue_age,alive,fractional_assump)
    else
        return SingleLife(mort,issue_age,alive,fractional_assump)
    end

end

"""
    JointAssumption()

An abstract type representing the different assumed relationship between the survival of the lives on a JointLife. Available options to use include:
- `Frasier()`
"""
abstract type JointAssumption end

"""
    Frasier()

The assumption of independnt lives in a joint life calculation.
Is a subtype of `JointAssumption`.
"""
struct Frasier <: JointAssumption end

"""
    Contingency()

An abstract type representing the different triggers for contingent benefits. Available options to use include:
- `LastSurvivor()`
"""
abstract type Contingency end

"""
    LastSurvivor()
The contingency whereupon benefits are payable upon both lives passing.
Is a subtype of `Contingency`
"""
struct LastSurvivor <: Contingency end

# TODO: Not Implemented
# """
#     FirstToDie()
# The contingency whereupon benefits are payable upon the first life passing.

# Is a subtype of `Contingency`
# """
# struct FirstToDie <: Contingency end

"""
    struct JointLife
        lives
        contingency
        joint_assumption
    end

    A `Life` object containing the necessary assumptions for contingent maths related to a joint life insurance. Use with a `LifeContingency` to do many actuarial present value calculations.

Keyword arguments:
- `lives` is a tuple of two `SingleLife`s
- `contingency` default is `LastSurvivor()`. It is the trigger for contingent benefits. See `?Contingency`.
- `joint_assumption` Default value is `Frasier()`. It is the assumed relationship between the mortality of the two lives. See `?JointAssumption`.

# Examples
    using MortalityTables
    mort = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")

    l1 = SingleLife(
        mort       = mort.select[30],
        issue_age  = 30
    )
    l2 = SingleLife(
        mort       = mort.select[30],
        issue_age  = 30
    )

    jl = JointLife(
        lives = (l1,l2),
        contingency = LastSurvivor(),
        joint_assumption = Frasier()
    )
"""
Base.@kwdef struct JointLife <: Life
    lives::Tuple{SingleLife,SingleLife}
    contingency::Contingency = LastSurvivor()
    joint_assumption::JointAssumption = Frasier()
end

"""
    struct LifeContingency
        life::Life
"""
struct LifeContingency
    life::Life
    int
end

Base.broadcastable(lc::LifeContingency) = Ref(lc)

"""
    omega(lc::LifeContingency)
    omega(l::Life)
    omega(i::InterestRate)

# `Life`s and `LifeContingency`s

Returns the last defined time_period for both the interest rate and mortality table.
Note that this is *different* than calling `omega` on a `MortalityTable`, which will give you the last `attained_age`.

Example: if the `LifeContingency` has issue age 60, and the last defined attained age for the `MortalityTable` is 100, then `omega` of the `MortalityTable` will be `100` and `omega` of the
`LifeContingency` will be `40`.

# `InterestRate`s

The last period that the interest rate is defined for. Assumed to be infinite (`Inf`) for
    functional and constant interest rate types. Returns the `lastindex` of the vector if
    a vector type.
"""
function mt.omega(lc::LifeContingency)
    # if one of the omegas is infinity, that's a Float so we need
    # to narrow the type with Int
    return Int(omega(lc.life))
end

function mt.omega(l::SingleLife)
    return mt.omega(l.mort) - l.issue_age + 1
end

function mt.omega(l::JointLife)
    return minimum( omega.(l.lives) )
end

###################
## COMMUTATIONS ###
###################

"""
    D(lc::LifeContingency, to_time)

``D_x`` is a retrospective actuarial commutation function which is the product of the survival and discount factor.
"""
function D(lc::LifeContingency, to_time)
    return discount(lc.int, to_time) * survival(lc,to_time)
end

"""
    l(lc::LifeContingency, to_time)

``l_x`` is a retrospective actuarial commutation function which is the survival up to a certain point in time. By default, will have a unitary basis (ie `1.0`), but you can specify `basis` keyword argument to use something different (e.g. `1000` is common in the literature.)
"""
function l(lc::LifeContingency, to_time; basis=1.0)
    return survival(lc.life,to_time) * basis
end

"""
    C(lc::LifeContingency, to_time)

``C_x`` is a retrospective actuarial commutation function which is the product of the discount factor and the difference in `l` (``l_x``).
"""
function C(lc::LifeContingency, to_time)
    discount(lc.int, to_time+1) * (l(lc,to_time) - l(lc, to_time+1))

end

"""
    N(lc::LifeContingency, from_time)

``N_x`` is a prospective actuarial commutation function which is the sum of the `D` (``D_x``) values from the given time to the end of the mortality table.
"""
function N(lc::LifeContingency, from_time)
    range = from_time:(omega(lc)-1)
    return foldxt(+,Map(from_time->D(lc, from_time)), range)
end

"""
    M(lc::LifeContingency, from_time)

The ``M_x`` actuarial commutation function where the `from_time` argument is `x`.
Issue age is based on the issue_age in the LifeContingency `lc`.
"""
function M(lc::LifeContingency, from_time)
    range = from_time:omega(lc)-1
    return foldxt(+,Map(from_time->C(lc, from_time)), range)
end

E(lc::LifeContingency, t, x) = D(lc,x + t) / D(lc,x)

##################
### Insurances ###
##################

abstract type Insurance end

LifeContingency(ins::Insurance) = LifeContingency(ins.life,ins.int)

struct WholeLife <: Insurance
    life
    int
end

struct Term <: Insurance
    life
    int
    n
end

"""
    Insurance(lc::LifeContingency; n=nothing)
    Insurance(life,interest; n=nothing)

Life insurance with a term period of `n`. If `n` is `nothing`, then whole life insurance.

Issue age is based on the `issue_age` in the LifeContingency `lc`.

# Examples

```
ins = Insurance(
    SingleLife(mort = UltimateMortality([0.5,0.5]),issue_age = 0),
    Yields.Constant(0.05),
    n = 1
)
```
"""
Insurance(lc::LifeContingency; n=nothing) = Insurance(lc.life,lc.int;n)

function Insurance(lc,int;n=nothing)
    if isnothing(n)
        return WholeLife(lc,int)
    elseif n < 1
        return ZeroBenefit(lc,int)
    else
        Term(lc,int,n)
    end
end

struct Due end
struct Immediate end

struct Annuity <: Insurance
    life
    int
    payable
    n
    start_time
    certain
    frequency
end

struct ZeroBenefit <: Insurance
    life
    int
end

function ZeroBenefit(lc::LifeContingency)
    return ZeroBenefit(lc.life,lc.int)
end

"""
    AnnuityDue(lc::LifeContingency; n=nothing, start_time=0; certain=nothing,frequency=1)
    AnnuityDue(life, interest; n=nothing, start_time=0; certain=nothing,frequency=1)

Annuity due with the benefit period starting at `start_time` and ending after `n` periods with `frequency` payments per year of `1/frequency` amount and a `certain` period with non-contingent payments.

# Examples

```
ins = AnnuityDue(
    SingleLife(mort = UltimateMortality([0.5,0.5]),issue_age = 0),
    Yields.Constant(0.05),
    n = 1
)
```

"""
function AnnuityDue(life, int; n=nothing,start_time=0,certain=nothing,frequency=1)
    if ~isnothing(n) && n < 1
        return ZeroBenefit(life,int)
    else
        Annuity(life,int,Due(),n,start_time,certain,frequency)
    end
end

function AnnuityDue(lc::LifeContingency; n=nothing,start_time=0,certain=nothing,frequency=1)
    return AnnuityDue(lc.life,lc.int;n,start_time,certain,frequency)
end

"""
    AnnuityImmediate(lc::LifeContingency; n=nothing, start_time=0; certain=nothing,frequency=1)
    AnnuityImmediate(life, interest; n=nothing, start_time=0; certain=nothing,frequency=1)

Annuity immediate with the benefit period starting at `start_time` and ending after `n` periods with `frequency` payments per year of `1/frequency` amount and a `certain` period with non-contingent payments.

# Examples

```
ins = AnnuityImmediate(
    SingleLife(mort = UltimateMortality([0.5,0.5]),issue_age = 0),
    Yields.Constant(0.05),
    n = 1
)
```

"""
function AnnuityImmediate(life, int; n=nothing,start_time=0,certain=nothing,frequency=1)
    if ~isnothing(n) && n < 1
        return ZeroBenefit(life,int)
    else
        return Annuity(life,int,Immediate(),n,start_time,certain,frequency)
    end
end

function AnnuityImmediate(lc::LifeContingency; n=nothing,start_time=0,certain=nothing,frequency=1)
    return AnnuityImmediate(lc.life,lc.int;n,start_time,certain,frequency)
end

"""
    survival(Insurance)

The survorship vector for the given insurance.
"""
function MortalityTables.survival(ins::Insurance)
    return [survival(ins.life,t-1) for t in timepoints(ins)]
end

function MortalityTables.survival(ins::Annuity)
    return [survival(ins.life,t) for t in timepoints(ins)]
end

"""
    discount(Insurance)

The discount vector for the given insurance.
"""
function Yields.discount(ins::Insurance)
    return Yields.discount.(ins.int,timepoints(ins))
end

"""
    benefit(Insurance)

The unit benefit vector for the given insurance.
"""
function benefit(ins::Insurance)
    return ones(length(timepoints(ins)))
end

function benefit(ins::ZeroBenefit)
    return zeros(length(timepoints(ins)))
end

function benefit(ins::Annuity)
    return ones(length(timepoints(ins))) ./ ins.frequency
end

"""
    probability(Insurance)

The vector of contingent benefit probabilities for the given insurance.
"""
function probability(ins::Insurance)
    return [survival(ins.life,t-1) * decrement(ins.life,t-1,t) for t in timepoints(ins)]
end

function probability(ins::ZeroBenefit)
    return ones(length(timepoints(ins)))
end

function probability(ins::Annuity)
    if isnothing(ins.certain)
        return  [survival(ins.life,t) for t in timepoints(ins)]
    else
        return [t <= ins.certain + ins.start_time ? 1.0 : survival(ins.life,t) for t in timepoints(ins)]
    end
end

"""
    cashflows(Insurance)

The vector of decremented benefit cashflows for the given insurance.
"""
function cashflows(ins::Insurance)
   return probability(ins) .* benefit(ins)
end

"""
    timepoints(Insurance)

The vector of times corresponding to the cashflow vector for the given insurance.
"""
function timepoints(ins::Insurance)
    return collect(1:omega(ins.life))
end

function timepoints(ins::Term)
    return collect(1:min(omega(ins.life),ins.n))
end

function timepoints(ins::ZeroBenefit)
    return [0.]
end

function timepoints(ins::Annuity)
    return timepoints(ins,ins.payable)
end

function timepoints(ins::Annuity,payable::Due)
    if isnothing(ins.n)
        end_time = omega(ins.life)
    else
        end_time = ins.n + ins.start_time - 1 / ins.frequency
    end
    timestep = 1 / ins.frequency
    collect(ins.start_time:timestep:end_time)
end

function timepoints(ins::Annuity,payable::Immediate)
    if isnothing(ins.n)
        end_time = omega(ins.life)
    else
        end_time = ins.n + ins.start_time
    end
    timestep = 1 / ins.frequency
    end_time = max(ins.start_time + timestep,end_time) # return at least one timepoint to avoid returning empty array
    collect((ins.start_time + timestep):timestep:end_time)
end

"""
    present_value(Insurance)

The actuarial present value of the given insurance.
"""
function ActuaryUtilities.present_value(ins)
    return present_value(ins.int,cashflows(ins),timepoints(ins))
end

"""
    premium_net(lc::LifeContingency)
    premium_net(lc::LifeContingency,to_time)

The net premium for a whole life insurance (without second argument) or a term life insurance through `to_time`.

The net premium is based on 1 unit of insurance with the death benfit payable at the end of the year and assuming annual net premiums.
"""
premium_net(lc::LifeContingency) = A(lc) / ä(lc)
premium_net(lc::LifeContingency,to_time) = A(lc,to_time) / ä(lc,to_time)

"""
     reserve_premium_net(lc::LifeContingency,time)

The net premium reserve at the end of year `time`.
"""
function  reserve_premium_net(lc::LifeContingency, time)
    PVFB = A(lc) - A(lc,n=time)
    PVFP = premium_net(lc) * (ä(lc) - ä(lc,n=time))
    return (PVFB - PVFP) / APV(lc,time)
end

"""
    APV(lc::LifeContingency,to_time)

The **actuarial present value** which is the survival times the discount factor for the life contingency.
"""
function APV(lc::LifeContingency,to_time)
    return survival(lc,to_time) * discount(lc.int,to_time)
end

"""
    decrement(lc::LifeContingency,to_time)
    decrement(lc::LifeContingency,from_time,to_time)

Return the probablity of death for the given LifeContingency.
"""
mt.decrement(lc::LifeContingency,from_time,to_time) = 1 - survival(lc.life,from_time,to_time)

"""
    survival(lc::LifeContingency,from_time,to_time)
    survival(lc::LifeContingency,to_time)

Return the probablity of survival for the given LifeContingency.
"""
mt.survival(lc::LifeContingency,to_time) = survival(lc.life, 0, to_time)
mt.survival(lc::LifeContingency,from_time,to_time) = survival(lc.life, from_time, to_time)

mt.survival(l::SingleLife,to_time) = survival(l,0,to_time)
mt.survival(l::SingleLife,from_time,to_time) =survival(l.mort,l.issue_age + from_time,l.issue_age + to_time, l.fractional_assump)

"""
    surival(life)

Return a survival vector for the given life.
"""
function mt.survival(l::Life)
    ω =   omega(l)
    return [survival(l,t) for t in 0:ω]
end

mt.survival(l::JointLife,to_time) = survival(l::JointLife,0,to_time)
function mt.survival(l::JointLife,from_time,to_time)
    return survival(l.contingency,l.joint_assumption,l::JointLife,from_time,to_time)
end

function mt.survival(ins::LastSurvivor,assump::JointAssumption,l::JointLife,from_time,to_time)
    to_time == 0 && return 1.0

    l1,l2 = l.lives
    ₜpₓ = survival(l1.mort,l1.issue_age + from_time,l1.issue_age + to_time,l1.fractional_assump)
    ₜpᵧ = survival(l2.mort,l2.issue_age + from_time,l2.issue_age + to_time,l2.fractional_assump)
    return ₜpₓ + ₜpᵧ - ₜpₓ * ₜpᵧ
end

Yields.discount(lc::LifeContingency,t) = discount(lc.int,t)
Yields.discount(lc::LifeContingency,t1,t2) = discount(lc.int,t1,t2)

# function compostion with kwargs.
# https://stackoverflow.com/questions/64740010/how-to-alias-composite-function-with-keyword-arguments
⋄(f, g) = (x...; kw...)->f(g(x...; kw...))

# unexported aliases
const V = reserve_premium_net
const v = Yields.discount
const A = present_value ⋄ Insurance
const a = present_value ⋄ AnnuityImmediate
const ä = present_value ⋄ AnnuityDue
const P = premium_net
const ω = omega

end # module


text:
macro check(clfunc)
    quote
        local err::CL_int
        err = $clfunc
        if err != CL_SUCCESS
            throw(CLError(err))
        end
        err
    end
end

macro check_release(clfunc)
    quote
        local err::CL_int
        err = $clfunc
        if err != CL_SUCCESS
            error("release! $clfunc failed with code $(err[1]))")
        end
    end
end

#TODO: these are just stubs for future expanded versions
macro ocl_v1_1_only(ex)
    quote
        $(esc(ex))
    end
end

macro ocl_v1_2_only(ex)
    quote
        $(esc(ex))
    end
end

macro return_event(evt)
    quote
        try
            return Event($(esc(evt)), retain=false)
        catch err
            @check api.clReleaseEvent($(esc(evt)))
            throw(err)
        end
    end
end

macro return_nanny_event(evt, obj)
    quote
        try
            return NannyEvent($(esc(evt)), $(esc(obj)))
        catch err
            @check api.clReleaseEvent($(esc(evt)))
            throw(err)
        end
    end
end

macro int_info(what, cl_obj_id, cl_obj_info, ret_type)
    local clFunc = Symbol("clGet$(what)Info")
    quote
        local result = Array($(esc(ret_type)), 1)
        local err::CL_int
        err = $clFunc($(esc(cl_obj)), $(esc(cl_obj_info)),
                      sizeof($(esc(ret_type))), result, C_NULL)
        if err != CL_SUCCESS
            throw(CLError(err))
        end
        result[1]
    end
end

macro vec_info(what, arg1, arg2, res_vec)
    local clFunc = Symbol("api.clGet$(what)Info")
    quote
        local size = Array(Csize_t, 1)
        @check clFunc($arg1, $arg2, 0, C_NULL, size)
        local n = size / sizeof(res_vec[1])
        resize!(res_vec, n)
        @check clFunc($arg1, $arg2, empty($res_vec) ? C_NULL : res_vec, size)
    end
end

macro str_info(what, arg1, arg2)
    local clFunc = Symbol("api.clGet$(what)Info")
    quote
        local size = Ref{Csize_t}()
        @check $(esc(clFunc))($(esc(arg1)), $(esc(arg2)), 0, C_NULL, size)
        local result = Array(CL_char, size[])
        @check $(esc(clFunc))($(esc(arg1)), $(esc(arg2)), size[], result, size)
        CLString(result)
    end
end

function _version_test(qm, elem :: Symbol, ex :: Expr, version :: VersionNumber)
    @assert qm == :?
    @assert ex.head == :(:)
    @assert length(ex.args) == 2

    esc(quote
        if cl.check_version($elem, $version)
            $(ex.args[1])
        else
            $(ex.args[2])
        end
    end)
end

macro min_v11(qm, elem, ex)
    _version_test(qm, elem, ex, v"1.1")
end

macro min_v12(qm, elem, ex)
    _version_test(qm, elem, ex, v"1.2")
end

macro min_v20(qm, elem, ex)
    _version_test(qm, elem, ex, v"2.0")
end


text:
import Mongoc

using Test
using Dates

@testset "BSON" begin

    @testset "as_json" begin
        @test_throws ErrorException Mongoc.BSON(""" { ajskdla sdjsafd } """)
        bson = Mongoc.BSON("""{"hey" : 1}""")
        @test bson["hey"] == 1
        @test Mongoc.as_json(bson) == """{ "hey" : 1 }"""
        @test Mongoc.as_json(bson, canonical=true) == """{ "hey" : { "\$numberInt" : "1" } }"""
    end

    @testset "oid conversion" begin
        local oid::Mongoc.BSONObjectId = "5b9fb22b3192e3fa155693a1"
        local oid_str::String = oid
    end

    @testset "oid compare" begin
        x = Mongoc.BSONObjectId()
        y = Mongoc.BSONObjectId()
        @test x != y
        @test x == x
        @test y == y
        @test hash(x) == hash(x)
        @test hash(y) == hash(y)
        @test hash(x) != hash(y)
        @test Mongoc.bson_oid_compare(x, y) < 0
        @test Mongoc.bson_oid_compare(y, x) > 0
        @test Mongoc.bson_oid_compare(x, x) == 0
        @test Mongoc.bson_oid_compare(y, y) == 0
    end

    @testset "oid string" begin
        x = Mongoc.BSONObjectId()
        y = Mongoc.BSONObjectId(string(x))
        @test x == y

        @test_throws ErrorException Mongoc.BSONObjectId("invalid_objectid")
    end

    @testset "oid time" begin
        Mongoc.get_time( Mongoc.BSONObjectId("5b9eaa2711c3dd0d6a46a5c4") ) == DateTime(2018, 9, 16, 19, 8, 23) # 2018-09-16T19:08:23
    end

    @testset "BSONObjectId segfault issue (#2)" begin
        io = IOBuffer()
        v = Vector{Mongoc.BSONObjectId}()

        for i in 1:5
            push!(v, Mongoc.BSONObjectId())
        end
        show(io, v)
    end

    @testset "oid sort" begin
        v = Vector{Mongoc.BSONObjectId}()
        for i in 1:10_000
            push!(v, Mongoc.BSONObjectId())
        end
        @test length(v) == length(unique(v))

        v_sorted = sort(v, lt = (a,b) -> Mongoc.bson_oid_compare(a,b) < 0)
        @test length(v_sorted) == length(v)
        for i in 1:length(v_sorted)
            @test v_sorted[i] == v[i]
        end
    end

    @testset "BSON Iterator" begin
        doc = Mongoc.BSON("""{ "a" : 1, "b" : 2.2, "str" : "my string", "bool_t" : true, "bool_f" : false, "array" : [1, 2, false, "inner_string"], "document" : { "a" : 1, "b" : "b_string"}, "null" : null  }""")

        new_id = Mongoc.BSONObjectId()
        doc["_id"] = new_id

        @test haskey(doc, "a")
        @test haskey(doc, "b")
        @test !haskey(doc, "c")
        @test haskey(doc, "str")
        @test haskey(doc, "_id")
        @test haskey(doc, "array")
        @test haskey(doc, "document")
        @test haskey(doc, "null")

        @test doc["a"] == 1
        @test doc["b"] == 2.2
        @test doc["str"] == "my string"
        @test doc["bool_t"]
        @test !doc["bool_f"]
        @test doc["_id"] == new_id
        @test doc["array"] == [1, 2, false, "inner_string"]
        @test doc["document"] == Dict("a"=>1, "b"=>"b_string")
        @test doc["null"] == nothing

        @test_throws ErrorException doc["invalid key"]

        doc_dict = Mongoc.as_dict(doc)
        @test doc_dict["a"] == 1
        @test doc_dict["b"] == 2.2
        @test doc_dict["str"] == "my string"
        @test doc_dict["bool_t"]
        @test !doc_dict["bool_f"]
        @test doc_dict["_id"] == new_id
        @test doc_dict["array"] == [1, 2, false, "inner_string"]
        @test doc_dict["document"] == Dict("a"=>1, "b"=>"b_string")
        @test doc_dict["null"] == nothing
    end

    @testset "BSON write" begin
        bson = Mongoc.BSON()
        new_oid = Mongoc.BSONObjectId()
        bson["_id"] = new_oid
        bson["int32"] = Int32(10)
        bson["int64"] = Int64(20)
        bson["string"] = "hey you"
        bson["bool_true"] = true
        bson["bool_false"] = false
        bson["double"] = 2.3
        bson["datetime"] = DateTime(2018, 2, 1, 10, 20, 35, 10)
        bson["vector"] = collect(1:10)
        bson["source"] = Mongoc.BSONCode("function() = 1")
        bson["null"] = nothing
        @test_throws ErrorException bson["key"] = Date(2018, 9, 18)

        let
            sub_bson = Mongoc.BSON()
            sub_bson["hey"] = "you"
            sub_bson["num"] = 10
            bson["sub_document"] = sub_bson
        end

        @test bson["_id"] == new_oid
        @test bson["int32"] == Int32(10)
        @test bson["int64"] == Int64(20)
        @test bson["string"] == "hey you"
        @test bson["bool_true"]
        @test !bson["bool_false"]
        @test bson["double"] == 2.3
        @test bson["datetime"] == DateTime(2018, 2, 1, 10, 20, 35, 10)
        @test bson["sub_document"]["hey"] == "you"
        @test bson["sub_document"]["num"] == 10
        @test bson["vector"] == collect(1:10)
        @test bson["source"] == Mongoc.BSONCode("function() = 1")
        @test bson["null"] == nothing

        let
            sub_bson = bson["sub_document"]
            @test sub_bson["hey"] == "you"
            @test sub_bson["num"] == 10
        end
    end

    @testset "BSON to Dict conversion" begin
        dict = Dict("a" => 1, "b" => false, "c" => "string", "d" => nothing)
        doc = Mongoc.BSON(dict)
        @test doc["a"] == 1
        @test doc["b"] == false
        @test doc["c"] == "string"
        @test doc["d"] == nothing
        @test dict == Mongoc.as_dict(doc)
    end

    @testset "BSON copy" begin
        @testset "exclude one key" begin
            src = Mongoc.BSON()
            src["hey"] = "you"
            src["out"] = 1
            dst = Mongoc.BSON()
            Mongoc.bson_copy_to_excluding_noinit(src.handle, dst.handle, "out")
            @test !haskey(dst, "out")
            @test dst["hey"] == "you"
        end

        @testset "no exclude keys" begin
            src = Mongoc.BSON()
            src["hey"] = "you"
            src["out"] = 1
            dst = Mongoc.BSON()
            Mongoc.bson_copy_to_excluding_noinit(src.handle, dst.handle)
            @test Mongoc.as_dict(src) == Mongoc.as_dict(dst)
        end
    end

    @testset "BSON Write to IO" begin

        @testset "write_bson closure" begin
            io = IOBuffer()

            Mongoc.bson_writer(io, initial_buffer_capacity=1) do writer
                Mongoc.write_bson(writer) do bson
                    bson["id"] = 10
                    bson["str"] = "aa"
                end
            end

            @test io.data == [0x1d,0x00,0x00,0x00,0x12,0x69,0x64,0x00,0x0a,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x02,0x73,0x74,0x72,0x00,0x03,0x00,0x00,0x00,0x61,0x61,0x00,0x00,0x00,0x00,0x00]
        end

        @testset "BSON copy single doc" begin
            src = Mongoc.BSON()
            src["id"] = 10
            src["str"] = "aa"

            io = IOBuffer()
            Mongoc.write_bson(io, src)

            @test io.data == [0x1d,0x00,0x00,0x00,0x12,0x69,0x64,0x00,0x0a,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x02,0x73,0x74,0x72,0x00,0x03,0x00,0x00,0x00,0x61,0x61,0x00,0x00,0x00,0x00,0x00]
        end

        @testset "BSON copy doc list" begin
            list = Vector{Mongoc.BSON}()

            let
                src = Mongoc.BSON()
                src["id"] = 1
                src["name"] = "1st"
                push!(list, src)
            end

            let
                src = Mongoc.BSON()
                src["id"] = 2
                src["name"] = "2nd"
                push!(list, src)
            end

            io = IOBuffer()
            Mongoc.write_bson(io, list)

            seekstart(io)
            vec_bson = Mongoc.read_bson(io)
            @test length(vec_bson) == 2

            let
                fst_bson = vec_bson[1]
                @test fst_bson["id"] == 1
                @test fst_bson["name"] == "1st"
            end

            let
                sec_bson = vec_bson[2]
                @test sec_bson["id"] == 2
                @test sec_bson["name"] == "2nd"
            end
        end

        @testset "read BSON from data" begin
            vec_bson = Mongoc.read_bson([0x1d,0x00,0x00,0x00,0x12,0x69,0x64,0x00,0x0a,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x02,0x73,0x74,0x72,0x00,0x03,0x00,0x00,0x00,0x61,0x61,0x00,0x00,0x00,0x00,0x00])
            @test length(vec_bson) == 1
            bson = vec_bson[1]
            @test bson["id"] == 10
            @test bson["str"] == "aa"
        end

        @testset "read/write BSON to/from File" begin
            filepath = joinpath(@__DIR__, "data.bson")
            isfile(filepath) && rm(filepath)

            list = Vector{Mongoc.BSON}()

            let
                src = Mongoc.BSON()
                src["id"] = 1
                src["name"] = "1st"
                push!(list, src)
            end

            let
                src = Mongoc.BSON()
                src["id"] = 2
                src["name"] = "2nd"
                push!(list, src)
            end

            try
                open(filepath, "w") do io
                    Mongoc.write_bson(io, list)
                end

                @test isfile(filepath)

                list_from_file = Mongoc.read_bson(filepath)
                @test length(list_from_file) == 2

                let
                    fst_bson = list_from_file[1]
                    @test fst_bson["id"] == 1
                    @test fst_bson["name"] == "1st"
                end

                let
                    sec_bson = list_from_file[2]
                    @test sec_bson["id"] == 2
                    @test sec_bson["name"] == "2nd"
                end

            finally
                isfile(filepath) && rm(filepath)
            end
        end
    end
end


text:
module MacrosTests

using ReverseDiff, ForwardDiff, Test, StaticArrays
using ForwardDiff: Dual, Partials, partials

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

tp = InstructionTape()
x, a, b = rand(3)

############
# @forward #
############

f0(x) = 1. / (1. + exp(-x))
f0(a, b) = sqrt(a^2 + b^2)

ReverseDiff.@forward f1(x::T) where {T<:Real} = 1. / (1. + exp(-x))
ReverseDiff.@forward f1(a::A, b::B) where {A,B<:Real} = sqrt(a^2 + b^2)

ReverseDiff.@forward f2(x) = 1. / (1. + exp(-x))
ReverseDiff.@forward f2(a, b) = sqrt(a^2 + b^2)

ReverseDiff.@forward function f3(x::T) where T<:Real
    return 1. / (1. + exp(-x))
end

ReverseDiff.@forward function f3(a::A, b::B) where {A,B<:Real}
    return sqrt(a^2 + b^2)
end

ReverseDiff.@forward function f4(x)
    return 1. / (1. + exp(-x))
end

ReverseDiff.@forward function f4(a, b)
    return sqrt(a^2 + b^2)
end

function test_forward(f, x, tp)
    xt = ReverseDiff.TrackedReal(x, zero(x), tp)

    y = f(x)
    @test isempty(tp)

    yt = f(xt)
    @test yt == y
    dual = f(Dual(x, one(x)))
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === xt
    @test instruction.output === yt
    @test instruction.cache[] === partials(dual, 1)
    empty!(tp)
end

function test_forward(f, a, b, tp)
    at = ReverseDiff.TrackedReal(a, zero(a), tp)
    bt = ReverseDiff.TrackedReal(b, zero(b), tp)

    c = f(a, b)
    dual = f(Dual(a, one(a), zero(a)), Dual(b, zero(b), one(b)))
    @test isempty(tp)

    tc = f(at, b)
    @test tc == c
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === (at, b)
    @test instruction.output === tc
    @test instruction.cache[] === SVector(partials(dual, 1), partials(dual, 1))
    empty!(tp)

    tc = f(a, bt)
    @test tc == c
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === (a, bt)
    @test instruction.output === tc
    @test instruction.cache[] === SVector(partials(dual, 2), partials(dual, 2))
    empty!(tp)

    tc = f(at, bt)
    @test tc == c
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === (at, bt)
    @test instruction.output === tc
    @test instruction.cache[] === SVector(partials(dual)...)
    empty!(tp)
end

for f in (ReverseDiff.@forward(f0), f1, f2, f3, f4, ReverseDiff.@forward(-))
    test_println("@forward named functions", f)
    test_forward(f, x, tp)
    test_forward(f, a, b, tp)
end

ReverseDiff.@forward f5 = (x) -> 1. / (1. + exp(-x))
test_println("@forward anonymous functions", f5)
test_forward(f5, x, tp)

ReverseDiff.@forward f6 = (a, b) -> sqrt(a^2 + b^2)
test_println("@forward anonymous functions", f6)
test_forward(f6, a, b, tp)

#########
# @skip #
#########

g0 = f0

ReverseDiff.@skip g1(x::T) where {T<:Real} = 1. / (1. + exp(-x))
ReverseDiff.@skip g1(a::A, b::B) where {A,B<:Real} = sqrt(a^2 + b^2)

ReverseDiff.@skip g2(x) = 1. / (1. + exp(-x))
ReverseDiff.@skip g2(a, b) = sqrt(a^2 + b^2)

ReverseDiff.@skip function g3(x::T) where T<:Real
    return 1. / (1. + exp(-x))
end

ReverseDiff.@skip function g3(a::A, b::B) where {A,B<:Real}
    return sqrt(a^2 + b^2)
end

ReverseDiff.@skip function g4(x)
    return 1. / (1. + exp(-x))
end

ReverseDiff.@skip function g4(a, b)
    return sqrt(a^2 + b^2)
end

function test_skip(g, x, tp)
    xt = ReverseDiff.TrackedReal(x, zero(x), tp)

    y = g(x)
    @test isempty(tp)

    yt = g(xt)
    @test yt === y
    @test isempty(tp)
end

function test_skip(g, a, b, tp)
    at = ReverseDiff.TrackedReal(a, zero(a), tp)
    bt = ReverseDiff.TrackedReal(b, zero(b), tp)

    c = g(a, b)
    @test isempty(tp)

    tc = g(at, b)
    @test tc === c
    @test isempty(tp)

    tc = g(a, bt)
    @test tc === c
    @test isempty(tp)

    tc = g(at, bt)
    @test tc === c
    @test isempty(tp)
end

for g in (ReverseDiff.@skip(g0), g1, g2, g3, g4)
    test_println("@skip named functions", g)
    test_skip(g, x, tp)
    test_skip(g, a, b, tp)
end

ReverseDiff.@skip g5 = (x) -> 1. / (1. + exp(-x))
test_println("@skip anonymous functions", g5)
test_skip(g5, x, tp)

ReverseDiff.@skip g6 = (a, b) -> sqrt(a^2 + b^2)
test_println("@skip anonymous functions", g6)
test_skip(g6, a, b, tp)

#########
# @grad #
#########

using LinearAlgebra
using ReverseDiff: @grad, TrackedReal, TrackedVector, TrackedMatrix, TrackedArray

@testset "@grad macro" begin
    x = rand(3);
    A = rand(3, 3);
    A_x = [vec(A); x];
    global custom_grad_called

    f1(x) = dot(x, x)
    f1(x::TrackedVector) = ReverseDiff.track(f1, x)
    @grad function f1(x::AbstractVector)
        global custom_grad_called = true
        xv = ReverseDiff.value(x)
        dot(xv, xv), Δ -> (Δ * 2 * xv,)
    end

    custom_grad_called = false
    g1 = ReverseDiff.gradient(f1, x)
    g2 = ReverseDiff.gradient(x -> dot(x, x), x)
    @test g1 == g2
    @test custom_grad_called

    f2(A, x) = A * x
    f2(A, x::TrackedVector) = ReverseDiff.track(f2, A, x)
    f2(A::TrackedMatrix, x) = ReverseDiff.track(f2, A, x)
    f2(A::TrackedMatrix, x::TrackedVector) = ReverseDiff.track(f2, A, x)
    @grad function f2(A::AbstractMatrix, x::AbstractVector)
        global custom_grad_called = true
        Av = ReverseDiff.value(A)
        xv = ReverseDiff.value(x)
        Av * xv, Δ -> (Δ * xv', Av' * Δ)
    end

    custom_grad_called = false
    g1 = ReverseDiff.gradient(x -> sum(f2(A, x)), x)
    g2 = ReverseDiff.gradient(x -> sum(A * x), x)
    @test g1 == g2
    @test custom_grad_called

    custom_grad_called = false
    g1 = ReverseDiff.gradient(A -> sum(f2(A, x)), A)
    g2 = ReverseDiff.gradient(A -> sum(A * x), A)
    @test g1 == g2
    @test custom_grad_called

    custom_grad_called = false
    g1 = ReverseDiff.gradient(A_x -> sum(f2(reshape(A_x[1:9], 3, 3), A_x[10:end])), A_x)
    g2 = ReverseDiff.gradient(A_x -> sum(reshape(A_x[1:9], 3, 3) * A_x[10:end]), A_x)
    @test g1 == g2
    @test custom_grad_called

    f3(A; dims) = sum(A, dims = dims)
    f3(A::TrackedMatrix; dims) = ReverseDiff.track(f3, A; dims = dims)
    @grad function f3(A::AbstractMatrix; dims)
        global custom_grad_called = true
        Av = ReverseDiff.value(A)
        sum(Av, dims = dims), Δ -> (zero(Av) .+ Δ,)
    end
    custom_grad_called = false
    g1 = ReverseDiff.gradient(A -> sum(f3(A, dims = 1)), A)
    g2 = ReverseDiff.gradient(A -> sum(sum(A, dims = 1)), A)
    @test g1 == g2
    @test custom_grad_called

    f4(::typeof(log), A; dims) = sum(log, A, dims = dims)
    f4(::typeof(log), A::TrackedMatrix; dims) = ReverseDiff.track(f4, log, A; dims = dims)
    @grad function f4(::typeof(log), A::AbstractMatrix; dims)
        global custom_grad_called = true
        Av = ReverseDiff.value(A)
        sum(log, Av, dims = dims), Δ -> (nothing, 1 ./ Av .* Δ)
    end
    custom_grad_called = false
    g1 = ReverseDiff.gradient(A -> sum(f4(log, A, dims = 1)), A)
    g2 = ReverseDiff.gradient(A -> sum(sum(log, A, dims = 1)), A)
    @test g1 == g2
    @test custom_grad_called

    f5(x) = log(x)
    f5(x::TrackedReal) = ReverseDiff.track(f5, x)
    @grad function f5(x::Real)
        global custom_grad_called = true
        xv = ReverseDiff.value(x)
        log(xv), Δ -> (1 / xv * Δ,)
    end
    custom_grad_called = false
    g1 = ReverseDiff.gradient(x -> f5(x[1]) * f5(x[2]) + exp(x[3]), x)
    g2 = ReverseDiff.gradient(x -> log(x[1]) * log(x[2]) + exp(x[3]), x)
    @test g1 == g2
    @test custom_grad_called

    f6(x) = sum(x)
    f6(x::TrackedArray{<:AbstractFloat}) = ReverseDiff.track(f6, x)
    @grad function f6(x::TrackedArray{T}) where {T <: AbstractFloat}
        global custom_grad_called = true
        xv = ReverseDiff.value(x)
        sum(xv), Δ -> (one.(xv) .* Δ,)
    end

    custom_grad_called = false
    g1 = ReverseDiff.gradient(f6, x)
    g2 = ReverseDiff.gradient(sum, x)
    @test g1 == g2
    @test custom_grad_called

    x2 = round.(Int, x)
    custom_grad_called = false
    g1 = ReverseDiff.gradient(f6, x2)
    g2 = ReverseDiff.gradient(sum, x2)
    @test g1 == g2
    @test !custom_grad_called
    f6(x::TrackedArray) = ReverseDiff.track(f6, x)
    @test_throws MethodError ReverseDiff.gradient(f6, x2)

    f7(x...) = +(x...)
    f7(x::TrackedReal{<:AbstractFloat}...) = ReverseDiff.track(f7, x...)
    @grad function f7(x::TrackedReal{T}...) where {T <: AbstractFloat}
        global custom_grad_called = true
        xv = ReverseDiff.value.(x)
        +(xv...), Δ -> one.(xv) .* Δ
    end
    custom_grad_called = false
    g1 = ReverseDiff.gradient(x -> f7(x...), x)
    g2 = ReverseDiff.gradient(sum, x)
    @test g1 == g2
    @test custom_grad_called

    f8(A; kwargs...) = sum(A, kwargs...)
    f8(A::TrackedMatrix; kwargs...) = ReverseDiff.track(f8, A; kwargs...)
    @grad function f8(A::AbstractMatrix; kwargs...)
        global custom_grad_called = true
        Av = ReverseDiff.value(A)
        sum(Av; kwargs...), Δ -> (zero(Av) .+ Δ,)
    end
    custom_grad_called = false
    g1 = ReverseDiff.gradient(A -> sum(f8(A, dims = 1)), A)
    g2 = ReverseDiff.gradient(A -> sum(sum(A, dims = 1)), A)
    @test g1 == g2
    @test custom_grad_called
end

end # module


text:
using StatsBase

"""
Compute `log(exp(a)+exp(b))`.
"""
function addExp(a::Float64, b::Float64)
	if a > b
		return a + log(1 + exp(b - a))
	end
	return b + log(1 + exp(a - b))
end

"""
    quantileMeans(x::AbstractVector, n::Integer)

Compute the n-quantile means of a vector `x`, i.e. compute the values
which partition `x` into `n` subsets of nearly equal size and
return the means of each subset.
"""
function quantileMeans(x::AbstractVector, n::Integer)
    quantiles = nquantile(x, n)
    means = zeros(Float64, n)
    for j in 1:n
        if j != n
            mask = map(e->(quantiles[j]<=e<quantiles[j+1]), x)
        else
            mask = map(e->(quantiles[j]<=e), x)
        end
        means[j] = mean(x[mask])
    end
    return means
end

"""
    dataLikelihood!(spn::SumProductNetwork, x::AbstractMatrix) -> Float64

Computes the log-likelihood of the data `x` given the current SPN model.
"""
function dataLikelihood!(spn::SumProductNetwork, x::AbstractMatrix)
    llh = 0.0
    for i in 1:size(x,1)
        xi = vec(x[i,:])
        setInput!(spn, xi)
        llhi = eval!(spn)
        llh += llhi
    end
    return llh
end


text:
±(a::Number, b::Number) = (a-b, a+b)

function find_shore_id(h_sea, River::Profile, Terrace::Profile, TerracePhysics, id::NamedTuple)
    # unpack
    zr, zt, n = River.z, Terrace.z, Terrace.nnod
    # indices
    ir, it = 0, 0
    midpoint = n÷2
    mr, mt = id.terrace.i1, id.river.i1
    # We start searching from the center of the profile and start moving with ±1 steps to both sides
    for j = 1:midpoint

        idx = ifelse(
            j == 1,
            mr ± 0, # ± 0 to make it type stable
            mr ± j
        )

        if ir == 0
            for i in idx
                # check river
                @inbounds ir = ifelse(
                    zr[i] <= h_sea,
                    i,
                    0
                )
            end
        end

        idx = ifelse(
            j == 1,
            mt ± 0, # ± 0 to make it type stable
            mt ± j
        )

        if it == 0
            for i in idx
                # check terrace
                @inbounds it = ifelse(
                    zt[i] <= h_sea,
                    i,
                    0
                )
            end
        end

        if (ir>0) && (it>0)
            infected_nodes_terrace = infected_indices(
                Terrace, TerracePhysics, it
            )
            infected_nodes_river = infected_indices(
                River, TerracePhysics, it
            )
            return (terrace = infected_nodes_terrace, river = infected_nodes_river)
        end

    end

end

function infected_indices(
    profile::Profile, TerracePhysics::TerraceParameters, id_shore::Int64
    )

    idx, i1 = id_shore, id_shore
    z = profile.z
    h_wb = 3/2*TerracePhysics.h_wb # influence of sea bed erosion (x1.5 as a safety buffer zone)
    zshore = z[idx] # height of sea level
    zᵢ = z[idx+1]

    while (zshore - zᵢ) < h_wb
        idx += 1
        @inbounds zᵢ = z[idx]
    end

    (i1 = i1, i2 = idx)

end

function find_shore_id(h_sea::Float64,sol::Array{Float64},n::Int64)
    for ii = 1:n-1
         # @inbounds if sol[ii] > h_sea && sol[ii+1] <= h_sea
        @inbounds if sol[ii] <= h_sea
            return ii
            break
        end
    end
end

function find_shore_id_terrace(h_sea::Float64,sol::Array{Float64},n::Int64)
    i1=1
    for ii = 1:n-1
         # @inbounds if sol[ii] > h_sea && sol[ii+1] <= h_sea
        @inbounds if sol[ii] <= h_sea
            i1 =  ii
        end
        @inbounds if h_sea-sol[ii] > 500
            return i1, ii
        end
    end
end


text:
"""
Build D matrix used in LSQ solver algorithm.
"""
function _buildD(func, β, ζ, Σ)

    Nf = func(β, ζ) |> length
    Nβ = length(β)
    Nζ = length(ζ)

    ∇fβ, ∇fζ = jacobian(func, β, ζ)

    return [
        zeros(Nβ, Nβ) zeros(Nβ, Nζ)      transpose(∇fβ);
        zeros(Nζ, Nβ)        inv(Σ)      transpose(∇fζ);
                  ∇fβ           ∇fζ       zeros(Nf, Nf)
    ]
end

"""
Apply solver algorithm to model `func` and data `input`.
"""
function solve(func::Function, input::LSQInput; iter_max=100, tol=1e-8)

    βₗ, z, Σ = input.β, input.z, input.Σ
    ζₗ = copy(z)
    Nf = func(βₗ, ζₗ) |> length
    Nβ = length(βₗ)
    Nζ = length(ζₗ)

    vₗ = [βₗ; ζₗ; zeros(Nf)] # iteratively updated vector

    # previous chi-square
    χ²ₗ₋₁ = χ² = nothing

    for l in 1:iter_max
        D = _buildD(func, βₗ, ζₗ, Σ)
        f = func(βₗ, ζₗ)
        vₗ += D \ [zeros(Nβ); Σ\(z - ζₗ); -f]
        βₗ = vₗ[1:Nβ]
        ζₗ = vₗ[Nβ+1:Nβ+Nζ]
        χ² = transpose(z - ζₗ)*inv(Σ)*(z - ζₗ)
        println("iteration $l: χ² = $(χ²)")
        if !isnothing(χ²ₗ₋₁) && abs((χ²ₗ₋₁ - χ²)/χ²) < tol
            break
        else
            χ²ₗ₋₁ = χ²
        end
    end

    invD = inv(_buildD(func, βₗ, ζₗ, Σ))
    Σβ = invD[1:Nβ, 1:Nβ]
    Σζ = invD[Nβ+1:Nβ+Nζ, Nβ+1:Nβ+Nζ]
    Σβζ= invD[1:Nβ, Nβ+1:Nβ+Nζ]

    return LSQOutput(βₗ, ζₗ, Σβ, Σζ, Σβζ, χ², Nf - Nβ, input.βindex, input.zindex)
end


text:
function __init__()

# THESE ARE THE THINGS TO MERGE!
# register(DataDep("VarianceComponentsHDFE",
#        "Downloading medium_nocontrols_pruned.csv, medium_controls_pruned.csv, large_nocontrols_pruned.csv, and large_controls_pruned.csv, huge_pruned_main.csv, large_controls_main.csv",
#        ["https://vchdfe.s3-us-west-2.amazonaws.com/medium_pruned_main.tar.gz",
#        "https://vchdfe.s3-us-west-2.amazonaws.com/medium_controls_pruned_main.tar.gz",
#        "https://vchdfe.s3-us-west-2.amazonaws.com/large_pruned_main.tar.gz",
#        "https://vchdfe.s3-us-west-2.amazonaws.com/large_controls_pruned_main.tar.gz",
#        "https://vchdfe.s3-us-west-2.amazonaws.com/huge_pruned_main.tar.gz",
#        "https://vchdfe.s3-us-west-2.amazonaws.com/large_controls_main.tar.gz"],
#        ["e01474ebd8d7df6d645b0c865079a108c9fd33e052b7ab163bb2f5749bc58387",
#        "a896b114eb952e60212196ecd1dffe2ea8bded263590a5b93e3b30b069feb0e1",
#        "583ec55959b16ef93db210736c93daa8c8eab9097d4492347fe442237f6b3c2c",
#        "e6ca435c7ac3bf7a8e425ab58581f36877845c099ad4da086d1990a3c270ac88",
#        "5743ed199dd6e6180f6fd721d3232bec0573a05d33f358c9a228a930e517ec12",
#        "82a95b368b2f6be381f10228026ebf2ffffc7e0d80bf5580ee905adf29b26128"];
#        post_fetch_method= [unpack, unpack, unpack, unpack, unpack, unpack]
# ))
    return nothing
end


text:
export DenseConvDims

"""
    DenseConvDims

Concrete subclass of `ConvDims` for a normal, dense, conv2d/conv3d.
"""
struct DenseConvDims{N, K, S, P, D} <: ConvDims{N}
    input_size::NTuple{N, Int}

    kernel_size::NTuple{K, Int}
    channels_in::Int
    channels_out::Int
    groupcount::Int

    stride::NTuple{S, Int}
    padding::NTuple{P, Int}
    dilation::NTuple{D, Int}
    flipkernel::Bool
end

function DenseConvDims(
    x_size::NTuple{M}, w_size::NTuple{M};
    stride = 1, padding = 0, dilation = 1, groups = 1,
    flipkernel::Bool = false,
) where {M}
    sstride, ppadding, ddilation = check_spdf(
        x_size, w_size, stride, padding, dilation)

    # Ensure channels are equal
    if x_size[end - 1] != w_size[end - 1] * groups
        xs = x_size[end - 1]
        ws = w_size[end - 1]
        throw(DimensionMismatch("Input channels must match! ($xs vs. $ws)"))
    end

    # Ensure groups are valid
    if x_size[end - 1] % w_size[end - 1] != 0 || w_size[end] % groups != 0
        throw(DimensionMismatch(
            "Group count should be divisble by input and output channels ($groups vs. $(w_size[end-1:end]))"))
    end

    DenseConvDims(
        x_size[1:(end - 2)],
        w_size[1:(end - 2)], x_size[end - 1], w_size[end], groups,
        sstride, ppadding, ddilation, flipkernel)
end

function DenseConvDims(x::AbstractArray, w::AbstractArray; kwargs...)
    if ndims(x) != ndims(w)
        throw(DimensionMismatch(
            "Rank of x and w must match! ($(ndims(x)) vs. $(ndims(w)))"))
    end
    return DenseConvDims(size(x), size(w); kwargs...)
end

# Useful for constructing a new DenseConvDims that has only a few elements different
# from the original progenitor object that it inherits shapes from.
@inline DenseConvDims(
    c::C; I=input_size(c), K=kernel_size(c),
    C_in=channels_in(c), C_out=channels_out(c), S=stride(c),
    P=padding(c), D=dilation(c), F=flipkernel(c), G=groupcount(c),
) where C <: ConvDims = DenseConvDims(
    I,
    K, C_in, C_out, G,
    S, P, D, F)

@inline groupcount(c::DenseConvDims) = c.groupcount
@inline channels_in(c::DenseConvDims) = c.channels_in
@inline channels_out(c::DenseConvDims) = c.channels_out

@inline input_size(c::DenseConvDims) = c.input_size
@inline kernel_size(c::DenseConvDims) = c.kernel_size

@inline stride(c::DenseConvDims) = c.stride
@inline padding(c::DenseConvDims) = c.padding
@inline dilation(c::DenseConvDims) = c.dilation
@inline flipkernel(c::DenseConvDims) = c.flipkernel

function check_dims(x::NTuple{M}, w::NTuple{M}, y::NTuple{M}, cdims::DenseConvDims) where {M}
    # First, check that channel counts are all correct:
    @assert x[M-1] * groupcount(cdims) == channels_in(cdims) DimensionMismatch("Data input channel count ($(x[M-1]) vs. $(channels_in(cdims)))")
    @assert y[M-1] == channels_out(cdims) ÷ groupcount(cdims)  DimensionMismatch("Data output channel count ($(y[M-1]) vs. $(channels_out(cdims)))")
    @assert w[M-1] * groupcount(cdims) == channels_in(cdims) DimensionMismatch("Kernel input channel count ($(w[M-1]) vs. $(channels_in(cdims)))")
    @assert w[M] * groupcount(cdims) == channels_out(cdims) DimensionMismatch("Kernel output channel count ($(w[M]) vs. $(channels_out(cdims)))")

    # Next, check that the spatial dimensions match up
    @assert x[1:M-2] == input_size(cdims) DimensionMismatch("Data input spatial size ($(x[1:M-2]) vs. $(input_size(cdims)))")
    @assert y[1:M-2] == output_size(cdims) DimensionMismatch("Data output spatial size ($(y[1:M-2]) vs. $(output_size(cdims)))")
    @assert w[1:M-2] == kernel_size(cdims) DimensionMismatch("Kernel spatial size ($(w[1:M-2]) vs. $(kernel_size(cdims)))")

    # Check the groups match
    @assert channels_in(cdims) % groupcount(cdims) == 0 DimensionMismatch("Groups ($(groupcount(cdims))) should be divisble by input channels $(channels_in(cdims))")

    # Finally, check that the batch size matches
    @assert x[M] == y[M] DimensionMismatch("Batch size ($(x[M]) vs. $(y[M]))")
end


text:
# Ref: https://github.com/JuliaSmoothOptimizers/Percival.jl/blob/master/src/AugLagModel.jl
# TODO:
# - check case when lb = ub for inequality constraints
# Two-sided Lagrangian
mutable struct AugLagEvaluator{T} <: AbstractNLPEvaluator
    inner::AbstractNLPEvaluator
    cons::AbstractVector{T}
    infeasibility::AbstractVector{T}
    ρ::T
    λ::AbstractVector{T}
    λc::AbstractVector{T}
    # Scaling
    scaler::AbstractScaler
    # Stats
    counter::AbstractCounter
end
function AugLagEvaluator(nlp::AbstractNLPEvaluator, u0;
                         penalties=Float64[],
                         scale=false,
                         c₀=0.1)

    if n_constraints(nlp) == 0
        @warn("Original model has no inequality constraint")
    end

    cons = similar(nlp.g_min)
    cx = similar(nlp.g_min)
    λc = similar(nlp.g_min) ; fill!(λc, 0)
    λ = similar(nlp.g_min) ; fill!(λ, 0)

    scaler = scale ?  MaxScaler(nlp, u0) : MaxScaler(nlp.g_min, nlp.g_max)
    return AugLagEvaluator(nlp, cons, cx, c₀, λ, λc, scaler, NLPCounter())
end

initial(ag::AugLagEvaluator) = initial(ag.inner)
bounds(ag::AugLagEvaluator, ::Variables) = ag.inner.u_min, ag.inner.u_max
bounds(ag::AugLagEvaluator, ::Constraints) = Float64[], Float64[]
n_variables(ag::AugLagEvaluator) = n_variables(ag.inner)
n_constraints(ag::AugLagEvaluator) = 0

function update!(ag::AugLagEvaluator, u)
    conv = update!(ag.inner, u)
    # Update constraints
    constraint!(ag.inner, ag.cons, u)
    # Rescale
    ag.cons .*= ag.scaler.scale_cons
    # Update (shifted) infeasibility error
    g♭ = ag.scaler.g_min
    g♯ = ag.scaler.g_max

    ag.λc = max.(0, ag.λ .+ ag.ρ .* (ag.cons .- g♯)) .+
            min.(0, ag.λ .+ ag.ρ .* (ag.cons .- g♭))
    ag.infeasibility .= max.(0, ag.cons .- g♯) .+ min.(0, ag.cons .- g♭)
    return conv
end

function update_penalty!(ag::AugLagEvaluator; η=10.0)
    ag.ρ = min(η * ag.ρ, 10e12)
end
function update_multipliers!(ag::AugLagEvaluator)
    ag.λ .= ag.λc
    return
end

function objective(ag::AugLagEvaluator, u)
    ag.counter.objective += 1
    base_nlp = ag.inner
    cx = ag.infeasibility
    # TODO: add multiplier
    obj = ag.scaler.scale_obj * objective(base_nlp, u) +
        0.5 * ag.ρ * dot(cx, cx) + dot(ag.λ, cx)
    return obj
end
function inner_objective(ag::AugLagEvaluator, u)
    return ag.scaler.scale_obj * objective(ag.inner, u)
end

function gradient!(ag::AugLagEvaluator, grad, u)
    ag.counter.gradient += 1
    base_nlp = ag.inner
    scaler = ag.scaler
    model = base_nlp.model
    λ = base_nlp.λ
    # Import buffer (has been updated previously in update!)
    buffer = base_nlp.buffer
    # Import AutoDiff objects
    ∇gᵤ = jacobian(model, base_nlp.autodiff.Jgᵤ, buffer)
    ∂obj = base_nlp.autodiff.∇f
    jvx = ∂obj.jvₓ ; fill!(jvx, 0)
    jvu = ∂obj.jvᵤ ; fill!(jvu, 0)

    # compute gradient of objective
    ∂cost(model, ∂obj, buffer)
    jvx .+= scaler.scale_obj .* ∂obj.∇fₓ
    jvu .+= scaler.scale_obj .* ∂obj.∇fᵤ

    # compute gradient of penalties
    fr_ = 0
    for cons in base_nlp.constraints
        n = size_constraint(base_nlp.model, cons)
        mask = fr_+1:fr_+n
        cx = @view ag.λc[mask]
        v = cx .* scaler.scale_cons[mask]
        jtprod(model, cons, ∂obj, buffer, v)
        jvx .+= ∂obj.∇fₓ
        jvu .+= ∂obj.∇fᵤ
        fr_ += n
    end
    # evaluate reduced gradient
    LinearSolvers.ldiv!(base_nlp.linear_solver, λ, base_nlp.∇gᵗ, jvx)
    grad .= jvu
    mul!(grad, transpose(∇gᵤ), λ, -1.0, 1.0)
end

function constraint!(ag::AugLagEvaluator, cons, u)
    @assert length(cons) == 0
    return
end

function jacobian!(ag::AugLagEvaluator, jac, u)
    @assert length(jac) == 0
    return
end

function jacobian_structure!(ag::AugLagEvaluator, rows, cols)
    @assert length(rows) == length(cols) == 0
end

primal_infeasibility!(ag::AugLagEvaluator, cons, u) = primal_infeasibility!(ag.inner, cons, u)
primal_infeasibility(ag::AugLagEvaluator, u) = primal_infeasibility(ag.inner, u)

function reset!(ag::AugLagEvaluator)
    reset!(ag.inner)
    empty!(ag.counter)
    fill!(ag.cons, 0)
    fill!(ag.infeasibility, 0)
    fill!(ag.λ, 0)
    fill!(ag.λc, 0)
end


text:
import Base: getindex, *, ndims, length, ^
import LinearAlgebra: ×

abstract type AbstractDomain end
abstract type AbstractBoundary end

struct DomainBoundary <: AbstractBoundary
    closure::Any
end

struct PointDomain{S} <: AbstractDomain
    point::S
end

struct IntervalDomain{AT, BT, PT} <: AbstractDomain
    min::AT
    max::BT
    periodic::PT
end

function IntervalDomain(min, max; periodic = false)
    @assert min < max
    return IntervalDomain(min, max, periodic)
end

function Periodic(min, max)
    @assert min < max
    return IntervalDomain(min, max, periodic = true)
end

S¹ = Periodic

function Interval(min, max)
    @assert min < max
    return IntervalDomain(min, max)
end

function Periodic(min, max)
    @assert min < max
    return IntervalDomain(min, max, periodic = true)
end

function Base.show(io::IO, Ω::IntervalDomain)
    min = Ω.min
    max = Ω.max
    printstyled(io, "[", color = 226)
    astring = @sprintf("%0.2f", min)
    bstring = @sprintf("%0.2f", max)
    printstyled(astring, ", ", bstring, color = 7)
    # printstyled("$min, $max", color = 7)
    Ω.periodic ? printstyled(io, ")", color = 226) :
    printstyled(io, "]", color = 226)
end

function Base.show(io::IO, o::PointDomain)
    printstyled("{", o.point, "}", color = 201)
end

# Product Domains
struct ProductDomain{DT} <: AbstractDomain
    domains::DT
end

function Base.show(io::IO, Ω::ProductDomain)
    for (i, domain) in enumerate(Ω.domains)
        print(domain)
        if i != length(Ω.domains)
            printstyled(io, "×", color = 118)
        end
    end
end

ndims(p::PointDomain) = 0
ndims(Ω::IntervalDomain) = 1
ndims(Ω::ProductDomain) = +(ndims.(Ω.domains)...)

length(Ω::IntervalDomain) = Ω.max - Ω.min
length(Ω::ProductDomain) = length.(Ω.domains)

×(arg1::AbstractDomain, arg2::AbstractDomain) = ProductDomain((arg1, arg2))
×(args::ProductDomain, arg2::AbstractDomain) =
    ProductDomain((args.domains..., arg2))
×(arg1::AbstractDomain, args::ProductDomain) =
    ProductDomain((arg1, args.domains...))
×(arg1::ProductDomain, args::ProductDomain) =
    ProductDomain((arg1.domains..., args.domains...))
×(args::AbstractDomain) = ProductDomain(args...)
*(arg1::AbstractDomain, arg2::AbstractDomain) = arg1 × arg2

function ^(Ω::IntervalDomain, T::Int)
    Ωᵀ = Ω
    for i in 1:(T - 1)
        Ωᵀ *= Ω
    end
    return Ωᵀ
end

function info(Ω::ProductDomain)
    println("This is a ", ndims(Ω), "-dimensional tensor product domain.")
    print("The domain is ")
    println(Ω, ".")
    for (i, domain) in enumerate(Ω.domains)
        domain_string = domain.periodic ? "periodic" : "wall-bounded"
        length = @sprintf("%.2f ", domain.max - domain.min)
        println(
            "The dimension $i domain is ",
            domain_string,
            " with length ≈ ",
            length,
        )
    end
    return nothing
end

function isperiodic(Ω::ProductDomain)
    max = [Ω.domains[i].periodic for i in eachindex(Ω.domains)]
    return prod(max)
end

function periodicityof(Ω::ProductDomain)
    periodicity = ones(Bool, ndims(Ω))
    for i in 1:ndims(Ω)
        periodicity[i] = Ω[i].periodic
    end
    return Tuple(periodicity)
end

getindex(Ω::ProductDomain, i::Int) = Ω.domains[i]

# Boundaries
struct Boundaries{S}
    boundaries::S
end

getindex(∂Ω::Boundaries, i) = ∂Ω.boundaries[i]

function Base.show(io::IO, ∂Ω::Boundaries)
    for (i, boundary) in enumerate(∂Ω.boundaries)
        printstyled("boundary ", i, ": ", color = 13)
        println(boundary)
    end
end

function ∂(Ω::IntervalDomain)
    if Ω.periodic
        return (nothing)
    else
        return Boundaries((PointDomain(Ω.min), PointDomain(Ω.max)))
    end
    return nothing
end

function ∂(Ω::ProductDomain)
    ∂Ω = []
    for domain in Ω.domains
        push!(∂Ω, ∂(domain))
    end
    splitb = []
    for (i, boundary) in enumerate(∂Ω)
        tmp = Any[]
        push!(tmp, Ω.domains...)
        if boundary != nothing
            tmp1 = copy(tmp)
            tmp2 = copy(tmp)
            tmp1[i] = boundary[1]
            push!(splitb, ProductDomain(Tuple(tmp1)))
            tmp2[i] = boundary[2]
            push!(splitb, ProductDomain(Tuple(tmp2)))
        end
    end
    return Boundaries(Tuple(splitb))
end

Base.@kwdef struct SphericalShellDomain{ℛ, ℋ, 𝒟} <: AbstractDomain
    radius::ℛ = 6378e3
    height::ℋ = 30e3
    depth::𝒟 = 3e3
end

length(Ω::SphericalShellDomain) = (Ω.radius, Ω.radius, Ω.height - Ω.depth)
ndims(Ω::SphericalShellDomain) = 3

function AtmosDomain(; radius = 6378e3, height = 30e3)
    depth = 0
    return SphericalShellDomain(; radius, height, depth)
end

function OceanDomain(; radius = 6378e3, depth = 3e3)
    height = 0
    return SphericalShellDomain(; radius, height, depth)
end


text:
struct BitsIntBool
    a::Int
    b::Bool
end

struct BitsCharFloat32Float64
    a::Char
    b::Float32
    c::Float64
end

@testset "Multi-field bits types" begin
    @test begin
        b = JlrsReflect.reflect([BitsIntBool])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsIntBool] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsIntBool")]
        pub struct BitsIntBool {
            pub a: i64,
            pub b: ::jlrs::wrappers::inline::bool::Bool,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsCharFloat32Float64])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsCharFloat32Float64] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsCharFloat32Float64")]
        pub struct BitsCharFloat32Float64 {
            pub a: ::jlrs::wrappers::inline::char::Char,
            pub b: f32,
            pub c: f64,
        }"""
    end
end


text:
using Test

supplyTemp = Supply("item00001", "each", "node001", "onhand", 100)
demandTemp = Demand("item00001", "each", "node001", "allocated", 10)

@testset "Model Tests" begin
    @testset "Supply model" begin
        @test supplyTemp.itemId == "item00001"
        @test supplyTemp.uom == "each"
        @test supplyTemp.nodeId == "node001"
        @test supplyTemp.type == "onhand"
        @test supplyTemp.quantity == 100
    end

    @testset "Demand model" begin
        @test demandTemp.itemId == "item00001"
        @test demandTemp.uom == "each"
        @test demandTemp.nodeId == "node001"
        @test demandTemp.type == "allocated"
        @test demandTemp.quantity == 10
    end
end

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
# There are two files in .gz with .data extension, need to pick the right one

dataset_path = joinpath(dirname(@__FILE__), "cpu-act.data.orig")
file_path = joinpath(dirname(@__FILE__), "cpu-act.data.gz-contents",
                     "cpu_act.data")

write("$dataset_path.custom", readstring(file_path))

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
@inline function allequal(x)
    length(x) < 2 && return true
    first_el = x[1]
    i = 2
    @inbounds for i=2:length(x)
        x[i] == first_el || return false
    end
    return true
end

"""
    vcat(A::SnpArray...;des="tmp_vcat_" * string(vcat_counter))

Concatenate SnpArray along dimension 1. (i.e. stack persons)
Only creates ".bed" file. Use vcat(A::SnpData...) for ".fam" and ".bim" files.
"""

arr_vcat_counter = 1
function Base.vcat(A::SnpArray...; des="tmp_vcat_arr_" * string(arr_vcat_counter))
    global arr_vcat_counter
    if des == "tmp_vcat_arr_" * string(arr_vcat_counter)
        arr_vcat_counter = arr_vcat_counter + 1
    end
    num_rows = [size(x)[1] for x in A]
    num_cols = [size(x)[2] for x in A]

    @assert allequal(num_cols) "number of columns are not the same"
    des_rows = sum(num_rows)
    des_cols = num_cols[1]
    cum_rows = Base.vcat([0], cumsum(num_rows))

    desbedfile = des * ".bed"

    makestream(desbedfile, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, Matrix{UInt8}(undef, (des_rows + 3) >> 2, des_cols))
    end
    des_arr = SnpArray(desbedfile, des_rows, "r+")

    for i in 1:length(A)
        des_arr[(cum_rows[i] + 1): cum_rows[i + 1], :] .= @view A[i][:, :]
    end
    des_arr
end

"""
    hcat(A...;des="tmp_hcat_" * string(hcat_counter))

Concatenate SnpData along dimension 2. (i.e. stack locations)
Only creates ".bed" file. Use hcat(A::SnpData...) for ".fam" and ".bim" files.
"""
arr_hcat_counter = 1
function Base.hcat(A::SnpArray...; des="tmp_hcat_arr_" * string(arr_hcat_counter))
    global arr_hcat_counter
    if des == "tmp_hcat_arr_"*string(arr_vcat_counter)
        arr_hcat_counter = arr_hcat_counter + 1
    end
    num_rows = [size(x)[1] for x in A]
    num_cols = [size(x)[2] for x in A]

    @assert allequal(num_rows) "number of rows are not the same"
    des_rows = num_rows[1]
    des_cols = sum(num_cols)
    cum_cols = Base.vcat([0], cumsum(num_cols))

    desbedfile = des * ".bed"

    makestream(desbedfile, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, Matrix{UInt8}(undef, (des_rows + 3) >> 2, des_cols))
    end
    des_arr = SnpArray(desbedfile, des_rows, "r+")

    des_arr.data .= hcat([s.data for s in A]...)
    des_arr
end

"""
   hvcat(rows::Tuple{Vararg{Int}}, values...; des="tmp_hvcat_arr_" * string(arr_hvcat_counter))

Horizontal and vertical concatenation in one call.
Only creates ".bed" file. Use hvcat(rows::Tuple{Vararg{Int}}, values::SnpData...) for ".fam" and ".bim" files.
"""
arr_hvcat_counter = 1
function Base.hvcat(rows::Tuple{Vararg{Int}}, A::SnpArray...; des="tmp_hvcat_arr_" * string(arr_hvcat_counter))
    global arr_hvcat_counter
    if des == "tmp_hvcat_arr_"*string(arr_hvcat_counter)
        arr_hvcat_counter = arr_hvcat_counter + 1
    end

    num_block_rows = length(rows)

    num_cols = 0
    for i = 1:rows[1]
        num_cols += size(A[i], 2)
    end

    num_rows = 0
    a = 1
    for i = 1:num_block_rows
        num_rows += size(A[a], 1)
        a += rows[i]
    end

    desbedfile = des * ".bed"

    makestream(desbedfile, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, Matrix{UInt8}(undef, (num_rows +3) >> 2, num_cols))
    end
    des_arr = SnpArray(desbedfile, num_rows, "r+")

    a = 1
    r = 1
    for i = 1:num_block_rows
        c = 1
        szi = size(A[a], 1)
        for j = 1:rows[i]
            Aj = A[a+j-1]
            szj = size(Aj, 2)
            if size(Aj, 1) != szi
                throw(ArgumentError("mismatched height in block row $(i) (expected $szi, got $(size(Aj, 1)))"))
            end
            if c-1+szj > num_cols
                throw(ArgumentError("block row $(i) has mismatched number of columns (expected $num_cols, got $(c-1+szj)"))
            end
            des_arr[r:r-1+szi, c:c-1+szj] .= @view Aj[:,:]
            c += szj
        end
        if c != num_cols+1
            throw(ArgumentError("block row $(i) has mismatched number of columns (expected $num_cols, got $(c-1))"))
        end
        r += szi
        a += rows[i]
    end
    des_arr
end


text:
using Pkg
packages = ["Documenter","Formatting", "Images", "TestImages",
"RDatasets", "InfoZIP", "ImageView", "ImageMagick","Mustache",
"StringEncodings", "TextAnalysis", "Latexify", "IJulia"]
for p in packages
    try
        if Pkg.installed(p) == nothing
            Pkg.add(p)
	end
    catch e
        Pkg.add(p)
    end
end
using Documenter
using SpectralClustering
try
   Pkg.installed("BibTeX")
catch
    Pkg.clone("https://github.com/JuliaTeX/BibTeX.jl.git")
end
try
    Pkg.installed("BibTeXFormat")
catch
    Pkg.clone("https://github.com/lucianolorenti/BibTeXFormat.jl.git")
end


text:
using Interpolations
using Test

@testset "IO" begin
    SPACE = " "

    @testset "BSpline" begin
        A = rand(8,20)

        itp = interpolate(A, BSpline(Constant()))
        @test summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant())) with element type Float64"

        itp = interpolate(A, BSpline(Constant()))
        @test summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant())) with element type Float64"

        itp = interpolate(A, BSpline(Linear()))
        @test summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Linear())) with element type Float64"

        itp = interpolate(A, BSpline(Quadratic(Reflect(OnCell()))))
        @test summary(itp) == "8×20 interpolate(OffsetArray(::Array{Float64,2}, 0:9, 0:21), BSpline(Quadratic(Reflect(OnCell())))) with element type Float64"

        itp = interpolate(A, (BSpline(Linear()), NoInterp()))
        @test summary(itp) == "8×20 interpolate(::Array{Float64,2}, (BSpline(Linear()), NoInterp())) with element type Float64"

        itp = interpolate!(copy(A), BSpline(Quadratic(InPlace(OnCell()))))
        @test summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Quadratic(InPlace(OnCell())))) with element type Float64"
    end

    @testset "Gridded" begin
        A = rand(20)
        A_x = collect(1.0:2.0:40.0)
        knots = (A_x,)
        itp = interpolate(knots, A, Gridded(Linear()))
        @test summary(itp) == "20-element interpolate((::Array{Float64,1},), ::Array{Float64,1}, Gridded(Linear())) with element type Float64"

        A = rand(8,20)
        knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
        itp = interpolate(knots, A, Gridded(Linear()))
        @test summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, Gridded(Linear())) with element type Float64"

        itp = interpolate(knots, A, (Gridded(Linear()),Gridded(Constant())))
        @test summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, (Gridded(Linear()), Gridded(Constant()))) with element type Float64"
    end

    @testset "scaled" begin
        itp = interpolate(1:1.0:10, BSpline(Linear()))
        sitp = scale(itp, -3:.5:1.5)
        @test summary(sitp) == "10-element scale(interpolate(::Array{Float64,1}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64"

        gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
        testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
        xs = -5:.5:5
        ys = -4:.2:4
        zs = Float64[testfunction(x,y) for x in xs, y in ys]
        itp2 = interpolate(zs, BSpline(Quadratic(Flat(OnGrid()))))
        sitp2 = scale(itp2, xs, ys)
        @test summary(sitp2) == "21×41 scale(interpolate(OffsetArray(::Array{Float64,2}, 0:22, 0:42), BSpline(Quadratic(Flat(OnGrid())))), (-5.0:0.5:5.0,$SPACE-4.0:0.2:4.0)) with element type Float64"
    end

    @testset "Monotonic" begin
        A = rand(20)
        A_x = collect(1.0:2.0:40.0)
        itp = interpolate(A_x, A, FritschButlandMonotonicInterpolation())
        @test summary(itp) == "20-element interpolate(::Array{Float64,1}, ::Array{Float64,1}, FritschButlandMonotonicInterpolation()) with element type Float64"
    end

    @testset "Extrapolation" begin
        A = rand(8,20)

        itpg = interpolate(A, BSpline(Linear()))
        etpg = extrapolate(itpg, Flat())
        @test summary(etpg) == "8×20 extrapolate(interpolate(::Array{Float64,2}, BSpline(Linear())), Flat()) with element type Float64"

        etpf = extrapolate(itpg, NaN)
        @test summary(etpf) == "8×20 extrapolate(interpolate(::Array{Float64,2}, BSpline(Linear())), NaN) with element type Float64"
    end
end # testset


text: print(2 x + 1)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using QuadGK
using StableRNGs
using Transits
using Test

const PLOT = get(ENV, "TEST_PLOTS", "false") == "true"
PLOT && include("plots.jl")

rng = StableRNG(2752)

@testset "Transits" begin
    include("elliptic.jl")
    include("Mn_integral.jl")
    include("poly.jl")
end


text:
using Distributions
using Random
using StaticArrays
using Test
using DocStringExtensions
using Printf

using CLIMA
using CLIMA.Atmos
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.Mesh.Filters
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using CLIMA.VariableTemplates

import CLIMA.DGmethods: boundary_state!
import CLIMA.Atmos: atmos_boundary_state!
import CLIMA.DGmethods.NumericalFluxes: boundary_flux_diffusive!

# ------------------- Description ---------------------------------------- #
# 1) Dry Rayleigh Benard Convection (re-entrant channel configuration)
# 2) Boundaries - `Sides` : Periodic (Default `bctuple` used to identify bot,top walls)
#                 `Top`   : Prescribed temperature, no-slip
#                 `Bottom`: Prescribed temperature, no-slip
# 3) Domain - 250m[horizontal] x 250m[horizontal] x 500m[vertical]
# 4) Timeend - 1000s
# 5) Mesh Aspect Ratio (Effective resolution) 1:1
# 6) Random seed in initial condition (Requires `forcecpu=true` argument)
# 7) Overrides defaults for
#               `C_smag`
#               `Courant_number`
#               `forcecpu`
#               `ref_state`
#               `solver_type`
#               `bc`
#               `sources`
# 8) Default settings can be found in src/Driver/Configurations.jl

# ------------------- Begin Boundary Conditions -------------------------- #
"""
  FixedTempNoSlip <: BoundaryCondition

Fixed temperature prescription at top and bottom walls
No slip velocity boundary conditions.
Y ≡ state
Σ ≡ diff
A ≡ aux
⁺ and ⁻ refer to exterior / interior faces

# Fields
$(DocStringExtensions.FIELDS)
"""
struct FixedTempNoSlip{FT} <: BoundaryCondition
  "Prescribed bottom wall temperature `[K]`"
  T_bot::FT
  "Prescribed top wall temperature `[K]`"
  T_top::FT
end
# Rayleigh-Benard problem with two fixed walls (prescribed temperatures)
function atmos_boundary_state!(::Union{NumericalFluxNonDiffusive, NumericalFluxGradient},
                               bc::FixedTempNoSlip,
                               m::AtmosModel,
                               Y⁺::Vars, A⁺::Vars,
                               n⁻,
                               Y⁻::Vars, A⁻::Vars,
                               bctype, t,_...)
  # Dry Rayleigh Benard Convection
  FT = eltype(Y⁺)
  @inbounds begin
    Y⁺.ρu = -Y⁻.ρu
    if bctype == 1
      E_int⁺ = Y⁺.ρ * cv_d * (bc.T_bot - T_0)
    else
      E_int⁺ = Y⁺.ρ * cv_d * (bc.T_top - T_0)
    end
    E_bc = (E_int⁺ + Y⁺.ρ * A⁺.coord[3] * grav)
    Y⁺.ρe = E_bc
  end
end
function atmos_boundary_flux_diffusive!(::CentralNumericalFluxDiffusive,
                                        bc::FixedTempNoSlip,
                                        m::AtmosModel, F,
                                        Y⁺::Vars, Σ⁺::Vars, A⁺::Vars,
                                        n⁻,
                                        Y⁻::Vars, Σ⁻::Vars, A⁻::Vars,
                                        bctype, t, _...)
  Σ⁺.∇h_tot = -Σ⁻.∇h_tot
end

boundary_state!(nf, m::AtmosModel, x...) =
  atmos_boundary_state!(nf, m.boundarycondition, m, x...)
boundary_flux_diffusive!(nf::NumericalFluxDiffusive,
                         atmos::AtmosModel,
                         F,
                         state⁺, diff⁺, aux⁺, n⁻,
                         state⁻, diff⁻, aux⁻,
                         bctype, t,
                         state1⁻, diff1⁻, aux1⁻) =
  atmos_boundary_flux_diffusive!(nf, atmos.boundarycondition, atmos,
                                 F,
                                 state⁺, diff⁺, aux⁺, n⁻,
                                 state⁻, diff⁻, aux⁻,
                                 bctype, t,
                                 state1⁻, diff1⁻, aux1⁻)
# ------------------- End Boundary Conditions -------------------------- #

const randomseed         = MersenneTwister(1)
const (xmin, ymin, zmin) = (0,0,0)
const (xmax, ymax, zmax) = (250,250,500)
const T_bot              = 299
const T_lapse            = grav/cp_d
const T_top              = T_bot - T_lapse*zmax

function init_problem!(bl, state, aux, (x,y,z), t)
  FT            = eltype(state)
  R_gas::FT     = R_d
  c_p::FT       = cp_d
  c_v::FT       = cv_d
  γ::FT         = c_p / c_v
  p0::FT        = MSLP
  δT            = sinpi(6*z/(zmax-zmin)) * cospi(6*z/(zmax-zmin)) + rand(randomseed)
  δw            = sinpi(6*z/(zmax-zmin)) * cospi(6*z/(zmax-zmin)) + rand(randomseed)
  ΔT            = grav/cp_d * z + δT
  T             = T_bot - ΔT
  P             = p0*(T/T_bot)^(grav/R_gas/T_lapse)
  ρ             = P / (R_gas * T)
  ρu, ρv, ρw    = FT(0) , FT(0) , ρ * δw
  E_int         = ρ * c_v * (T-T_0)
  E_pot         = ρ * grav * z
  E_kin         = ρ * FT(1/2) * δw^2
  ρe_tot        = E_int + E_pot + E_kin
  state.ρ       = ρ
  state.ρu      = SVector(ρu, ρv, ρw)
  state.ρe      = ρe_tot
  state.moisture.ρq_tot = FT(0)
end

function config_problem(FT, N, resolution, xmax, ymax, zmax)

    # Boundary conditions
    bc = FixedTempNoSlip{FT}(T_bot, T_top)

    # Turbulence
    C_smag = FT(0.23)

    # Set up the model
    model = AtmosModel{FT}(AtmosLESConfiguration;
                           turbulence=SmagorinskyLilly{FT}(C_smag),
                               source=(Gravity(),),
                    boundarycondition=bc,
                           init_state=init_problem!)
    ode_solver = CLIMA.ExplicitSolverType(solver_method=LSRK144NiegemannDiehlBusch)
    config = CLIMA.Atmos_LES_Configuration("DryRayleighBenardConvection",
                                           N, resolution, xmax, ymax, zmax,
                                           init_problem!,
                                           solver_type=ode_solver,
                                           model=model)
    return config
end

function main()
    CLIMA.init()
    FT = Float64
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(10)
    # Time integrator setup
    t0 = FT(0)
    CFLmax = FT(0.90)
    timeend = FT(1000)

    @testset "DryRayleighBenardTest" begin
      for Δh in Δh
        Δv = Δh
        resolution = (Δh, Δh, Δv)
        driver_config = config_problem(FT, N, resolution, xmax, ymax, zmax)
        solver_config = CLIMA.setup_solver(t0, timeend, driver_config,
                                           forcecpu=true, Courant_number=CFLmax)
        # User defined callbacks (TMAR positivity preserving filter)
        cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do (init=false)
            Filters.apply!(solver_config.Q, 6, solver_config.dg.grid, TMARFilter())
            nothing
        end
        result = CLIMA.invoke!(solver_config;
                              user_callbacks=(cbtmarfilter,),
                              check_euclidean_distance=true)
        # result == engf/eng0
        @test isapprox(result,FT(1); atol=1.5e-2)
      end
    end
end

main()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# Tracing charged particle in the structured SWMF outputs from MHD-EPIC simulation.
# For more details about the format, please checkout Batsrus.jl.
#
# Hongyang Zhou, hyzhou@umich.edu

using Statistics: mean
using Batsrus
using TestParticle
using OrdinaryDiffEq
using Meshes
using FieldTracer
using PyPlot

## Utility functions

function initial_conditions(i)
   j = i - 1
   [x[xc_], y[yc_+j], z[zc_], Uix[xc_,yc_+j,zc_], Uiy[xc_,yc_+j,zc_], Uiz[xc_,yc_+j,zc_]]
end

"Initial state perturbation for EnsembleProblem."
function prob_func(prob, i, repeat)
   remake(prob,u0=initial_conditions(i))
end

## Data processing

filename = "3d_var_region0_0_t00002500_n00106468.out"
data = readdata(filename)

var = getvars(data, ["Bx", "By", "Bz", "Ex", "Ey", "Ez", "uxs0", "uys0", "uzs0", "uxs1", "uys1", "uzs1"])

const RG = 2634e3 # [m]

x = range(extrema(data.x[:,:,:,1])..., length=size(data.x, 1)) .* RG
y = range(extrema(data.x[:,:,:,2])..., length=size(data.x, 2)) .* RG
z = range(extrema(data.x[:,:,:,3])..., length=size(data.x, 3)) .* RG
B = zeros(Float32, 3, length(x), length(y), length(z)) # [T]
E = zeros(Float32, 3, length(x), length(y), length(z)) # [V/m]

# Convert into SI units
B[1,:,:,:] .= var["Bx"] .* 1e-9
B[2,:,:,:] .= var["By"] .* 1e-9
B[3,:,:,:] .= var["Bz"] .* 1e-9
E[1,:,:,:] .= var["Ex"] .* 1e-6
E[2,:,:,:] .= var["Ey"] .* 1e-6
E[3,:,:,:] .= var["Ez"] .* 1e-6

Δx = x[2] - x[1]
Δy = y[2] - y[1]
Δz = z[2] - z[1]

grid = Meshes.CartesianGrid((length(x)-1, length(y)-1, length(z)-1),
   (x[1], y[1], z[1]),
   (Δx, Δy, Δz))

## Initial conditions

Uex, Uey, Uez = var["uxs0"] .* 1e3, var["uys0"] .* 1e3, var["uzs0"] .* 1e3
Uix, Uiy, Uiz = var["uxs1"] .* 1e3, var["uys1"] .* 1e3, var["uzs1"] .* 1e3

xc_ = floor(Int, length(x)/2) + 1
yc_ = floor(Int, length(y)/2) + 1
zc_ = floor(Int, length(z)/2) + 1

stateinit_e = [x[xc_], y[yc_], z[zc_], Uex[xc_,yc_,zc_], Uey[xc_,yc_,zc_], Uez[xc_,yc_,zc_]]
stateinit_p = [x[xc_], y[yc_], z[zc_], Uix[xc_,yc_,zc_], Uiy[xc_,yc_,zc_], Uiz[xc_,yc_,zc_]]

param_electron = prepare(grid, E, B, species=Electron)
tspan_electron = (0.0, 0.1)

param_proton = prepare(grid, E, B, species=Proton)
tspan_proton = (0.0, 10.0)
trajectories = 5

#=
prob_e = ODEProblem(trace_numeric!, stateinit_e, tspan_electron, param_electron)
ensemble_prob = EnsembleProblem(prob_e, prob_func=prob_func)
sol_e = solve(prob_e, Tsit5(); save_idxs=[1,2,3], alg_hints=[:nonstiff])
sol_e = solve(ensemble_prob, Tsit5(), EnsembleThreads();
   trajectories, save_idxs=[1,2,3])
=#

prob_p = ODEProblem(trace_numeric!, stateinit_p, tspan_proton, param_proton)
ensemble_prob = EnsembleProblem(prob_p, prob_func=prob_func)
#sol_p = solve(prob_p; save_idxs=[1,2,3], alg_hints=[:nonstiff])
sol_p = solve(ensemble_prob, AutoTsit5(Rosenbrock23()), EnsembleThreads();
   trajectories, save_idxs=[1,2,3])

## Visualization

using3D()
fig = plt.figure(figsize=(10,6))
ax = fig.gca(projection="3d")

## Field tracing

for i in 1:10:length(x)
   xs, ys, zs = x[i], 0.0, 0.0
   x1, y1, z1 = trace(B[1,:,:,:], B[2,:,:,:], B[3,:,:,:], xs, ys, zs, x, y, z, ds=0.2, maxstep=1000)
   line = ax.plot(x1 ./ RG, y1 ./RG, z1 ./ RG, "k-", alpha=0.3)
end

for i in 1:10:length(y)
   xs, ys, zs = x[xc_], y[i], z[zc_]
   x1, y1, z1 = trace(B[1,:,:,:], B[2,:,:,:], B[3,:,:,:], xs, ys, zs, x, y, z, ds=0.2, maxstep=1000)
   line = ax.plot(x1 ./ RG, y1 ./RG, z1 ./ RG, "k-", alpha=0.3)
end

#xyz = plot(sol_e, vars=(1,2,3), xlabel="x", ylabel="y", zlabel="z",
#   aspect_ratio=:equal, label="electron", lw=2)

#ax.plot(sol_e[1,:], sol_e[2,:], sol_e[3,:], label="electron")

n = 200 # number of timepoints
#ts = range(0, stop=tspan_electron[2], length=n)
#ax.plot(sol_e(ts,idxs=1) ./ RG, sol_e(ts,idxs=2) ./ RG, sol_e(ts,idxs=3) ./ RG, label="electron")

ts = range(0, stop=tspan_proton[2], length=n)
for i = 1:trajectories
   if sol_p[i].t[end] < tspan_proton[2]
      ts⁺ = range(0, stop=sol_p[i].t[end], length=n)
      ax.plot(sol_p[i](ts⁺,idxs=1) ./ RG, sol_p[i](ts⁺,idxs=2) ./ RG, sol_p[i](ts⁺,idxs=3) ./ RG, label="proton $i", lw=1.5)
   else
      ax.plot(sol_p[i](ts,idxs=1) ./ RG, sol_p[i](ts,idxs=2) ./ RG, sol_p[i](ts,idxs=3) ./ RG, label="proton $i", lw=1.5)
   end
end

#ax.plot(sol_p[1,:], sol_p[2,:], sol_p[3,:], label="proton")

ax.legend()
title("Particle trajectory near Ganymede's magnetopause from PIC")
xlabel("x [Rg]")
ylabel("y [Rg]")
zlabel("z [Rg]")

ax.set_box_aspect([1.17,4,4])
TestParticle.set_axes_equal(ax)

#plotdata(data, "bz", plotmode="contbar", cut="y", cutPlaneIndex=40, level=20)


text:
# test complexify.jl

"""
  This function tests whether compexify.jl functions work correctly
"""
function test_complexify()

  @testset "----- Testing Complexify.jl ------" begin

    a = 1.0
    b = -1.0
    @test isapprox( absvalue(a), abs(a)) atol=1e-15
    @test isapprox( absvalue(b), abs(b)) atol=1e-15
    @test isapprox( absvalue_rev(a, 2)[1], absvalue(a)) atol=1e-15
    @test isapprox( absvalue_rev(a, 2)[2], 2) atol=1e-15
    @test isapprox( absvalue_rev(b, 2)[1], absvalue(b)) atol=1e-15
    @test isapprox( absvalue_rev(b, 2)[2], -2) atol=1e-15

    c = complex(1.0, 1.0)
    d = complex(1.0, -1.0)
    f = complex(-1.0, 1.0)
    g = complex(-1.0, -1.0)
    j = complex(0.0, 0.0)

    @test isapprox( absvalue(c), complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue_rev(c, 2)[1], complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue_rev(c, 2)[2], 2) atol=1e-15

    @test isapprox( absvalue(d), complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue_rev(d, 2)[1], complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue_rev(d, 2)[2], 2) atol=1e-15

    @test isapprox( absvalue(f), complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue_rev(f, 2)[1], complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue_rev(f, 2)[2], -2) atol=1e-15

    @test isapprox( absvalue(g), complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue_rev(g, 2)[1], complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue_rev(g, 2)[2], -2) atol=1e-15

    # absvalue2
    @test isapprox( absvalue2(c), complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue2(d), complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue2(f), complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue2(g), complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue2(j), complex(0.0, 0.0)) atol=1e-15

    # absvalue3
    @test isapprox( absvalue3(c), complex(1.0, 1.0)) atol=1e-15
    @test isapprox( absvalue3(d), complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue3(f), complex(1.0, -1.0)) atol=1e-15
    @test isapprox( absvalue3(g), complex(1.0, 1.0)) atol=1e-15
    #@test isapprox( absvalue3(j), complex(0.0, 0.0)) atol=1e-15

    # the spline is always below the y=x line
    @test absvalue2(1e-14) < 1e-14
    @test absvalue2(1e-14) > 0
    @test absvalue2(-1e-14) < 1e-14
    @test absvalue2(-1e-14) > 0

    # test derivative
    h = 1e-20; pert = Complex128(0, h)
    vals = [2.0, 2e-13,  1e-14, -2e-13, -1e-14, 0.0, 0 + eps(), 0 - eps()]
    for val in vals
      @test isapprox(imag(absvalue2(val + pert))/h, absvalue2_deriv(val)) atol=2e-14
      # this function is discontinuous at 0 when delta <= 0
      if val != 0
        @test isapprox(imag(absvalue3(val + pert))/h, absvalue3_deriv(val)) atol=2e-14
      end
    end

    x = 1.0; y = 1.0
    checkdiff(x, y)
    checkderiv(x, y)
    check_atan2_rev(x, y)

    x = 0.1; y = 1.0
    checkdiff(x, y)
    checkderiv(x, y)
    check_atan2_rev(x, y)

    x = 0.01; y = 1.0
    checkdiff(x, y)
    checkderiv(x, y)
    check_atan2_rev(x, y)

    x = 0.1; y = -1.0
    checkdiff(x, y)
    checkderiv(x, y)
    check_atan2_rev(x, y)

    x = -0.001; y = -1.0
    checkdiff(x, y)
    checkderiv(x, y)
    check_atan2_rev(x, y)

  end
end

"""
  Helper function for testing atan2 value
"""
function checkdiff(x::Real, y::Real)

  x2 = Complex(x, 0)
  y2 = Complex(y, 0)

  theta1 = atan2(y, x)
  theta2 = atan2(y2, x2)

  @test isapprox(theta1, theta2) atol=2*eps()

  return theta1 - theta2
end

"""
  Helper function for testing atan2 derivative (finite difference)
"""
function checkderiv(x::Real, y::Real)

  h = 1e-20
  hfd = 1e-6
  x2 = Complex128(x, h)

  d1 = imag(atan2(y, x2))/h
  d1fd = (atan2(y, x + hfd) - atan2(y, x))/hfd

  @test isapprox(d1, d1fd) atol=1e-6

  y2 = Complex128(y, h)
  d1 = imag(atan2(y2, x))/h
  d1fd = (atan2(y + hfd, x) - atan2(y, x))/hfd

  @test isapprox(d1, d1fd) atol=1e-6

  return nothing
end

"""
  Helper function for checking atan2_rev
"""
function check_atan2_rev(x::Real, y::Real)

  h = 1e-20
  pert = Complex128(0, h)

  x2 = Complex128(x, 0)
  y2 = Complex128(y, 0)

  theta = atan2(y, x)

  y2 += pert
  d1 = imag(atan2(y2, x2))/h
  y2 -= pert

  x2 += pert
  d2 = imag(atan2(y2, x2))/h
  x2 -= pert

  theta_bar = 1
  theta2, y_bar, x_bar = atan2_rev(y, x, theta_bar)

  @test isapprox(d1, y_bar) atol=1e-14
  @test isapprox(d2, x_bar) atol=1e-14
  @test isapprox(theta, theta2) atol=1e-14

  return nothing
end

#test_complexify()
add_func1!(EulerTests, test_complexify, [TAG_COMPLEX, TAG_SHORTTEST])