text:
"""
Implementation of Evolution Strategy: (μ/ρ(+/,)λ)-ES

The constructor takes following keyword arguments:

- `initStrategy`: an initial strategy description, (default: empty)
- `recombination`: ES recombination function for population (default: `first`), see [Crossover](@ref)
- `srecombination`: ES recombination function for strategies (default: `first`), see [Crossover](@ref)
- `mutation`: [Mutation](@ref) function for population (default: `first`)
- `smutation`: [Mutation](@ref) function for strategies (default: `identity`)
- `μ`/`mu`: the number of parents
- `ρ`/`rho`: the mixing number, ρ ≤ μ, (i.e., the number of parents involved in the procreation of an offspring)
- `λ`/`lambda`: the number of offspring
- `selection`: the selection strategy `:plus` or `:comma` (default: `:plus`)
"""
struct ES <: AbstractOptimizer
    initStrategy::AbstractStrategy
    recombination::Function
    srecombination::Function
    mutation::Function
    smutation::Function
    μ::Integer
    ρ::Integer
    λ::Integer
    selection::Symbol

    ES(; initStrategy::AbstractStrategy = NoStrategy(),
        recombination::Function = first,
        srecombination::Function = first,
        mutation::Function = ((r,s) -> r),
        smutation::Function = identity,
        μ::Integer = 1,
        mu::Integer = μ,
        ρ::Integer = μ,
        rho::Integer = ρ,
        λ::Integer = 1,
        lambda::Integer = λ,
        selection::Symbol = :plus) =
        new(initStrategy, recombination, srecombination, mutation, smutation,
            mu, rho, lambda, selection)
end
population_size(method::ES) = method.μ
default_options(method::ES) = (iterations=1000, abstol=1e-10)
summary(m::ES) = "($(m.μ)/$(m.ρ)$(m.selection == :plus ? '+' : ',')$(m.λ))-ES"
show(io::IO,m::ES) = print(io, summary(m))

mutable struct ESState{T,IT,ST} <: AbstractOptimizerState
    N::Int
    fitness::Vector{T}
    strategies::Vector{ST}
    fittest::IT
end
value(s::ESState) = first(s.fitness)
minimizer(s::ESState) = s.fittest
strategy(s::ESState) = first(s.strategies)

"""Initialization of ES algorithm state"""
function initial_state(method::ES, options, objfun, population)
    T = typeof(value(objfun))
    individual = first(population)
    N = length(individual)

    # Populate fitness and strategies
    fitness = map(i -> value(objfun, i), population)
    strategies = Array{AbstractStrategy}(undef, method.μ)
    for i in 1:method.μ
        strategies[i] = copy(method.initStrategy)
    end

    # setup initial state
    return ESState(N, fitness, strategies, copy(individual))
end

function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::ES, itr) where {IT}
    @unpack initStrategy,recombination,srecombination,mutation,smutation,μ,ρ,λ,selection = method

    @assert ρ <= μ "Number of parents involved in the procreation of an offspring should be no more then total number of parents"
    if selection == :comma
        @assert μ < λ "Offspring population must be larger then parent population"
    end
    offspring = Array{IT}(undef, λ)
    fitoff = fill(Inf, λ)
    stgoff = Array{AbstractStrategy}(undef, λ)

    for i in 1:λ
        # Recombine the ρ selected parents to form a recombinant individual
        if ρ == 1
            j = rand(1:μ)
            recombinantStrategy = state.strategies[j]
            recombinant = copy(population[j])
        else
            idx = randperm(μ)[1:ρ]
            recombinantStrategy = srecombination(state.strategies[idx])
            recombinant = recombination(population[idx])
        end

        # Mutate the strategy parameter set of the recombinant
        stgoff[i] = smutation(recombinantStrategy)

        # Mutate the objective parameter set of the recombinant using the mutated strategy parameter set
        # to control the statistical properties of the object parameter mutation
        off = mutation(recombinant, stgoff[i])

        # Apply constraints
        offspring[i] = apply!(constraints, off)

        # Evaluate fitness
        fitoff[i] = value(objfun, offspring[i])
    end
    # apply penalty to fitness
    penalty!(fitoff, constraints, offspring)

    # Select new parent population
    if selection == :plus
        idxs = sortperm(vcat(state.fitness, fitoff))[1:μ]
        skip = idxs[idxs.<=μ]
        for i = 1:μ
            if idxs[i] ∉ skip
                ii = idxs[i] - μ
                population[i] = offspring[ii]
                state.strategies[i] = stgoff[ii]
                state.fitness[i] = fitoff[ii]
            end
        end
    else
        idxs = sortperm(fitoff)[1:μ]
        for (i,j) in enumerate(idxs)
            population[i] = offspring[j]
            state.strategies[i] = stgoff[j]
            state.fitness[i] = fitoff[j]
        end
    end

    # indicate a fittest individual
    state.fittest = first(population)

    return false
end


text:
module PhononAction

using ..Models: HolsteinModel, SSHModel
using ..Utilities: get_index

export calc_Sb, calc_dSbdxdx!

"""
Calculates the pure phonon action Sb such that exp{-Sb}.
"""
function calc_Sb(holstein::HolsteinModel{T1,T2}, shifted::Bool=false) where {T1,T2}

    x   = holstein.x::Vector{T1}
    N   = holstein.Nsites::Int
    Lτ  = holstein.Lτ::Int
    Δτ  = holstein.Δτ::T1
    ω   = holstein.ω::Vector{T1}
    ω₄  = holstein.ω₄::Vector{T1}
    λ   = holstein.λ::Vector{T1}
    Sb = 0.0::T1

    # iterate over sites in lattice
    for i in 1:N
        # iterate over time slice in lattice
        for τ in 1:Lτ
            # get τ-1 accounting for periodic boundary conditions
            τm1 = mod1(τ-1,Lτ)
            # xᵢ(τ)
            xᵢτ = x[get_index(τ,i,Lτ)]
            # xᵢ(τ-1)
            xᵢτm1 = x[get_index(τm1,i,Lτ)]
            # calculate potential energy
            Sb += ω[i]^2*xᵢτ^2/2 + ω₄[i]*xᵢτ^4 - λ[i]*xᵢτ*shifted
            # calculate kintetic energy
            Sb += (xᵢτ-xᵢτm1)^2/Δτ^2/2
        end
    end

    # calculate phonon potential energy associated with phonon dispersion
    if length(holstein.ωᵢⱼ)>0
        ωᵢⱼ                = holstein.ωᵢⱼ::Vector{T2}
        sign_ωᵢⱼ           = holstein.sign_ωᵢⱼ::Vector{Int}
        neighbor_table_ωᵢⱼ = holstein.neighbor_table_ωᵢⱼ::Matrix{Int}
        # iterate over dispersive phonon modes
        for m in 1:length(ωᵢⱼ)
            ωᵢⱼ² = ωᵢⱼ[m] * ωᵢⱼ[m]
            sgn = sign_ωᵢⱼ[m]
            # getting pair of neighboring sites
            i = neighbor_table_ωᵢⱼ[1,m]
            j = neighbor_table_ωᵢⱼ[2,m]
            # iterating over time slices
            for τ in 1:L
                # indexing offset into vectors associated with τ time slice
                indx_i = get_index(τ,i,Lτ)
                indx_j = get_index(τ,j,Lτ)
                # updating partial derivative
                Sb += ωᵢⱼ² * ( x[indx_i] + sgn*x[indx_j] )^2/2
            end
        end
    end

    # necessary scaling factor as defintiion of Sb includes a Δτ out front
    Sb = Δτ * Sb

    return Sb
end

function calc_Sb(ssh::SSHModel{T1,T2,T3}, shifted::Bool=false) where {T1,T2,T3}

    x  = ssh.x::Vector{T1}
    N  = ssh.Nph::Int
    Lτ = ssh.Lτ::Int
    Δτ = ssh.Δτ::T1
    ω  = ssh.ω::Vector{T1}
    ω₄ = ssh.ω₄::Vector{T1}
    Sb = 0.0

    # iterate over phonons in lattice
    for i in 1:N
        # τ=0 field associated with phonon
        field_τ0 = get_index(1,i,Lτ)
        # if primary phonon
        if ssh.primary_field[field_τ0] == field_τ0
            # action associated with current phonon
            sb = 0.0
            # iterate over time slice in lattice
            for τ in 1:Lτ
                # get τ-1 accounting for periodic boundary conditions
                τm1   = mod1(τ-1,Lτ)
                # get field
                field = get_index(τ,i,Lτ)
                # xᵢ(τ)
                xᵢτ   = x[field]
                # xᵢ(τ-1)
                xᵢτm1 = x[get_index(τm1,i,Lτ)]
                # calculate potential energy
                sb   += Δτ*ω[i]^2*xᵢτ^2/2 + Δτ*ω₄[i]*xᵢτ^4
                # calculate kintetic energy
                sb   += (xᵢτ-xᵢτm1)^2/Δτ/2
            end
            # add to Sb total, normalizing by the number of equivalent fields there are
            Sb += sb
        end
    end

    return Sb
end

"""
Calculates the dervative phonon action with respect to each phonon field and adds that value in place
to the vector dSbdx.
"""
function calc_dSbdx!(dSbdx::Vector{T2}, holstein::HolsteinModel{T1,T2}, shifted::Bool=false)  where {T1<:AbstractFloat,T2<:Number}

    @assert length(dSbdx)==holstein.Ndof

    x          = holstein.x::Vector{T1}
    nsites     = holstein.Nsites::Int
    Lτ         = holstein.Lτ::Int
    Δτ         = holstein.Δτ::T1
    ω          = holstein.ω::Vector{T1}
    ω₄         = holstein.ω₄::Vector{T1}
    λ          = holstein.λ::Vector{T1}

    #####################################################
    ## Calculating Derivative Phonon Action Associated ##
    ## With Local Phonon Frequency And Phonon Momentum ##
    #####################################################]

    # iterating over site
    @fastmath @inbounds for site in 1:nsites
        Δτω²  = Δτ * ω[site] * ω[site]
        Δτ4ω₄ = Δτ * 4 * ω₄[site]
        Δτλ   = Δτ * λ[site]
        # iterating over time slices
        for τ in 1:Lτ
            # get τ+1 accounting for periodic boundary conditions
            τp1 = τ%Lτ+1
            # get τ-1 accounting for periodic boundary conditions
            τm1 = (τ-2+Lτ)%Lτ+1
            # indexing offset into vectors associated with τ time slice
            indx_τ = get_index(τ,site,Lτ)
            # indexing offset into vectors associated with τ+1 time slice
            indx_τp1 = get_index(τp1,site,Lτ)
            # indexing offset into vectors associated with τ-1 time slice
            indx_τm1 = get_index(τm1,site,Lτ)
            # phonon field at current time slice
            xτ = x[indx_τ]
            # updating partial derivative
            dSbdx[indx_τ] += Δτω² * xτ - Δτλ * shifted # derivative of Δτ⋅ω²/2⋅x² term
            dSbdx[indx_τ] += Δτ4ω₄ * xτ * xτ * xτ      # derivative of Δτ⋅ω₄⋅x⁴ term
            dSbdx[indx_τ] -= ( x[indx_τp1] + x[indx_τm1] - 2.0*xτ )/Δτ
        end
    end

    #############################################################
    ## Calculating Derivative Of Phonon Action Associated With ##
    ##               Dispersive Phonon Modes                   ##
    #############################################################

    if length(holstein.ωᵢⱼ)>0
        ωᵢⱼ                = holstein.ωᵢⱼ::Vector{T2}
        sign_ωᵢⱼ           = holstein.sign_ωᵢⱼ::Vector{Int}
        neighbor_table_ωᵢⱼ = holstein.neighbor_table_ωᵢⱼ::Matrix{Int}
        # iterate over dispersive phonon modes
        for m in 1:length(ωᵢⱼ)
            Δτωᵢⱼ² = Δτ * ωᵢⱼ[m] * ωᵢⱼ[m]
            sgn = sign_ωᵢⱼ[m]
            # getting pair of neighboring sites
            i = neighbor_table_ωᵢⱼ[1,m]
            j = neighbor_table_ωᵢⱼ[2,m]
            # iterating over time slices
            for τ in 1:Lτ
                # indexing offset into vectors associated with τ time slice
                indx_i = get_index(τ,i,Lτ)
                indx_j = get_index(τ,j,Lτ)
                # updating partial derivative
                Δ = Δτωᵢⱼ² * ( x[indx_i] + sgn*x[indx_j] )
                dSbdx[indx_i] += Δ
                dSbdx[indx_j] += sgn*Δ
            end
        end
    end

    return nothing
end

function calc_dSbdx!(dSbdx::Vector{T2}, ssh::SSHModel{T1,T2,T3}, shifted::Bool=false)  where {T1,T2,T3}

    @assert length(dSbdx)==ssh.Ndof

    x   = ssh.x::Vector{T1}
    Nph = ssh.Nph::Int
    Lτ  = ssh.Lτ::Int
    Δτ  = ssh.Δτ::T1
    ω   = ssh.ω::Vector{T1}
    ω₄  = ssh.ω₄::Vector{T1}

    #####################################################
    ## Calculating Derivative Phonon Action Associated ##
    ## With Local Phonon Frequency And Phonon Momentum ##
    #####################################################

    # iterating over phonons in lattice
    @fastmath @inbounds for i in 1:Nph
        Δτω²  = Δτ * ω[i] * ω[i]
        Δτ4ω₄ = Δτ * 4 * ω₄[i]
        # iterating over time slices
        for τ in 1:Lτ
            # get τ+1 accounting for periodic boundary conditions
            τp1       = mod1(τ+1,Lτ)
            # get τ-1 accounting for periodic boundary conditions
            τm1       = mod1(τ-1,Lτ)
            # indexing offset into vectors associated with τ time slice
            field_τ   = get_index(τ,i,Lτ)
            # indexing offset into vectors associated with τ+1 time slice
            field_τp1 = get_index(τp1,i,Lτ)
            # indexing offset into vectors associated with τ-1 time slice
            field_τm1 = get_index(τm1,i,Lτ)
            # phonon field at current time slice
            xτ        = x[field_τ]
            # updating partial derivative
            dsbdx     = Δτω² * xτ # derivative of Δτ⋅ω²/2⋅x² term
            dsbdx    += Δτ4ω₄ * xτ * xτ * xτ # derivative of Δτ⋅ω₄⋅x⁴ term.
            dsbdx    -= ( x[field_τp1] + x[field_τm1] - 2.0*xτ )/Δτ # kinetic energy term
            # increment total derivative value
            dSbdx[field_τ] += dsbdx
        end
    end

    return nothing
end

end

ext: jl

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
function build_sysimg(
    m::Module,
    options::ComoniconOptions.Comonicon;
    # allow override these two options
    incremental = options.sysimg.incremental,
    filter_stdlibs = options.sysimg.filter_stdlibs,
    cpu_target = options.sysimg.cpu_target,
)
    lib = pkgdir(m, options.sysimg.path, "lib")
    if !ispath(lib)
        @info "creating library path: $lib"
        mkpath(lib)
    end

    @info "compile under project: $(pkgdir(m))"
    @info options.sysimg

    if incremental != options.sysimg.incremental
        @info "incremental override to $incremental"
    end

    if filter_stdlibs != options.sysimg.filter_stdlibs
        @info "filter_stdlibs override to $filter_stdlibs"
    end

    if cpu_target != options.sysimg.cpu_target
        @info "cpu_target override to $cpu_target"
    end

    exec_file = String[pkgdir(m, x) for x in options.sysimg.precompile.execution_file]
    stmt_file = String[pkgdir(m, x) for x in options.sysimg.precompile.statements_file]

    create_sysimage(
        nameof(m);
        sysimage_path = joinpath(lib, sysimg(options.name)),
        incremental = incremental,
        filter_stdlibs = filter_stdlibs,
        project = pkgdir(m),
        precompile_execution_file = exec_file,
        precompile_statements_file = stmt_file,
        cpu_target = cpu_target,
    )

    return
end

function build_application(m::Module, options::ComoniconOptions.Comonicon)
    build_dir = pkgdir(m, options.application.path, options.name)
    if !ispath(build_dir)
        @info "creating build path: $build_dir"
        mkpath(build_dir)
    end

    @info "application options: " options.application

    exec_file = String[pkgdir(m, x) for x in options.application.precompile.execution_file]
    stmt_file = String[pkgdir(m, x) for x in options.application.precompile.statements_file]

    create_app(
        pkgdir(m),
        build_dir;
        app_name = options.name,
        precompile_execution_file = exec_file,
        precompile_statements_file = stmt_file,
        incremental = options.application.incremental,
        filter_stdlibs = options.application.filter_stdlibs,
        force = true,
        cpu_target = options.application.cpu_target,
        c_driver_program = get_c_driver_program(options.application),
    )

    if options.install.completion
        @info "generating completion scripts"
        build_completion(m, options)
    end
    return
end

function get_c_driver_program(options::ComoniconOptions.Application)
    default_c_driver_program = pkgdir(PackageCompiler, "src", "embedding_wrapper.c")
    return if isnothing(options.c_driver_program)
        default_c_driver_program
    else
        options.c_driver_program
    end
end

function build_tarball(m::Module, options::ComoniconOptions.Comonicon)
    build_tarball_app(m, options)
    build_tarball_sysimg(m, options)
    return
end

function build_tarball_app(m::Module, options::ComoniconOptions.Comonicon)
    isnothing(options.application) && return
    @info "building application"
    build_application(m, options)
    # pack tarball
    tarball = tarball_name(m, options.name; application = true)
    @info "creating application tarball $tarball"
    cd(pkgdir(m, options.application.path)) do
        run(`tar -czvf $tarball $(options.name)`)
    end
    return
end

function build_tarball_sysimg(m::Module, options::ComoniconOptions.Comonicon)
    isnothing(options.sysimg) && return

    @info "building system image"
    build_sysimg(m, options)
    # pack tarball
    tarball = tarball_name(m, options.name)
    @info "creating system image tarball $tarball"
    cd(pkgdir(m, options.sysimg.path)) do
        run(`tar -czvf $tarball lib`)
    end
    return
end

function download_sysimg(m::Module, options::ComoniconOptions.Comonicon)
    url = sysimg_url(m, options)

    try
        if !isnothing(url)
            tarball = download(url)
            path = pkgdir(m, options.sysimg.path)
            unpack(tarball, path)
            # NOTE: sysimg won't be shared, so we can just remove it
            isfile(tarball) && rm(tarball)
            return
        end
    catch e
        @warn "fail to download $url, building the system image locally"
        # force incremental build
    end
    build_sysimg(m, options; incremental = true, filter_stdlibs = false, cpu_target = "native")
    return
end

function build_completion(m::Module, options::ComoniconOptions.Comonicon)
    completion_dir = pkgdir(m, options.application.path, options.name, "completions")
    if !ispath(completion_dir)
        @info "creating path: $completion_dir"
        mkpath(completion_dir)
    end

    for sh in ["zsh"]
        script = completion_script(sh, m)
        script === nothing && continue
        write(joinpath(completion_dir, "$sh.completion"), script)
    end
    return
end

function sysimg_url(mod::Module, options::ComoniconOptions.Comonicon)
    name = options.name
    isnothing(options.download) && return
    host = options.download.host

    if host == "github.com"
        url =
            "https://github.com/" *
            options.download.user *
            "/" *
            options.download.repo *
            "/releases/download/"
    else
        error("host $host is not supported, please open an issue at $COMONICON_URL")
    end

    tarball = tarball_name(mod, name)
    url *= "v$(Comonicon.get_version(mod))/$tarball"
    return url
end

function tarball_name(m::Module, name::String; application::Bool = false)
    if application
        return "$name-application-$(Comonicon.get_version(m))-$(osname())-$(Sys.ARCH).tar.gz"
    else
        return "$name-sysimg-$(Comonicon.get_version(m))-julia-$VERSION-$(osname())-$(Sys.ARCH).tar.gz"
    end
end

"""
    osname()

Return the name of OS, will be used in building tarball.
"""
function osname()
    return Sys.isapple() ? "darwin" :
           Sys.islinux() ? "linux" :
           error("unsupported OS, please open an issue to request support at $COMONICON_URL")
end

"""
    cmd_script(mod, shadow; kwargs...)

Generates a shell script that can be use as the entry of
`mod.command_main`.

# Arguments

- `mod`: a module that contains the commands and the entry.
- `shadow`: location of a Julia script that calls the actual `mod.command_main`.

# Keywords

- `exename`: The julia executable name, default is [`default_exename`](@ref).
- `sysimg`: System image to use, default is `nothing`.
- `project`: the project path of the CLI.
- `compile`: julia compile level, can be [:yes, :no, :all, :min]
- `optimize`: julia optimization level, default is 2.
"""
function cmd_script(
    mod::Module,
    shadow::String;
    project::String = pkgdir(mod),
    exename::String = default_exename(),
    sysimg = nothing,
    compile = nothing,
    optimize = 2,
)

    head = "#!/bin/sh\n"
    if (project !== nothing) && ispath(project)
        head *= "JULIA_PROJECT=$project "
    end
    head *= exename
    script = String[head]

    if sysimg !== nothing
        push!(script, "-J$sysimg")
    end

    if compile in [:yes, :no, :all, :min]
        push!(script, "--compile=$compile")
    end

    push!(script, "-O$optimize")
    push!(script, "--startup-file=no")
    push!(script, "-- $shadow \"\$@\"")

    return join(script, " \\\n    ")
end

function completion_script(sh::String, m::Module)
    isdefined(m, :CASTED_COMMANDS) || error("cannot find Comonicon CLI entry")
    haskey(m.CASTED_COMMANDS, "main") || error("cannot find Comonicon CLI entry")
    main = m.CASTED_COMMANDS["main"]

    if sh == "zsh"
        return emit_zshcompletion(main)
    else
        @warn(
            "$sh autocompletion is not supported, " *
            "please open an issue at $COMONICON_URL for feature request."
        )
    end
    return
end

"""
    detect_shell()

Detect shell type via `SHELL` environment variable.
"""
function detect_shell()
    haskey(ENV, "SHELL") || error("cannot find available shell command")
    return basename(ENV["SHELL"])
end

relhome(path) = joinpath("\$HOME", relpath(path, homedir()))

function contains_path(rcfile::String, path::String, env = ENV)
    bin = joinpath(path, "bin")
    return contains_shell_var(rcfile, "PATH", bin, env)
end

function contains_fpath(rcfile::String, path::String, env = ENV)
    completions = joinpath(path, "completions")
    return contains_shell_var(rcfile, "FPATH", completions, env)
end

function contains_shell_var(rcfile, var, value, env)
    if !haskey(env, var)
        for line in readlines(rcfile)
            if strip(line) == "export $var=\"$(relhome(value)):\$$var\"" ||
               strip(line) == "export $var=\"$value:\$$var\""
                return true
            end
        end
        return false
    end

    for each in split(env[var], ":")
        (each == relhome(value) || each == value) && return true
    end
    return false
end

function export_shell_var(rcfile, var, path)
    isempty(rcfile) && return
    script = """
    # generated by Comonicon
    export $var="$(relhome(path)):\$$var"
    """
    open(rcfile, "a") do io
        write(io, "\n" * script * "\n")
    end
end

function write_path(rcfile::String, path::String)
    bin = joinpath(path, "bin")
    return export_shell_var(rcfile, "PATH", bin)
end

function write_fpath(rcfile::String, path::String)
    completions = joinpath(path, "completions")
    return export_shell_var(rcfile, "FPATH", completions)
end

function detect_rcfile(shell::String=detect_shell(), home::String=homedir())
    if shell == "zsh"
        rcfile = joinpath((haskey(ENV, "ZDOTDIR") ? ENV["ZDOTDIR"] : home), ".zshrc")
    elseif shell == "bash"
        rcfile = joinpath(home, ".bashrc")
        if !isfile(rcfile)
            rcfile = joinpath(home, ".bash_profile")
        end
        return rcfile
    else
        @warn "auto installation for $shell is not supported, please open an issue under Comonicon.jl"
        return
    end
end

function install_env_path(install_path::String, shell::String=detect_shell(), home::String=homedir(), env=ENV, yes::Bool = false)
    rcfile = detect_rcfile(shell, home)
    isnothing(rcfile) && return

    msg = "cannot detect $install_path/bin in PATH, do you want to add it in PATH?"
    if !contains_path(rcfile, install_path, env) && prompt(msg, yes)
        write_path(rcfile, install_path)
    end

    msg = "cannot detect $install_path/completions in FPATH, do you want to add it in FPATH?"
    if !contains_fpath(rcfile, install_path, env) && prompt(msg, yes)
        write_fpath(rcfile, install_path)
    end

    @info "open a new terminal, or source $rcfile to enable the new PATH."
    return
end


text:
struct DBFSComment
    database::SQLite.DB
    pkey::Int
    campaign::DBCampaign
    person::DBPerson
    datetime::DateTime
    comment::String
end

function Base.write(db::SQLite.DB, ::Type{DBFSComment}, campaign::Int, dbp::DBPerson, comment::String)
    stmt1 = SQLite.Stmt(db, "INSERT INTO CAMPAIGNCOMMENT ( CAMPAIGN, PERSON, DATETIME, COMMENT ) VALUES ( ?, ?, ?, ? );")
    q = DBInterface.execute(stmt1, ( campaign, dbp.pkey, Dates.datetime2julian(now()), comment ))
    return DBInterface.lastrowid(q)
end
function Base.write(db::SQLite.DB, ::Type{DBFSComment}, campaign::DBCampaign, dbp::DBPerson, comment::String)
    @assert db == campaign.database "Must write comments to the same database as the campaign."
    write(db, DBFSComment, campaign.pkey, dbp, comment)
end

function Base.read(db::SQLite.DB, ::Type{DBFSComment}, campaign::DBCampaign)::Vector{DBFSComment}
    stmt1 = SQLite.Stmt(db, "SELECT PKEY, CAMPAIGN, PERSON, DATETIME, COMMENT FROM CAMPAIGNCOMMENT WHERE CAMPAIGN=?;")
    q = DBInterface.execute(stmt1, ( campaign.pkey, ))
    res = DBFSComment[]
    while !SQLite.done(q)
        r = SQLite.Row(q)
        @assert r[:CAMPAIGN] == campaign.pkey
        c = DBFSComment(db, r[:PKEY], campaign, read(db, DBPerson, r[:PERSON]), Dates.julian2datetime(r[:DATETIME]), r[:COMMENT])
        push!(res, c)
    end
    return res
end


text:
# Provides a framework on how to translate between LAI coordinates and physical Z-dimensions
#=

using Distributions:Beta,cdf,pdf
using Parameters

abstract type AbstractVerticalCanopyStructure end

"The beta distribution probability for LAI as function of z/h"
Base.@kwdef struct BetaCanopyStructure{FT} <: AbstractVerticalCanopyStructure
    p::FT = 2.5
    q::FT = 2.5
    B = Beta(p,q)
end

# Bonan ML canopy distribution (Beta function):
function convert_z_to_dLAI!(mod::BetaCanopyStructure, canopy, z)
    @unpack LAI, nlayers, h = canopy
    # LAI = Leaf Area Index
    # h = Canopy height
    LAI * cdf.(mod.B,z/h)
end # functi
=#

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function cubemapping(a, b, c=nothing)
# a, b, c are three vectors describing the simulation box
# a = [a1, a2, a3]
# b = [b1, b2, b3]
# c = [c1, c2, c3]

if c === nothing
    C2B = [a[:] b[:]]; # transform the unit square to the simulation box (a, b)
    B2C = inv(C2B); # transform the simulation box (a, b) to the unit square
else
    C2B = [a[:] b[:] c[:]]; # transform the unit cube to the simulation box (a, b, c)
    B2C = inv(C2B); # transform the simulation box (a, b, c) to the unit cube
end

return B2C, C2B;

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# # A vector of SVectors, denoting the edge ID each face is connected to.
# function face_edge_connectivity(mesh::Quadrilateral8Mesh_2D)
#     if length(mesh.edges) === 0
#         @error "Mesh does not have edges!"
#     end
#     # A vector of MVectors of zeros for each face
#     # Each MVector is the length of the number of edges
#     face_edge = [MVector{4, UInt32}(0, 0, 0, 0) for _ in eachindex(mesh.faces)]
#     # For each face in the mesh, generate the edges.
#     # Search for the index of the edge in the mesh.edges vector
#     # Insert the index of the edge into the face_edge connectivity vector
#     for i in eachindex(mesh.faces)
#         for (j, edge) in enumerate(edges(mesh.faces[i]))
#             face_edge[i][j] = searchsortedfirst(mesh.edges, SVector(edge.data))
#         end
#         if any(x->x === 0x00000000, face_edge[i])
#             @error "Could not determine the face/edge connectivity of face $i"
#         end
#     end
#     return [SVector(sort(conn).data) for conn in face_edge]
# end
#
# # A vector of SVectors, denoting the edge ID each face is connected to.
# function face_edge_connectivity(mesh::Triangle6Mesh_2D)
#     if length(mesh.edges) === 0
#         @error "Mesh does not have edges!"
#     end
#     # A vector of MVectors of zeros for each face
#     # Each MVector is the length of the number of edges
#     face_edge = [MVector{3, UInt32}(0, 0, 0) for _ in eachindex(mesh.faces)]
#     # For each face in the mesh, generate the edges.
#     # Search for the index of the edge in the mesh.edges vector
#     # Insert the index of the edge into the face_edge connectivity vector
#     for i in eachindex(mesh.faces)
#         for (j, edge) in enumerate(edges(mesh.faces[i]))
#             face_edge[i][j] = searchsortedfirst(mesh.edges, SVector(edge.data))
#         end
#         if any(x->x === 0x00000000, face_edge[i])
#             @error "Could not determine the face/edge connectivity of face $i"
#         end
#     end
#     return [SVector(sort(conn).data) for conn in face_edge]
# end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Mamba

## Data
bones = Dict{Symbol, Any}(
  :gamma => reshape(
    [ 0.7425,     NaN,     NaN,    NaN,  10.2670,     NaN,     NaN,     NaN,
     10.5215,     NaN,     NaN,    NaN,   9.3877,     NaN,     NaN,     NaN,
      0.2593,     NaN,     NaN,    NaN,  -0.5998,     NaN,     NaN,     NaN,
     10.5891,     NaN,     NaN,    NaN,   6.6701,     NaN,     NaN,     NaN,
      8.8921,     NaN,     NaN,    NaN,  12.4275,     NaN,     NaN,     NaN,
     12.4788,     NaN,     NaN,    NaN,  13.7778,     NaN,     NaN,     NaN,
      5.8374,     NaN,     NaN,    NaN,   6.9485,     NaN,     NaN,     NaN,
     13.7184,     NaN,     NaN,    NaN,  14.3476,     NaN,     NaN,     NaN,
      4.8066,     NaN,     NaN,    NaN,   9.1037,     NaN,     NaN,     NaN,
     10.7483,     NaN,     NaN,    NaN,   0.3887,  1.0153,     NaN,     NaN,
      3.2573,  7.0421,     NaN,    NaN,  11.6273, 14.4242,     NaN,     NaN,
     15.8842, 17.4685,     NaN,    NaN,  14.8926, 16.7409,     NaN,     NaN,
     15.5487, 16.8720,     NaN,    NaN,  15.4091, 17.0061,     NaN,     NaN,
      3.9216,  5.2099,     NaN,    NaN,  15.4750, 16.9406, 17.4944,     NaN,
      0.4927,  1.3556,  2.3016,  3.2535,  1.3059,  1.8793,  2.4970,  3.2306,
      1.5012,  1.8902,  2.3689,  2.9495,  0.8021,  2.3873,  3.9525,  5.3198,
      5.0022,  6.3704,  8.2832, 10.4988,  4.0168,  5.1537,  7.1053, 10.3038],
     4, 34)',
   :delta =>
     [2.9541, 0.6603, 0.7965, 1.0495, 5.7874, 3.8376, 0.6324,
      0.8272, 0.6968, 0.8747, 0.8136, 0.8246, 0.6711, 0.978,
      1.1528, 1.6923, 1.0331, 0.5381, 1.0688, 8.1123, 0.9974,
      1.2656, 1.1802, 1.368,  1.5435, 1.5006, 1.6766, 1.4297,
      3.385,  3.3085, 3.4007, 2.0906, 1.0954, 1.5329],
   :ncat =>
     [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
      3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 5, 5],
   :grade => reshape(
     [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,2,1,1,2,1,1,
      2,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,3,1,1,2,1,1,
      2,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,4,3,3,3,1,1,
      2,1,1,1,2,2,1,1,1,1,1,1,NaN,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,4,5,4,3,1,1,
      2,1,1,1,2,2,1,1,2,1,1,1,1,1,1,1,2,1,1,3,2,1,1,1,1,1,3,1,5,5,5,4,2,3,
      2,1,1,1,2,2,1,2,1,1,1,1,1,2,1,1,2,NaN,1,3,2,1,1,1,1,1,3,1,5,5,5,5,3,3,
      2,1,1,1,2,2,1,1,1,NaN,NaN,1,1,1,1,1,2,NaN,1,3,3,1,1,1,1,1,3,1,5,5,5,5,3,3,
      2,1,2,2,2,2,2,2,1,NaN,NaN,1,2,2,1,1,2,2,1,3,2,1,1,1,1,1,3,1,5,5,5,5,3,4,
      2,1,1,2,2,2,2,2,2,1,1,1,2,1,1,1,2,1,1,3,3,1,1,1,1,1,3,1,5,5,5,5,4,4,
      2,1,2,2,2,2,2,2,2,1,1,1,2,2,2,1,2,NaN,2,3,3,1,1,1,1,1,3,1,5,5,5,5,5,5,
      2,1,NaN,2,2,2,NaN,2,2,1,NaN,NaN,2,2,NaN,NaN,2,1,2,3,3,NaN,1,NaN,1,1,3,1,5,5,5,5,5,5,
      2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,1,NaN,2,1,3,2,5,5,5,5,5,5,
      2,2,2,2,2,2,2,2,2,2,NaN,2,2,2,2,2,2,2,2,3,3,3,NaN,2,NaN,2,3,4,5,5,5,5,5,5],
     34, 13)',
  :nChild => 13,
  :nInd => 34
)

## Model Specification
model = Model(

  grade = Stochastic(2,
    (ncat, delta, theta, gamma, nChild, nInd) ->
      begin
        p = Array{Float64}(5)
        UnivariateDistribution[
          begin
            n = ncat[j]
            p[1] = 1.0
            for k in 1:(n - 1)
              Q = invlogit(delta[j] * (theta[i] - gamma[j, k]))
              p[k] -= Q
              p[k + 1] = Q
            end
            Categorical(p[1:n])
          end
          for i in 1:nChild, j in 1:nInd
        ]
      end,
    false
  ),

  theta = Stochastic(1,
    () -> Normal(0, 100)
  )

)

## Initial Values
inits = [
  Dict(:grade => bones[:grade], :theta => [0.5,1,2,3,5,6,7,8,9,12,13,16,18]),
  Dict(:grade => bones[:grade], :theta => [1,2,3,4,5,6,7,8,9,10,11,12,13])
]

## Sampling Scheme
scheme = [MISS(:grade),
          AMWG(:theta, 0.1)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, bones, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
using PyCall

# TODO: refactor code here and in math_functions.jl.
# Wrap a sympy function with an Symata "function"
# Pass keyword arguments. In Symata, these are expressions with head 'Rule'
# We use this for more than just simplification functions.
# TODO: check number of args, etc.
# TODO: some do not take kwargs. don't waste time looking for them.
macro make_simplify_func(mxprsym, sympyfunc)
    smxprsym = string(mxprsym)[2:end]     # Symata symbol
    ssympyfunc = string(sympyfunc)        # SymPy function
    qsympyfunc = QuoteNode(sympyfunc)
    esc(quote
        function apprules(mx::Mxpr{$mxprsym})
              kws = Dict()
              nargs = sjtopy_kw(mx,kws) # extract keywords from args to mx into kws, return positional args
              (isempty(kws) ? sympy[$qsympyfunc](nargs...) : sympy[$qsympyfunc](nargs...; kws...)) |> pytosj
        end
        set_sysattributes($smxprsym)
        register_sjfunc_pyfunc($smxprsym,$ssympyfunc)
    end)
end

macro make_simplify_func_postp(mxprsym, sympyfunc, postfunc)
    smxprsym = string(mxprsym)[2:end]     # Symata symbol
    ssympyfunc = string(sympyfunc)        # SymPy function
    qsympyfunc = QuoteNode(sympyfunc)
    esc(quote
        function apprules(mx::Mxpr{$mxprsym})
              kws = Dict()
              nargs = sjtopy_kw(mx,kws)
              sres = (isempty(kws) ? sympy[$qsympyfunc](nargs...) : sympy[$qsympyfunc](nargs...; kws...)) |> pytosj
              $postfunc(sres)
        end
        set_sysattributes($smxprsym)
        register_sjfunc_pyfunc($smxprsym,$ssympyfunc)
    end)
end

### Factor

@sjdoc Factor """
    Factor(expr)

factor `expr`. Options are,  `modulus`, `gaussian`, `extension`, and `domain`.

```
modulus => n, gaussian => False, extention => Sqrt(2).
```
"""

@make_simplify_func :Factor factor

#### Expand

@sjdoc Expand """
    Expand(expr)

expand powers and products in `expr`.
"""

@make_simplify_func :Expand expand

#### Limit

@sjdoc Limit """
    Limit(expr, var => lim)

give the `limit` of `expr` as `var` approaches `lim`.
"""

function apprules(mx::Mxpr{:Limit})
    (pymx,var,lim) = map(sjtopy, (mx[1],mx[2][1],mx[2][2]))
    pylimit = sympy[:limit](pymx,var,lim)
    return pytosj(pylimit)
end

register_sjfunc_pyfunc("Limit", "limit")

#### Integrate

@sjdoc Integrate """
    Integrate(expr, x)

gives the indefinite integral of `expr` with respect to `x`.

    Integrate(expr, [x,a,b])

gives the definite integral.
"""

#apprules(mx::Mxpr{:Integrate}) = do_Integrate(mx,margs(mx)...)
# Works for exp with one variable. Is supposed to integrate wrt all vars., but gives error instead.
function do_Integrate(mx::Mxpr{:Integrate},expr)
    pymx = sjtopy(expr)
    pyintegral = sympy[:integrate](pymx)
    return pytosj(pyintegral)
end

function do_Integrate(mx::Mxpr{:Integrate}, expr, varspecs...)
    pymx = sjtopy(expr)
    pyvarspecs = varspecs_to_tuples_of_sympy(collect(varspecs))
    pyintegral = sympy[:integrate](pymx,pyvarspecs...)
    sjres = pytosj(pyintegral)
    sjres
end

function do_Integrate_kws(mx::Mxpr{:Integrate}, kws, expr)
    pymx = sjtopy(expr)
    pyintegral = sympy[:integrate](pymx, kws)
    return pytosj(pyintegral)
end

# We annotate Dict here to fix BoundsError bug in Integrate(x)
function do_Integrate_kws(mx::Mxpr{:Integrate}, kws::T, expr, varspecs...) where T<:Dict
    pymx = sjtopy(expr)
    pyvarspecs = varspecs_to_tuples_of_sympy(collect(varspecs))
    pyintegral = sympy[:integrate](pymx,pyvarspecs...; kws...)
    sjres = pytosj(pyintegral)
    sjres
end

# FIXME: we do deepsetfixed and a symbol Int is returned. If we pull it out,
# it is evaluated to Infinity[1] somehow. Maybe this is positive float Inf
function apprules(mx::Mxpr{:Integrate})
    kws = Dict()
    nargs = separate_rules(mx,kws)
    res = isempty(kws) ? do_Integrate(mx,margs(mx)...) : do_Integrate_kws(mx,kws,nargs...)
    res = list_to_conditional_expression(res)
    fix_integrate_piecewise(typeof(mx),res)
end

register_sjfunc_pyfunc("Integrate", "integrate")

## Sympy returns the unevaluated form as the last member of Piecewise. Mma does not do this.
## In Symata it causes an infinite evaluation loop. So we remove these final forms for Integrate and Sum
fix_integrate_piecewise(y,x) = x
function fix_integrate_piecewise(mxprtype, mx::Mxpr{:Piecewise})
    isempty(mx) && return mx
    length(mx) < 2 && return mx
    last = mx[end]
    isempty(last) && return mx
    if isa(last,Mxpr{:ConditionalExpression})
        isa(last[1], mxprtype) && return mx[1]
    end
    return length(mx) == 1 ? mx[1] : mx
end

list_to_conditional_expression(mx::ListT) = mxpra(:ConditionalExpression, margs(mx))
list_to_conditional_expression(mx) = mx

### LaplaceTransform

@sjdoc LaplaceTransform """
    LaplaceTransform(expr, t, s)

give the Laplace transform of `expr`.

This function returns `[F, a, cond]` where `F` is the Laplace transform of `f`,
`Re(s)>a` is the half-plane of convergence, and `cond` are auxiliary convergence conditions.
(This additional information is currently disabled.)
"""

function apprules(mx::Mxpr{:LaplaceTransform})
    kws = Dict( :noconds => true )
    nargs = sjtopy_kw(mx,kws)
    pyres = @try_sympyfunc sympy[:laplace_transform](nargs...; kws...)  "LaplaceTransform: unknown error."  mx
    res = pyres |> pytosj
    res = list_to_conditional_expression(res)
    fix_integrate_piecewise(typeof(mx),res)
end

#### InverseLaplaceTransform

@sjdoc InverseLaplaceTransform """
   InverseLaplaceTransform(expr, s, t)

gives the inverse Laplace transform of `expr`.
"""

function apprules(mx::Mxpr{:InverseLaplaceTransform})
    result = sympy[:inverse_laplace_transform](map(sjtopy, margs(mx))...)
    sjresult = pytosj(result)
    if isa(sjresult,Mxpr) && mhead(sjresult) == :InverseLaplaceTransform
        setfixed(sjresult)
        if mhead(margs(sjresult)[end]) == :Dummy
            pop!(margs(sjresult)) # we may also want to strip the Dummy()
        end
    end
    sjresult = list_to_conditional_expression(sjresult)
    fix_integrate_piecewise(typeof(mx),sjresult)
end

#### FourierTransform
# TODO, pass options (rules)

@sjdoc FourierTransform """
    FourierTransform(expr, x, k)

gives the Fourier transform of `expr`.

This function returns `[F, cond]` where `F` is the Fourier transform of `f`,
and `cond` are auxiliary convergence conditions.
"""

apprules(mx::Mxpr{:FourierTransform}) = sympy[:fourier_transform](sjtopy(margs(mx))...) |> pytosj

#### InverseFourierTransform

function apprules(mx::Mxpr{:InverseFourierTransform})
    result = sympy[:inverse_fourier_transform](map(sjtopy, margs(mx))...)
    sjresult = pytosj(result)
    if mhead(sjresult) == :InverseFourierTransform
        setfixed(sjresult)
        if mhead(margs(sjresult)[end]) == :Dummy
            pop!(margs(sjresult)) # we may also want to strip the Dummy()
        end
    end
    sjresult = list_to_conditional_expression(sjresult)
    sjresult
end

#### Sum

@sjdoc Sum """
    Sum(expr, [x,a,b])

sums over `x` from `a` to `b`.
"""

apprules(mx::Mxpr{:Sum}) = do_Sum(mx,margs(mx)...)

function do_Sum(mx::Mxpr{:Sum}, expr, varspecs...)
    ## FIXME: only evaluate expr here if the range of the sum is infinite.
    pymx = sjtopy(doeval(expr))
    pyvarspecs = varspecs_to_tuples_of_sympy(reverse(collect(varspecs))) # Symata and Mma use the same convention for position of inner loop
    pysum = sympy[:summation](pymx,pyvarspecs...)
    res = pytosj(pysum)
    if mhead(res) == :Sum
        summand = res[1]
        specs = margs(res)[2:end]
        return mxpr(:Sum,summand,reverse(specs)...)
    end
    fix_integrate_piecewise(typeof(mx),res)
end

#### Product

@sjdoc Product """
    Product(expr, [x,a,b])

the product of `expr` over `x` from `a` to `b`.
"""

apprules(mx::Mxpr{:Product}) = do_Product(mx,margs(mx)...)

function do_Product(mx::Mxpr{:Product}, expr, varspecs...)
    pymx = sjtopy(expr)
    pyvarspecs = varspecs_to_tuples_of_sympy(collect(varspecs))
    pysum = sympy[:product](pymx,pyvarspecs...)
    return pytosj(pysum)
end

register_sjfunc_pyfunc("Product", "product")

#### Series

@sjdoc Series """
    Series(expr,[x,x0,n])

give the Taylor series expansion of `expr`.
"""

apprules(mx::Mxpr{:Series}) = do_Series(mx,margs(mx)...)

function do_Series(mx::Mxpr{:Series}, expr, varspecs...)
    pymx = sjtopy(expr)
    pyspec = []
    for dspec in margs(mx)[2:end]    # Following is more than neccessary. Also, maybe we could use tuples instead of lists
        if isa(dspec,ListT)
            for xdspec in margs(dspec)
                push!(pyspec,sjtopy(xdspec))
            end
        else
            push!(pyspec,sjtopy(dspec))
        end
    end
    pyseries = sympy[:series](pymx,pyspec...)
    return pytosj(pyseries)
end

register_sjfunc_pyfunc("Series", "series")

#### D  (derivative)

## TODO: Implement D and Derivative in Symata. This will still fall through to SymPy
## TODO: In Mma the derivative of a function returns a pure function.

@sjdoc D """
    D(expr, x)

give the partial derivative of `expr` with respect to `x`.

    D(expr,[x,n])

give the `n`th partial derivative.

    D(expr,[x,n1],y,[z,n2])

give the mixed derivative.
"""

function apprules(mx::Mxpr{:D})
    pymx = sjtopy(mx[1])
    pyspec = []
    for dspec in margs(mx)[2:end]  # D(expr, [x,1], y, ...) --> diff(expr,x,1,y,...)
        if isa(dspec,ListT)
            for xdspec in margs(dspec)
                push!(pyspec,sjtopy(xdspec))
            end
        else
            push!(pyspec,sjtopy(dspec))
        end
    end
    pyderivative = sympy[:diff](pymx,pyspec...)
    return pytosj(pyderivative)
end

register_sjfunc_pyfunc("D", "diff")

#### Together

@sjdoc Together """
    Together(sum)

rewrite a sum of terms as a product.
"""

@make_simplify_func :Together together

#### Apart

@sjdoc Apart """
    Apart(product)

compute a partial fraction decomposition of `product`.
"""

apprules(mx::Mxpr{:Apart}) = sympy[:apart](map(sjtopy, margs(mx))...)[:doit]() |> pytosj

register_sjfunc_pyfunc("Apart", "apart")

#### Simplify

@sjdoc Simplify """
    Simplify(expr, kw1 => v1, ...)

rewrites `expr` in a simpler form using keyword options `kw1, ...`.
"""

@make_simplify_func :Simplify simplify
@make_simplify_func :TrigSimp trigsimp
@make_simplify_func :ExpTrigSimp exptrigsimp
@make_simplify_func :RatSimp ratsimp
@make_simplify_func :RadSimp radsimp
@make_simplify_func :PowSimp powsimp
@make_simplify_func :PowDenest powdenest
@make_simplify_func :LogCombine logcombine
@make_simplify_func :ExpandTrig expand_trig
@make_simplify_func :ExpandLog expand_log
@make_simplify_func :SeparateVars separatevars
@make_simplify_func :BesselSimp besselsimp
@make_simplify_func :HyperSimp hypersimp
@make_simplify_func :HyperExpand hyperexpand
@make_simplify_func :NSimplify nsimplify
@make_simplify_func :CombSimp combsimp
@make_simplify_func :SqrtDenest sqrtdenest
@make_simplify_func :Div div

@make_simplify_func_postp :Cse cse _cse_post

function _cse_post(lists)
    (rulelist,expr) = (reverse(margs(lists[1])),lists[2]) # have we defined splat for Mxpr ? we should do it.
    nargs = newargs()
    for x in rulelist
        push!(nargs, mxpr(:Rule, margs(x)))
    end
    rules = mxpr(:List,nargs)
    mxpr(:List, expr, rules)
end

# These apparently have been removed from SymPy
#@make_simplify_func :Separate separate
#@make_simplify_func :OptCse opt_cse  maybe renamed cse_opts
#@make_simplify_func :CollectSqrt collectsqrt apparently gone

@sjdoc TrigSimp """
    TrigSimp(expr)

does trigonometric simplification.
"""

@sjdoc RatSimp """
    RatSimp(expr)

rewrite `expr` with a common denominator, cancel and reduce.
"""

@sjdoc RadSimp """
    RadSimp(expr)

rationalize the denominator.
"""

#### FullSimplify

@sjdoc FullSimplify """
    FullSimplify(expr)

rewrite `expr` in a simpler form, algorithm is more extensive than `Simplify(expr)`,
but likely to be slower.
"""

apprules(mx::Mxpr{:FullSimplify}) = do_FullSimplify(mx)

function do_FullSimplify(mx::Mxpr{:FullSimplify})
    funcs = [sympy[:simplify], sympy[:expand], sympy[:fu], sympy[:powsimp], sympy[:sqrtdenest]]
#    objective = pyeval("lambda x: len(str(x))")
    objective = py"lambda x: len(str(x))"  ## FIXME: this will break now on v0.5. Compatibility
    megasimp = sympy[:strategies][:tree][:greedy]((funcs, funcs), objective)
    mx[1] |> sjtopy |> megasimp |> pytosj
end

@sjdoc Cancel """
    Cancel(expr)

cancel common factors in the numerator and denominator.
"""

apprules(mx::Mxpr{:Cancel}) = mx[1] |> sjtopy |> sympy[:cancel] |> pytosj

register_sjfunc_pyfunc("Cancel", "cancel")

@sjdoc Collect """
    Collect(expr,x)

collect terms involving the same power of `x`.

    Collect(expr,[x,y])

collect terms involving first `x`, then `y`.
"""

@mkapprule Collect

@doap Collect(expr,x) = sympy[:collect](expr |> sjtopy, x |> sjtopy ) |> pytosj
@doap Collect(expr,x,lst::Mxpr{:List}) = sympy[:collect](expr |> sjtopy, x |> sjtopy , list |> sjtopy) |> pytosj

register_sjfunc_pyfunc("Collect", "collect")

#### Solve

@sjdoc Solve """
    Solve(expr)

solve `expr == 0` for one variable.

    Solve(expr,var)

solve `expr == 0` for `var`.

    Solve([expr1,expr2,...], [var1,var2,...])

solve a system of equations.
"""

## TODO: This takes keyword args, in fact a system for sympy keyword args.
##  options => Dict(:Rational => True,
@mkapprule Solve

# TODO: find free symbols in expr and return rules as in following methods
@doap function Solve(expr)
    sres = (expr |> sjtopy |> sympy[:solve] |> pytosj)
#    mxpr(:List, (map(t -> mxpr(:List, mxpr(:Rule,sym,t)), margs(sres)))...)
end

@doap function Solve(expr, var::Symbol)
    pyexpr = expr |> sjtopy
    pyvar = var |> sjtopy
    res =  sympy[:solve](pyexpr,pyvar)
    sres = res |>  pytosj
    mxpr(:List, (mapmargs(t -> mxpr(:List, mxpr(:Rule,var,t)), margs(sres)))...)
end

@doap function Solve(eqs::Mxpr{:List}, vars::Mxpr{:List})
    peqs = eqs |> sjtopy
    pyvars = vars |> sjtopy
    sres = sympy[:solve](peqs,pyvars) |>  pytosj
    if isa(sres,Dict)  # why does pytosj sometimes return Dict and sometimes List ?
        nargs = newargs()
        for (k,v) in sres
            push!(nargs, mxpr(:Rule, k, v))
        end
        mxpr(:List, mxpr(:List, nargs...))
    else
        return sres  # TODO. convert these to Rules
    end
end

register_sjfunc_pyfunc("Solve", "solve")

# This is broken
# apprules(mx::Mxpr{:DSolve}) = do_DSolve(mx,margs(mx)...)
# do_DSolve(mx, expr) = expr |> sjtopy |> sympy[:dsolve] |> pytosj
@mkapprule DSolve

@doap function DSolve(expr, fcn)
    pyexpr = expr |> sjtopy
    pyfcn = fcn |> sjtopy
    res =  sympy[:dsolve](pyexpr,pyfcn)
    sres = res |>  pytosj
#    mxpr(:List, (map(t -> mxpr(:List, mxpr(:Rule,var,t)), margs(sres)))...)
end

#### Roots
@sjdoc Roots """
    Roots(expr)

solves for the roots of `expr`.

returns a `List` of `Lists`. The two elements of each sublist give the root and its multiplicity.
"""

apprules(mx::Mxpr{:Roots}) = mx[1] |> sjtopy |> sympy[:roots] |> pytosj  |> Symata.unpacktoList

register_sjfunc_pyfunc("Roots", "roots")

### RealRoots
@sjdoc RealRoots """
    RealRoots(expr)

solves for the real roots of `expr`.
"""

apprules(mx::Mxpr{:RealRoots}) = mx[1] |> sjtopy |> sympy[:real_roots] |> pytosj

register_sjfunc_pyfunc("RealRoots", "real_roots")

### Singularities

@mkapprule Singularities :nargs => 2

@doap function Singularities(expr,var)
    res = pytosj(sympy[:singularities](sjtopy(expr),sjtopy(var)))
    if isa(res,Mxpr{:EmptySet})
        return List()
    elseif isa(res,Mxpr{:FiniteSet})
        mxpr(:List, mapmargs( x -> mxpr(:Rule, var, x), margs(res))...)
    else
        return res
    end
end

### TODO: IsIncreasing, etc.

#### ToSymPy

@sjdoc ToSymPy """
    ToSymPy(expr)

convert `expr` to a (python) PyObject.
"""

function apprules(mx::Mxpr{:ToSymPy})
    res = sjtopy(margs(mx)...)
end

#### ToSymata

@sjdoc ToSymata """
    ToSymata(expr)

convert the python PyObject `expr` to a Symata expression. Normally, expressions computed
by SymPy are automatically converted to Symata expressions.
"""

apprules(mx::Mxpr{:ToSymata}) = do_ToSymata(mx,margs(mx)...)

do_ToSymata(mx::Mxpr{:ToSymata}, s::Symbol) = setfixed(_pytosj(symval(s)))
do_ToSymata(mx::Mxpr{:ToSymata}, e::PyCall.PyObject) = setfixed(_pytosj(e))
do_ToSymata(mx::Mxpr{:ToSymata}, x) = x

#### PossibleClosedForm

apprules(mx::Mxpr{:PossibleClosedForm}) = do_PossibleClosedForm(mx,margs(mx)...)

@sjdoc PossibleClosedForm """
    PossibleClosedForm(x)

attempts to find an exact formula for the floating point number `x`.
"""

do_PossibleClosedForm(mx::Mxpr{:PossibleClosedForm},x::AbstractFloat) =  x |>  mpmath.identify |> sympy[:sympify] |> pytosj

function do_PossibleClosedForm(mx::Mxpr{:PossibleClosedForm}, args...)
    kws = Dict()
    nargs = sjtopy_kw(mx,kws)
    pyargs = sjtopy(nargs...)
    setdps::Bool = false
    if haskey(kws,"dps")
        set_mpmath_dps(kws["dps"])
        delete!(kws, "dps")
        setdps = true
    end
    pyres = mpmath.identify(pyargs...; kws...)
    if setdps restore_mpmath_dps() end
    pyres |> sympy[:sympify] |> pytosj
end

##### ConditionalExpression

# There are several other properties to implement here.

@sjdoc ConditionalExpression """
    ConditionalExpression(expr,cond)

returns `expr` only when `cond` is true. If `cond` is false,
return `Undefined`.
"""

@mkapprule ConditionalExpression :nargs => 2
@doap ConditionalExpression(expr, cond::Bool) = cond ? expr : Undefined
@doap ConditionalExpression(expr, cond) = mx

### Refine

@mkapprule Refine  :nodefault => true

@sjdoc Refine """
    Refine(expr)

simplifies expr using assumptions. For instance, `Assume(x,positive)`.
"""

@doap function Refine(args...)
    result = sympy[:refine](map(sjtopy, args)...)
    result |> pytosj
end

### CoefficientList

@mkapprule CoefficientList

@doap function CoefficientList(expr,var)
    (spexpr,spvar) = (sjtopy(expr),sjtopy(var))
    p = sympy[:Poly](spexpr,spvar)
#    cs = p[:coeffs]()  # lists only non-zero coefficients
    cs = reverse!(p[:all_coeffs]())
    List(mapmargs(pytosj,cs))
end

### CoefficientSympy

@mkapprule CoefficientSympy

@doap function CoefficientSympy(expr,subexpr)
    (spexpr,spsubexpr) = (sjtopy(expr),sjtopy(subexpr))
    cs = spexpr[:coeff](spsubexpr)
    pytosj(cs)
end

@doap function CoefficientSympy(expr,subexpr,pow)
    (spexpr,spsubexpr) = (sjtopy(expr),sjtopy(mpow(subexpr,pow)))
    cs = spexpr[:coeff](spsubexpr)
    pytosj(cs)
end

## utility

# input -- Array of Symata Lists and/or Symbols
# output -- Array of tuples (from Lists) of SymPy objects, or single SymPy objects
# Eg: For translating Integrate(expr,[x,a,b],y) --> integrate(expr,(x,a,b),y)
function varspecs_to_tuples_of_sympy(args::Array)
    oarr = []
    for x in args
        push!(oarr, isa(x,ListT) ? tuple(map(sjtopy, margs(x))...) : sjtopy(x))
        # if isa(x,ListT)
        #     push!(oarr, tuple(map(sjtopy, margs(x))...))
        # else
        #     push!(oarr,sjtopy(x))
        # end
    end
    return oarr
end


text:
using Test

using CompScienceMeshes
using StaticArrays

T = Float64
P = SVector{3,T}

width, height = 1.0, 1.0
mesh = meshrectangle(width, height, 0.5)

# This predicate tests whether a simplex is in the (x==0) plane
function pred(simplex)
    verts = mesh.vertices[simplex]
    for i in 1 : length(verts)
        if abs(verts[i][1]) > eps(Float64)
            return false
        end
    end
    return true
end

# Test if meshes can be intersected
line = meshsegment(width, 1/3, 3)
translate!(line, P(0.0, height, 0.0))

edges = skeleton(mesh, 1)
γ1 = submesh(line, edges)
@test numcells(γ1) == 2
for c in cells(γ1)
    ctr = cartesian(center(chart(edges, c)))
    @test ctr[2] == height
    @test ctr[3] == 0
end
# @test γ1.faces == [
#     index(3,6),
#     index(6,9)]

pred2 = interior_tpredicate(mesh)
γ2 = submesh(pred2, edges)
@test numcells(γ2) == 8

overlaps = overlap_gpredicate(line)
pred3 = x -> overlaps(simplex(vertices(mesh, x)))
pred4 = x -> pred2(x) || pred3(x)
γ3 = submesh(pred4, edges)
@test numcells(γ3) == 10

#test interior_vpredicate
function boundaryvertices(mesh::Mesh)
    vpred = interior_vpredicate(mesh)
    bidx = findall(!vpred,1:numvertices(mesh))
    bnd = mesh.vertices[bidx]
end

tr = meshrectangle(1.0,1.0,1/4)
bv = boundaryvertices(tr)
btr = boundary(tr)
@test length(bv) == numcells(skeleton(btr,0))


text:
module FSimPlots

using FSimZoo
using Plots
import Plots: plot, plot!
using Transducers
using UnPack
using LinearAlgebra

include("multicopters.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ d8a922a0-e5b0-11eb-3290-75b867a2f730
begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
end

# ╔═╡ 5642cd29-c2ab-4943-985e-64c5999f7bae
using
    DiffEqFlux, DifferentialEquations, GalacticOptim,
    Flux, Optim, BlackBoxOptim,
    Zygote, ForwardDiff,
    Distributions, Statistics,
    Distances,
    KernelDensity,
    CUDA,
    Logging, TerminalLoggers,
    StatsPlots,
    DataFrames

# ╔═╡ 0085e83f-a3df-42f8-8a0e-658ec2afca50
begin
    global_logger(TerminalLogger())
    plotly()
end

# ╔═╡ 3f8520e1-89f4-4b21-8e52-7bf48c760db5
data_dist = Beta(2, 4)

# ╔═╡ 59e766f3-c011-4aa3-a42d-f0cf03d1c00a
plot(data_dist, label=string(data_dist))

# ╔═╡ 7fb8118e-26a2-4d23-adb0-d20d3d024ff6
train_data = rand(data_dist, 1, 100)

# ╔═╡ 5cf5412d-7eb0-410d-897a-a8d0e273b21f
test_data = rand(data_dist, 1, 1000)

# ╔═╡ 74801cce-a47b-4f60-a8c8-f56bcdc2534d
mdl_dist_f_kde = KernelDensity.kernel_dist(Normal, KernelDensity.default_bandwidth(reshape(train_data, size(train_data, 2))))

# ╔═╡ 1cc754a6-385d-4f87-9f06-9d3cb8bb1a3c
mdl_dist_f_data = Normal(mean(train_data), std(train_data))

# ╔═╡ d8869f7a-988c-4dc1-9f8f-6467ff1657c5
mdl_dist_f_dist = Normal(mean(data_dist), std(data_dist))

# ╔═╡ f7563269-4749-41c7-b4a8-1fd8f9b73c14
let
    p = plot(data_dist, label=string(data_dist))
    p = plot!(p, mdl_dist_f_dist, label=string(mdl_dist_f_dist))
    p = plot!(p, mdl_dist_f_data, label=string(mdl_dist_f_data))
    p = plot!(p, mdl_dist_f_kde, label=string(mdl_dist_f_kde))
end

# ╔═╡ d7175898-e654-4f0f-885b-d0e104e62507
n_hidden = 3

# ╔═╡ 929715b5-3e57-4403-a081-3fd6499d916d
actv = tanh

# ╔═╡ fa1edbc5-a05c-463f-83d9-8a4cc492f731
nn = Chain(
    Dense(1, n_hidden, actv),
    # Dense(n_hidden, n_hidden, actv),
    Dense(n_hidden, 1, actv),
) |> f64

# ╔═╡ 35769a13-0537-4413-b5f9-4a7a126bda3c
tspan = (0.0, 1.0)

# ╔═╡ a65312c6-a9c3-4a0f-ba88-afe122e30395
ffjord_mdl = FFJORD(nn, tspan, Tsit5(), basedist=mdl_dist_f_dist)

# ╔═╡ 595b19d9-de04-4038-a87e-316ada76e070
ffjord_mdl.p

# ╔═╡ daab25dd-6a54-46fa-a26c-d2f05183659e
ffjord_mdl.basedist

# ╔═╡ 4bb079bb-fc0a-476e-bb42-af1caa1805f2
plot(ffjord_mdl.basedist)

# ╔═╡ 763966b9-a42c-4787-8e5f-c925a4b44f15
regularize = false

# ╔═╡ ac7cccd8-e564-4b4a-bf10-bcea49d2ff18
monte_carlo = false

# ╔═╡ 30a83fcc-cb24-4033-a665-92869215da78
function loss(θ)
    logpx, λ₁, λ₂ = ffjord_mdl(train_data, θ; regularize, monte_carlo)
    -mean(logpx)
end

# begin
# 	actual_pdf = pdf.(data_dist, train_data)
# 	function loss(p)
# 		learned_pdf = exp.(first(ffjord_mdl(train_data, p; regularize, monte_carlo)))
# 		dis = totalvariation(learned_pdf, actual_pdf) / size(train_data, 2)
# 		-dis
# 	end
# end

# ╔═╡ 842fc843-2a89-445c-8319-50a5fafec73b
adtype = GalacticOptim.AutoZygote()
# adtype = GalacticOptim.AutoForwardDiff()

# ╔═╡ 4b77653f-1cbb-4d7c-8dc8-6ec0e91b4b0d
optm1 = Flux.ADAM(0.1)
# optm1 = Optim.GradientDescent()

# ╔═╡ 9686dc0b-1f5e-4ce0-958c-7059caf93e70
optm2 = Optim.BFGS()

# ╔═╡ 32fc54a0-3bd7-4ab6-a5de-93cb6d1c3e2e
optm3 = BBO(:adaptive_de_rand_1_bin)

# ╔═╡ f56bfb88-8ac1-4227-a0b1-60eb613aec67
maxiters_m1 = 100

# ╔═╡ 114ce7be-ceaa-4efd-8890-95cc95cf3219
maxiters_m2 = 25

# ╔═╡ 4e42aa35-dfae-422e-ac8c-13c6b80efba4
maxiters_m3 = 10

# ╔═╡ bb213e1d-3ece-40de-8007-4fd12a814823
progress = true

# ╔═╡ 7a6d80ad-2d15-42f0-8012-6316f9905ecc
res1 = DiffEqFlux.sciml_train(loss, ffjord_mdl.p, optm1, adtype; progress, maxiters=maxiters_m1)

# ╔═╡ 83018b3a-27e6-42ae-9dc3-0f684022c386
res1.u

# ╔═╡ 987e0cee-36ba-4130-8b7d-65a8bcfd9bb3
res2 = DiffEqFlux.sciml_train(loss, res1.u, optm2, adtype; progress, maxiters=maxiters_m2)
# res2 = DiffEqFlux.sciml_train(loss, res1.u, optm2, adtype; progress)

# ╔═╡ fd432df9-d795-4732-806f-e0e910c3599d
res2.u

# ╔═╡ 7b40ebee-3559-4cfc-a823-d3772177acc3
# res3 = DiffEqFlux.sciml_train(loss, res2.u, optm3, adtype; progress, maxiters=maxiters_m3)

# ╔═╡ 5ed15a42-3ace-4850-99ac-f4ce705034cf
final_res = res2

# ╔═╡ 60e40d80-eee0-4a08-a679-6727f9a91852
loss(final_res.u)

# ╔═╡ 4f4b8781-5eca-442b-8950-b94f866351db
train_dis = let data = train_data
    actual_pdf = pdf.(data_dist, data)
    learned_pdf = exp.(first(ffjord_mdl(data, final_res.u; regularize, monte_carlo)))
    dis = totalvariation(learned_pdf, actual_pdf) / size(data, 2)
end

# ╔═╡ f97042ce-980d-433b-a293-63b84b01c3e9
test_dis = let data = test_data
    actual_pdf = pdf.(data_dist, data)
    learned_pdf = exp.(first(ffjord_mdl(data, final_res.u; regularize, monte_carlo)))
    dis = totalvariation(learned_pdf, actual_pdf) / size(data, 2)
end

# ╔═╡ d9225bf7-d4b6-4462-bd51-0639ff90d0df
begin
    hist_range = range(0, 1, length=1000 + 1)
    hist_data = reshape(collect(hist_range), (1, size(hist_range, 1)))
end

# ╔═╡ ebf17afb-9e1d-4465-a337-603309ef345d
begin
    hist_pdf_0 = exp.(first(ffjord_mdl(hist_data; regularize, monte_carlo)))
    hist_pdf_1 = exp.(first(ffjord_mdl(hist_data, res1.u; regularize, monte_carlo)))
    hist_pdf_2 = exp.(first(ffjord_mdl(hist_data, res2.u; regularize, monte_carlo)))
end

# ╔═╡ faf10042-af69-46d2-a891-7e774fbcea63
kde_mdl = kde(reshape(train_data, size(train_data, 2)))

# ╔═╡ afbc4978-2eee-41fb-a91a-f89fc9989b4a
    kde_train_dis = let data = train_data
    actual_pdf = pdf.(data_dist, data)
    learned_pdf = reshape(pdf(kde_mdl, reshape(data, size(data, 2))), size(data))
    dis = totalvariation(learned_pdf, actual_pdf) / size(data, 2)
end

# ╔═╡ d8c2056b-6af3-42f5-8ce6-dfb3cd8b638f
kde_test_dis = let data = test_data
    actual_pdf = pdf.(data_dist, data)
    learned_pdf = reshape(pdf(kde_mdl, reshape(data, size(data, 2))), size(data))
    dis = totalvariation(learned_pdf, actual_pdf) / size(data, 2)
end

# ╔═╡ d9afd296-82a9-4989-bcf8-e08ec01e25c2
kde_hist_pdf = reshape(pdf(kde_mdl, reshape(hist_data, size(hist_data, 2))), size(hist_data))

# ╔═╡ 38961acf-e62d-4abf-9e21-0f2e5646506a
let data = hist_data
    p = plot(data_dist, label=string(data_dist))
    p = plot!(p, reshape(data, size(data, 2)), reshape(hist_pdf_0, size(data, 2)), label="p(x) - init")
    p = plot!(p, reshape(data, size(data, 2)), reshape(hist_pdf_1, size(data, 2)), label="p(x) - after adam")
    p = plot!(p, reshape(data, size(data, 2)), reshape(hist_pdf_2, size(data, 2)), label="p(x) - after bfgs")
    p = plot!(p, reshape(data, size(data, 2)), reshape(kde_hist_pdf, size(data, 2)), label="p(x) - kde")
end

# ╔═╡ 3e187938-77e0-4fef-b0ee-b04dd0fee97b
result = DataFrame(
    Method=["FFJORD", "FFJORD", "KDE", "KDE"],
    Data=["train_data", "test_data", "train_data", "test_data"],
    TV_Error=[train_dis, test_dis, kde_train_dis, kde_test_dis],
)

# ╔═╡ 537a31bc-360c-4231-9d35-8e899d642126
let data = train_data
    p = plot(data_dist, label=string(data_dist))
    p = scatter!(p, reshape(data, size(data, 2)), reshape(pdf.(data_dist, data), size(data, 2)), label="sample for training")
    p = plot!(p, reshape(hist_data, size(hist_data, 2)), reshape(hist_pdf_2, size(hist_data, 2)), label="p(x) - final")
end

# ╔═╡ Cell order:
# ╠═d8a922a0-e5b0-11eb-3290-75b867a2f730
# ╠═5642cd29-c2ab-4943-985e-64c5999f7bae
# ╠═0085e83f-a3df-42f8-8a0e-658ec2afca50
# ╠═3f8520e1-89f4-4b21-8e52-7bf48c760db5
# ╠═59e766f3-c011-4aa3-a42d-f0cf03d1c00a
# ╠═7fb8118e-26a2-4d23-adb0-d20d3d024ff6
# ╠═5cf5412d-7eb0-410d-897a-a8d0e273b21f
# ╠═74801cce-a47b-4f60-a8c8-f56bcdc2534d
# ╠═1cc754a6-385d-4f87-9f06-9d3cb8bb1a3c
# ╠═d8869f7a-988c-4dc1-9f8f-6467ff1657c5
# ╠═f7563269-4749-41c7-b4a8-1fd8f9b73c14
# ╠═d7175898-e654-4f0f-885b-d0e104e62507
# ╠═929715b5-3e57-4403-a081-3fd6499d916d
# ╠═fa1edbc5-a05c-463f-83d9-8a4cc492f731
# ╠═35769a13-0537-4413-b5f9-4a7a126bda3c
# ╠═a65312c6-a9c3-4a0f-ba88-afe122e30395
# ╠═595b19d9-de04-4038-a87e-316ada76e070
# ╠═daab25dd-6a54-46fa-a26c-d2f05183659e
# ╠═4bb079bb-fc0a-476e-bb42-af1caa1805f2
# ╠═763966b9-a42c-4787-8e5f-c925a4b44f15
# ╠═ac7cccd8-e564-4b4a-bf10-bcea49d2ff18
# ╠═30a83fcc-cb24-4033-a665-92869215da78
# ╠═842fc843-2a89-445c-8319-50a5fafec73b
# ╠═4b77653f-1cbb-4d7c-8dc8-6ec0e91b4b0d
# ╠═9686dc0b-1f5e-4ce0-958c-7059caf93e70
# ╠═32fc54a0-3bd7-4ab6-a5de-93cb6d1c3e2e
# ╠═f56bfb88-8ac1-4227-a0b1-60eb613aec67
# ╠═114ce7be-ceaa-4efd-8890-95cc95cf3219
# ╠═4e42aa35-dfae-422e-ac8c-13c6b80efba4
# ╠═bb213e1d-3ece-40de-8007-4fd12a814823
# ╠═7a6d80ad-2d15-42f0-8012-6316f9905ecc
# ╠═83018b3a-27e6-42ae-9dc3-0f684022c386
# ╠═987e0cee-36ba-4130-8b7d-65a8bcfd9bb3
# ╠═fd432df9-d795-4732-806f-e0e910c3599d
# ╠═7b40ebee-3559-4cfc-a823-d3772177acc3
# ╠═5ed15a42-3ace-4850-99ac-f4ce705034cf
# ╠═60e40d80-eee0-4a08-a679-6727f9a91852
# ╠═4f4b8781-5eca-442b-8950-b94f866351db
# ╠═f97042ce-980d-433b-a293-63b84b01c3e9
# ╠═d9225bf7-d4b6-4462-bd51-0639ff90d0df
# ╠═ebf17afb-9e1d-4465-a337-603309ef345d
# ╠═faf10042-af69-46d2-a891-7e774fbcea63
# ╠═afbc4978-2eee-41fb-a91a-f89fc9989b4a
# ╠═d8c2056b-6af3-42f5-8ce6-dfb3cd8b638f
# ╠═d9afd296-82a9-4989-bcf8-e08ec01e25c2
# ╠═38961acf-e62d-4abf-9e21-0f2e5646506a
# ╠═3e187938-77e0-4fef-b0ee-b04dd0fee97b
# ╠═537a31bc-360c-4231-9d35-8e899d642126

max_stars_repo_name: prbzrg/FFJORD-Test

max_issues_repo_name: prbzrg/FFJORD-Test

max_forks_repo_name: prbzrg/FFJORD-Test


text:
@testset "Eigen_3D" begin
    all_types = (RotMatrix{3}, AngleAxis, RotationVec,
                 UnitQuaternion, RodriguesParam, MRP,
                 RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
                 RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ,
                 RotX, RotY, RotZ,
                 RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY)
    oneaxis_types = (RotX, RotY, RotZ)

    @testset "$(T)" for T in all_types, F in (one, rand)
        R = F(T)
        λs = eigvals(R)
        vs = eigvecs(R)
        E = eigen(R)
        v1 = vs[:,1]
        v2 = vs[:,2]
        v3 = vs[:,3]
        @test R * vs ≈ transpose(λs) .* vs
        @test norm(v1) ≈ 1
        @test norm(v2) ≈ 1
        @test norm(v3) ≈ 1
        @test E.values == λs
        @test E.vectors == vs
        if !(T in oneaxis_types) && VERSION ≥ v"1.2"
            # If the rotation angle is in [0°, 180°], then the eigvals will be equal.
            # Note that the randomized RotX (and etc.) have rotation angle in [0°, 360°].
            # This needs Julia(≥1.2) to get sorted eigenvalues in a canonical order
            # See https://github.com/JuliaLang/julia/pull/21598
            @test eigvals(R) ≈ eigvals(collect(R))
        end
    end
end

@testset "Eigen_2D" begin
    all_types = (RotMatrix{2}, Angle2d)

    @testset "$(T)" for T in all_types, θ in 0.0:0.1:π
        R = T(Angle2d(θ))
        λs = eigvals(R)
        vs = eigvecs(R)
        E = eigen(R)
        v1 = vs[:,1]
        v2 = vs[:,2]
        @test R * vs ≈ transpose(λs) .* vs
        @test norm(v1) ≈ 1
        @test norm(v2) ≈ 1
        @test E.values == λs
        @test E.vectors == vs
        if VERSION ≥ v"1.2"
            # If the rotation angle θ is in [0°, 180°], then the eigvals will be equal.
            # Note that the randomized RotX (and etc.) have rotation angle in [0°, 360°].
            # This needs Julia(≥1.2) to get sorted eigenvalues in a canonical order
            # See https://github.com/JuliaLang/julia/pull/21598
            @test eigvals(R) ≈ eigvals(collect(R))
        end
    end
end


text:
using FastTransforms, Test

@testset "Chebyshev transform"  begin
    @testset "Chebyshev points" begin
        @test chebyshevpoints(10; kind=1) == chebyshevpoints(Float64, 10; kind=1)
        @test chebyshevpoints(10; kind=2) == chebyshevpoints(Float64, 10; kind=2)
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            @test chebyshevpoints(T, 0, kind=1) == chebyshevpoints(T, 0, kind=2) == T[]
            @test chebyshevpoints(T, 1, kind=1) == chebyshevpoints(T, 1, kind=2) == T[0]

            n = 20
            @inferred(chebyshevpoints(T, n, kind=1))
            @inferred(chebyshevpoints(T, n, kind=2))
            @test_throws ArgumentError chebyshevpoints(n, kind=-1)
        end
    end

    @testset "Chebyshev first kind points <-> first kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_1 = chebyshevpoints(T, n, kind=1)
            f = exp.(p_1)
            f̌ = chebyshevtransform(f; kind=1)

            f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = plan_chebyshevtransform(f; kind=1)
            @test P*f == f̌
            @test f == f̃
            P = plan_chebyshevtransform!(f; kind=1)
            @test P*f == f̌
            @test f == f̌
            Pi = plan_ichebyshevtransform(f̌; kind=1)
            @test Pi*f̌ ≈ f̃
            @test f̌ == f̄
            Pi = plan_ichebyshevtransform!(f̌; kind=1)
            @test Pi*f̌ ≈ f̃
            @test f̌ ≈ f̃

            @test chebyshevtransform(T[1]; kind=1) == T[1]
            @test ichebyshevtransform(T[1]; kind=1) == T[1]
            @test chebyshevtransform(T[]; kind=1) == T[]
            @test ichebyshevtransform(T[]; kind=1) == T[]
        end
    end
    @testset "Chebyshev second kind points <-> first kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_2 = chebyshevpoints(T, n, kind=2)
            f = exp.(p_2)
            f̌ = chebyshevtransform(f; kind=2)

            f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)

            P = plan_chebyshevtransform!(f; kind=2)
            Pi = plan_ichebyshevtransform!(f; kind=2)
            @test all((P \ copy(f)) .=== Pi * copy(f))
            @test all((Pi \ copy(f̌)) .=== P * copy(f̌))
            @test f ≈ P \ (P*copy(f)) ≈ P * (P\copy(f)) ≈ Pi \ (Pi*copy(f)) ≈ Pi * (Pi \ copy(f))

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = plan_chebyshevtransform(f; kind=2)
            @test P*f == f̌
            @test f == f̃
            P = plan_chebyshevtransform!(f; kind=2)
            @test P*f == f̌
            @test f == f̌
            Pi = plan_ichebyshevtransform(f̌; kind=2)
            @test Pi*f̌ ≈ f̃
            @test f̌ == f̄
            Pi = plan_ichebyshevtransform!(f̌; kind=2)
            @test Pi*f̌ ≈ f̃
            @test f̌ ≈ f̃

            @test chebyshevtransform(T[1]; kind=2) == T[1]
            @test ichebyshevtransform(T[1]; kind=2) == T[1]
            @test chebyshevtransform(T[]; kind=2) == T[]
            @test ichebyshevtransform(T[]; kind=2) == T[]
        end
    end

    @testset "Chebyshev first kind points <-> second kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_1 = chebyshevpoints(T, n, kind=1)
            f = exp.(p_1)
            f̌ = chebyshevutransform(f; kind=1)

            f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test ichebyshevutransform(f̌; kind=1) ≈ exp.(p_1)

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = plan_chebyshevutransform(f; kind=1)
            @test P*f == f̌
            @test f == f̃
            P = plan_chebyshevutransform!(f; kind=1)
            @test P*f == f̌
            @test f == f̌
            Pi = plan_ichebyshevutransform(f̌; kind=1)
            @test Pi*f̌ ≈ f̃
            @test f̌ == f̄
            Pi = plan_ichebyshevutransform!(f̌; kind=1)
            @test Pi*f̌ ≈ f̃
            @test f̌ ≈ f̃

            @test chebyshevutransform(T[1]; kind=1) == T[1]
            @test ichebyshevutransform(T[1]; kind=1) == T[1]
            @test chebyshevutransform(T[]; kind=1) == T[]
            @test ichebyshevutransform(T[]; kind=1) == T[]
        end
    end

    @testset "Chebyshev second kind points <-> second kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_2 = chebyshevpoints(T, n, kind=2)[2:end-1]
            f = exp.(p_2)
            f̌ = chebyshevutransform(f; kind=2)

            f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test ichebyshevutransform(f̌; kind=2) ≈ exp.(p_2)

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = plan_chebyshevutransform(f; kind=2)
            @test P*f == f̌
            @test f == f̃
            P = plan_chebyshevutransform!(f; kind=2)
            @test P*f == f̌
            @test f == f̌
            Pi = plan_ichebyshevutransform(f̌; kind=2)
            @test Pi*f̌ ≈ f̃
            @test f̌ == f̄
            Pi = plan_ichebyshevutransform!(f̌; kind=2)
            @test Pi*f̌ ≈ f̃
            @test f̌ ≈ f̃

            @test chebyshevutransform(T[1]; kind=2) == T[1]
            @test ichebyshevutransform(T[1]; kind=2) == T[1]
            @test chebyshevutransform(T[]; kind=2) == T[]
            @test ichebyshevutransform(T[]; kind=2) == T[]
        end
    end
end


text:
struct TaylorGreenVortex <: FluidFlowProblem
    rho_0::Float64
    u_max::Float64
    ν::Float64
    NX::Int64
    NY::Int64
    domain_size::Tuple{Float64, Float64}
    static::Bool
    A::Float64
    B::Float64
    a::Float64
    b::Float64
end
function TaylorGreenVortex(
    ν = 1.0 / 6.0,
    scale = 2,
    NX = 16 * scale,
    NY = NX,
    domain_size = (2pi, 2pi);
    static = true,
    A = 1,
    B = -1,
    a = 1,
    b = 1,
)
    u_max = 0.01 / scale
    Re = NX * u_max / ν

    p = TaylorGreenVortex(1.0, u_max, ν, NX, NY, domain_size, static, A, B, a, b)
    if (p.u_max > sqrt(2 / 3) * delta_x(p) / delta_t(p))
        @warn p.u_max, sqrt(2 / 3) * delta_x(p) / delta_t(p)
    end
    p
end

function density(
    q::Quadrature,
    problem::TaylorGreenVortex,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    return pressure(q, problem, x, y, timestep)

    # If not athermal
    return 1.0
end

function pressure(
    q::Quadrature,
    problem::TaylorGreenVortex,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    # return 1.0
    a = problem.a
    A = problem.A
    b = problem.b
    B = problem.B

    P =
        -(1 / 4) *
        problem.rho_0 *
        decay(problem, x, y, timestep)^2 *
        (A^2 * cos(2 * a * x) + B^2 * cos(2 * b * y))

    return 1.0 + q.speed_of_sound_squared * problem.u_max^2 * P
end

function velocity(
    problem::TaylorGreenVortex,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    a = problem.a
    A = problem.A
    b = problem.b
    B = problem.B
    return decay(problem, x, y, timestep) * [
        A * cos(a * x) * sin(b * y)
        B * sin(a * x) * cos(b * y)
    ]
end

function velocity_gradient(
    problem::TaylorGreenVortex,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    a = problem.a
    A = problem.A
    b = problem.b
    B = problem.B

    u_x = -a * A * sin(a * x) * sin(b * y)
    v_y = -b * B * sin(a * x) * sin(b * y)
    u_y = b * A * cos(a * x) * cos(b * y)
    v_x = a * B * cos(a * x) * cos(b * y)

    return decay(problem, x, y, timestep) * [u_x v_x; u_y v_y]
end

function decay(problem::TaylorGreenVortex, x::Float64, y::Float64, timestep::Float64)
    a = problem.a
    b = problem.b

    return problem.static ? 1.0 : exp(-(a^2 + b^2) * viscosity(problem) * timestep)
end

function force(problem::TaylorGreenVortex, x::Float64, y::Float64, time::Float64 = 0.0)
    return problem.static ? 2 * viscosity(problem) * velocity(problem, x, y, 0.0) :
           [0.0 0.0]
end

has_external_force(problem::TaylorGreenVortex) = problem.static

# t_c(problem) = 1
# t_c(problem::TaylorGreenVortex) = ln(2) / (problem.ν * (problem.k_x^2 + problem.k_y^2))


text:
# A BonesaTuner has state for optimizing a certain algorithm on a set of
# parameters described by a ParameterSet. The parameter vectors that are tuned
# only take values in [0.0, 1.0] and it is up to the ParameterSet to translate
# that to meaningful ranges.
type BonesaTuner
  numvectors::Int64                     # Current num of param vectors in archive
  parameter_vectors::Array{Float64, 2}  # Archive of param vectors, each column is one vector
  utilities::Array{Float64, 2}          # Each column has Nu * Np utility values for one and the same parameter vector
  num_utility_values
end

# The predicted utility of a parameter vector x on problem is a weigthed, kernel
# smoothed average of the utilities of the vectors in the archive.
function predicted_utility(ba::BonesaArchive, x, problem)
end

function time_is_up(bt::BonesaTuner)
  time() - bt.start_time >= bt.max_time_seconds
end

# Start tuning parameters given a maximum time budget in hours.
function tune_parameters(bt::BonesaTuner, max_hours = 1.0)
  # Set up for measuring time left.
  bt.start_time = time()
  bt.max_time_seconds = max_hours * 60.0 * 60.0

  # Create a random set of vectors and then evaluate them on the problem set.
  create_random_parameter_vectors(bt, bt.M)
  bt.utilities = zeros(num_utilities(bt) * num_problems(bt.problem_set), bt.M)
  for(i in 1:num_param_vectors(bt))
    if time_is_up(bt)
      # Cut the size of the archive since we did not have time to eval them all... :(
      bt.parameter_vectors = bt.parameter_vectors[:,1:(i-1)]
      bt.utilities = bt.utilities[:,1:(i-1)]
      break
    end
    bt.utilities[:, i] = evaluate_param_vector(bt, bt.parameter_vectors[:,i])
  end

  # Iteratively add to archive as we learn more about which param vectors are good.
  while( !time_is_up(bt) )

    # Update constants and temp calcs needed to calc pareto strenght below

  end

  # Select the terminal set (pareto front) of the param vectors.

end

# Calc predicted utilities, variance and support as well as pareto dominance
# and strengths for all param vectors in the archive.
function calc_temp_values_for_archive(bt::BonesaTuner)
  # Update magic constant c based on size of archive
  bt.c = calc_approximate_magic_constant_c(size_of_archive(bt), numdims(bt))

  # Calc the good vectors (saved as (boolean) indicator values)
  mean_actual_utilities = mean(bt.utilities, 2) # Mean per row => mean for each utility value

end

function support_of_vector(bt::BonesaTuner, pvector)

end

# For now we use random sampling but use lhs here instead, long-term.
function create_random_parameter_vectors(bt::BonesaTuner, num_vectors)
  bt.parameter_vectors = rand(num_params(bt), num_vectors)
end

function num_utility_values(bt::BonesaTuner)
  if bt.num_utility_values == false
    # Evaluate a param vector once to know how many utility values are returned.
    res = evaluate_param_vector_on_problem(bt, bt.problem_set[1], rand(num_params(bt)))
    bt.num_utility_values = length(res) + 1 # Add one for time
  else
    bt.num_utility_values
  end
end

function evaluate_param_vector(bt, pvector)
  parameters = map_to_parameter_values(bt.parameter_set, pvector)
  nu = num_utility_values(bt)
  np = num_problems(bt.problem_set)
  utilities = zeros(np*nu)

  for(i in 1:np)
    p = bt.problem_set[i]
    start = 1 + (i-1) * nu
    tic()
    utilities[ustart:(ustart+nu-2)] = evaluate_param_vector_on_problem(bt, p, parameters)
    utilities[ustart+nu-1] = toq()
  end

  utilities
end

# Find an approximation to the BONESA magic constant c. It is the value
# that makes sum(w) == 50, where w is exp(c * normalized_distance(x, y)) for
# m uniformly, randomly sampled on the l-dimensional hypercube.
function calc_approximate_magic_constant_c(m, l, num_repeats = 5, tolerance = 1e-3)
  # Do a number of random approximations and then average
  mean([find_approximate_c_given_random_sample_of_points(m, l, tolerance) for i in 1:num_repeats])
end

function find_approximate_c_given_random_sample_of_points(m, l, tolerance = 1e-3)
  # Randomly sample m l-dimensional points
  points = rand(l, m)

  # To calc the distances we need the std dev per dimension
  std_devs = std(points, 2)

  # Now calc the distances between points. We only sample them if m is large.
  distances = zeros(Float64, m, m)
  for(i in 1:m)
    x = points[:,i]
    for(j in i:m)
      y = points[:,j]
      distances[i, j] = distances[i, j] = std_distance(x, y, std_devs)
    end
  end

  # Find a value of c by using bisection so that we are close to 50.0
  f = (c) -> average_sum_of_distances_for_c(c, distances) - 50.0
  bisection(f, -100.0, 100.0, tolerance)
end

# Find a value v that makes abs(f(v)) <= tolerance where a and b are guesses
# for v where b > a. This implementation assumes that we can find a negative
# (or positive) value of the function by stepping down (or up) => will not
# work for any function but throws an ArgumentError after many iterations of trying.
function bisection(f, a, b = 2*a, tolerance = 1e-5, max_iterations = 1e4)
  if f(a) > 0 && f(b) > 0
    a = step_value_until_changes_sign(f, a, -1.0, max_iterations, max_iterations)
  elseif f(a) < 0 && f(b) < 0
    b = step_value_until_changes_sign(f, b, 1.0, max_iterations, max_iterations)
  end

  p = (a + b)/2
  error = abs(f(p))
  iterations = 0
  while iterations < max_iterations && error > tolerance
    iterations += 1
    (f(a)*f(p) < 0) ? (b = p) : (a = p)
    p = (a + b)/2
    error = abs(f(p))
  end
  if iterations >= max_iterations
    throw(ArgumentError("Could not find a good approximation after $(iterations) iterations."))
  else
    return p
  end
end

function step_value_until_changes_sign(func, value::Float64, step_size::Float64 = 1.0, max_iterations = 1e3)
  wanted_sign = -1 * sign(func(value))
  iterations = 0
  while iterations < max_iterations && sign(func(value)) != wanted_sign
    iterations += 1
    value = value + step_size
    step_size *= 2
  end
  if iterations >= max_iterations
    throw(ArgumentError("Could not invert the sign after $(iterations) iterations."))
  else
    return value
  end
end

function average_sum_of_distances_for_c(c, distances)
  mean( exp(c * collect(triu(distances))) )
end

function std_distance(x, y, std_devs)
  sqrt(1/l * sum(((x .- y) ./ std_devs).^2))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

include("helpdb/Base.jl")


text:
import Pkg
# import TOML
import UUIDs

@inline function _delete_yanked_versions!(versions)
    versions_deepcopy = deepcopy(versions)
    for (version, info) in pairs(versions_deepcopy)
        if get(info, "yanked", false) == "true"
            delete!(versions, version)
        end
    end
    return nothing
end

"""
    parse_registry(registry::LocalFolderRegistry)

Parse a registry located at the filesystem path `path_to_registry`.

Example usage:
```julia
julia> registry = LocalFolderRegistry("/path/to/my/registry/")
julia> parseresult = parse_registry(registry)
```
"""
@inline function parse_registry(registry::LocalFolderRegistry)
    path_to_registry = registry.path
    base_path = abspath(path_to_registry)
    registry = Pkg.TOML.parsefile(joinpath(base_path, "Registry.toml"))
    packages = registry["packages"]
    parseresult = ParsedRegistry()
    for (uuid_string, info) in packages
        uuid = UUIDs.UUID(uuid_string)
        name = info["name"]
        path = joinpath(base_path, info["path"])
        package = FullUnversionedPackage(name, uuid)
        versions = Pkg.TOML.parsefile(joinpath(path, "Versions.toml"))
        _delete_yanked_versions!(versions)
        if !isempty(versions)
            all_version_strings = collect(keys(versions))
            all_version_numbers = VersionNumber.(all_version_strings)
            i = argmax(all_version_numbers)
            max_version_string = all_version_strings[i]
            max_version_tree = versions[max_version_string]["git-tree-sha1"]
            parseresult[package] = VersionAndTree(max_version_string, max_version_tree)
        end
    end
    return parseresult
end

@inline function parse_registries!(cache, registry_list::AbstractVector{<:AbstractRegistry})
    for registry in registry_list
        if !haskey(cache, registry)
            parseresult = parse_registry(registry)
            cache[registry] = parseresult
        end
    end
    return nothing
end


text:
struct OrderStatistic{n,r,B} <: AbstractMeasure
    base::B
end

function basemeasure(d::OrderStatistic{n,r})
    WeightedMeasure(logbeta(n+1, r+1), d.base)
end

function logdensity(d::OrderStatistic{n,r}, x)
    logF(x) = logcdf(d.base, x)
    logFᶜ(x) = logccdf(d.base, x)
    return (r - 1) * logF(x) + (n - r) * logFᶜ(x)
end


text:
## Test Laplace transform pairs

"""
    TransformPair{T,V}

A pair of functions for analyzing an inverse Laplace transform method.
Field `ft` is the real-space function. Field `fs` is the Laplace-space function.
"""
struct TransformPair{T,V}
    ft::T
    fs::V
end

"""
    ILtPair(ilt::AbstractILt, ft::Function)

return an object of type `ILtPair` that associates `ilt` the inverse Laplace transform of
a function with it "exact" numerical inverse `ft`. Calling `abserr(p, t)` returns the
absolute error between the inverse transform and the exact value.

# Example

This example compares the inversion using the Weeks algorithm of the Laplace transform
of `cos(t)` to its exact value at `t=1.0`.
```
julia> p = ILtPair(Weeks(s -> s/(1+s^2)), cos);
julia> abserr(p, 1.0)

0.0
```
"""
struct ILtPair{T,V}
    ilt::T
    ft::V
end

# type ILtPair{T,V} <: AbstractILt
"""
    iltpair_power(n)

Return a `TransformPair` for the power function `x^n`. This can be used
to construct an `ILTPair`.

```julia-repl
julia> p  = Talbot(iltpair_power(3));

julia> Float64(abserr(p, 1)) # test Talbot method for Laplace transform of x^3, evaluated at 1
2.0820366247539812e-26
```
"""
function iltpair_power(n)
    fs = s -> s^(-n - 1) * SpecialFunctions.gamma(1 + n)
    ft = t -> t^n
    return TransformPair(ft, fs)
end

"""
    abserr(p::ILtPair, t)

Compute the absolute error between the estimated inverse Laplace transform and
"exact" numerical solution contained in `p` at the point `t`. See `ILtPair`.
"""
abserr(p::ILtPair, t) = abs(p.ilt(t) - p.ft(t))

## FIXME: Consider forwarding from MacroTools
# Create ILtPair with ilt type
for ilttype in (:Talbot, :GWR, :Weeks)
    @eval $(ilttype)(p::TransformPair, args...) = ILtPair($(ilttype)(p.fs, args...), p.ft)
end

for f in (:optimize,)
    @eval $(f)(p::ILtPair, args...) = $(f)(p.ilt, args...)
end


text:
#systems parameterized in the standard "generalized" form
"""Generalized Lotka-Volterra ecosystem with intrisic growth rate r and community matrix A."""
type GeneralizedLotkaVolterra{S<:AbstractVector{Float64},T<:AbstractMatrix{Float64}} <: AbstractLotkaVolterra
  r::S
  A::T
  function GeneralizedLotkaVolterra{S,T}(r::S,A::T)
    size(A,1) == size(A,2) == length(r) ? new(r,A) : error("dimension mismatch")
  end
end
function GeneralizedLotkaVolterra(r::AbstractVector{Float64},A::AbstractMatrix{Float64})
    GeneralizedLotkaVolterra{typeof(r),typeof(A)}(r,A)
end
function GeneralizedLotkaVolterra(A::AbstractMatrix{Float64})
    N = size(A,1)
    GeneralizedLotkaVolterra(ones(N),A)
end

"""The number of species in the ecosystem"""
nspecies(Z::GeneralizedLotkaVolterra) = length(Z.r)
"""The ecosystem community matrix"""
communitymatrix(Z::GeneralizedLotkaVolterra) = Z.A
"""The intrinsic gorwth rate of the ecosystem"""
intrinsicrate(Z::GeneralizedLotkaVolterra) = Z.r


text:
module mcrosssection1
using FinEtools
using FinEtoolsFlexBeams.CrossSectionModule
using Test
function test()
    cs = CrossSectionModule.CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
    par = cs.parameters(0.0)
    for (c, r) in zip((par.A, par.J, par.I1, par.I2, par.I3), (112.75829799164978, 2023.5649824692157, 2023.5649824692157, 1011.7824912346078, 1011.7824912346078))
        @test c ≈ r
    end
    true
end
end
using .mcrosssection1
mcrosssection1.test()

module mcrosssection2
using FinEtools
using FinEtoolsFlexBeams.CrossSectionModule
using Test
function test()
    cs = CrossSectionModule.CrossSectionRectangle(s -> 42.0, s -> 4.2, s -> [0.0, 0.0, 1.0])
    par = cs.parameters(0.0)
    for (c, r) in zip((par.A, par.J, par.I1, par.I2, par.I3), (176.4, 970.849152, 26190.108000000004, 259.30800000000005, 25930.800000000003))
        @test c ≈ r
    end
    true
end
end
using .mcrosssection2
mcrosssection2.test()

module mcrosssection3
using FinEtools
using FinEtoolsFlexBeams.CrossSectionModule
using Test
function test()
    cs = CrossSectionModule.CrossSectionRectangle(s -> 1.3 * 4.2, s -> 4.2, s -> [0.0, 0.0, 1.0])
    par = cs.parameters(0.0)
    for (c, r) in zip((par.A, par.J, par.I1, par.I2, par.I3), (22.932000000000006, 71.66370760423138, 90.68000760000004, 33.71004000000001, 56.969967600000025))
         @test c ≈ r
    end
 true
end
end
using .mcrosssection3
mcrosssection3.test()

module mcrosssection4
using FinEtools
using FinEtoolsFlexBeams.CrossSectionModule
using Test
function test()
    R = 0.5
    cs = CrossSectionModule.CrossSectionCircle(s -> (1/2+2*s)*R, s -> [0.0, 0.0, 1.0])
    par = cs.parameters(0.0)
    for (c, r) in zip((par.A, par.J, par.I1, par.I2, par.I3), (0.19634954084936207, 0.006135923151542565, 0.006135923151542565, 0.0030679615757712823, 0.0030679615757712823))
        @test c ≈ r
    end
    true
end
end
using .mcrosssection4
mcrosssection4.test()


text:
begin
    rng = Random.MersenneTwister(1)
    n = 1000
    p = 0.02
    q = 1 - sqrt(1-p)
    Areal = sprandn(rng, n, n, p)
    Breal = randn(rng, n)
    Bcomplex = Breal + randn(rng, n) * im
    Acomplex = sprandn(rng, n, n, q) + sprandn(rng, n, n, q) * im
    @testset "symmetric/Hermitian sparse multiply with $S($U)" for S in (Symmetric, Hermitian), U in (:U, :L), (A, B) in ((Areal,Breal), (Acomplex,Bcomplex))
        Asym = S(A, U)
        As = sparse(Asym) # takes most time
        @test which(mul!, (typeof(B), typeof(Asym), typeof(B))).module in
            (SparseArrays, SparseWrappers)
        @test norm(Asym * B - As * B, Inf) <= eps() * n * p * 2
    end
end


text:
using Test
using Base.BinaryPlatforms
import Libdl
using BinaryBuilderBase
using BinaryBuilderBase: template, dlopen_flags_str

# The platform we're running on
const platform = HostPlatform()

@testset "Products" begin
    @test template(raw"$libdir/foo-$arch/$nbits/bar-$target", Platform("x86_64", "windows")) ==
        "bin/foo-x86_64/64/bar-x86_64-w64-mingw32"
    @test template(raw"$target/$nbits/$arch/$libdir", Platform("x86_64", "linux"; libc = "musl")) ==
        "x86_64-linux-musl/64/x86_64/lib"

    lp = LibraryProduct("libfakechroot", :libfakechroot, "lib/fakechroot")
    @test lp.libnames ==  ["libfakechroot"]
    @test lp.dir_paths == ["lib/fakechroot"]
    ep = ExecutableProduct("fooify", :fooify, "bin/foo_inc")
    @test ep.binnames ==  ["fooify"]

    @test_throws ErrorException LibraryProduct("sin", :sin)
    @test_throws ErrorException ExecutableProduct("convert", :convert)
    @test_throws ErrorException FileProduct("open", :open)

    # Test sorting of products....
    @test sort([LibraryProduct("libbar", :libbar), ExecutableProduct("foo", :foo), FrameworkProduct("buzz", :buzz)]) ==
        [FrameworkProduct("buzz", :buzz), ExecutableProduct("foo", :foo), LibraryProduct("libbar", :libbar)]
    # ...and products info
    p1 = LibraryProduct(["libchafa"], :libchafa, ) => Dict("soname" => "libchafa.so.0","path" => "lib/libchafa.so")
    p2 = ExecutableProduct(["chafa"], :chafa, ) => Dict("path" => "bin/chafa")
    products_info = Dict{Product,Any}(p1, p2)
    @test sort(products_info) == [p2, p1]

    temp_prefix() do prefix
        # Test that basic satisfication is not guaranteed
        e_path = joinpath(bindir(prefix), "fooifier")
        l_path = joinpath(last(libdirs(prefix)), "libfoo.$(Libdl.dlext)")
        e = ExecutableProduct("fooifier", :fooifier)
        ef = FileProduct(joinpath("bin", "fooifier"), :fooifier)
        l = LibraryProduct("libfoo", :libfoo)
        lf = FileProduct(l_path, :libfoo)

        @test @test_logs (:info, r"does not exist") !satisfied(e, prefix; verbose=true)
        @test @test_logs (:info, r"not found") !satisfied(ef, prefix; verbose=true)
        @test @test_logs (:info, r"^Could not locate") !satisfied(l, prefix; verbose=true)
        @test @test_logs (:info, r"^Could not locate") !satisfied(l, prefix; verbose=true, isolate=true)
        @test @test_logs (:info, r"^FileProduct .* not found") !satisfied(lf, prefix; verbose=true)

        # Test that simply creating a file that is not executable doesn't
        # satisfy an Executable Product (and say it's on Linux so it doesn't
        # complain about the lack of an .exe extension)
        mkpath(bindir(prefix))
        touch(e_path)
        @test @test_logs (:info, r"^FileProduct .* found at") satisfied(ef, prefix; verbose=true)
        @static if !Sys.iswindows()
            # Windows doesn't care about executable bit, grumble grumble
            @test_logs (:info, r"is not executable") (:info, r"does not exist") begin
                @test !satisfied(e, prefix; verbose=true, platform=Platform("x86_64", "linux"))
            end
        end

        # Make it executable and ensure this does satisfy the Executable
        chmod(e_path, 0o777)
        @test_logs (:info, r"matches our search criteria") begin
            @test satisfied(e, prefix; verbose=true, platform=Platform("x86_64", "linux"))
    end

        # Remove it and add a `$(path).exe` version to check again, this
        # time saying it's a Windows executable
        Base.rm(e_path; force=true)
        touch("$(e_path).exe")
        chmod("$(e_path).exe", 0o777)
        @test locate(e, prefix; platform=Platform("x86_64", "windows")) == "$(e_path).exe"

        # Test that simply creating a library file doesn't satisfy it if we are
        # testing something that matches the current platform's dynamic library
        # naming scheme, because it must be `dlopen()`able.
        mkpath(last(libdirs(prefix)))
        touch(l_path)
        @test @test_logs (:info, r"^FileProduct .* found at") satisfied(lf, prefix; verbose=true)
        @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") (:info, r"cannot be dlopen'ed") (:info, r"^Could not locate") !satisfied(l, prefix; verbose=true)
        @test @test_logs (:info, r"^FileProduct .* found at") satisfied(lf, prefix; verbose=true, isolate=true)
        @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") (:info, r"cannot be dlopen'ed") (:info, r"^Could not locate") !satisfied(l, prefix; verbose=true, isolate=true)

        # But if it is from a different platform, simple existence will be
        # enough to satisfy a LibraryProduct
        @static if Sys.iswindows()
            p = Platform("x86_64", "linux")
            mkpath(last(libdirs(prefix, p)))
            l_path = joinpath(last(libdirs(prefix, p)), "libfoo.so")
            touch(l_path)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(l, prefix; verbose=true, platform=p)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(l, prefix; verbose=true, platform=p, isolate=true)

            # Check LibraryProduct objects with explicit directory paths
            ld = LibraryProduct("libfoo", :libfoo)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(ld, prefix; verbose=true, platform=p)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(ld, prefix; verbose=true, platform=p, isolate=true)
        else
            p = Platform("x86_64", "windows")
            mkpath(last(libdirs(prefix, p)))
            l_path = joinpath(last(libdirs(prefix, p)), "libfoo.dll")
            touch(l_path)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(l, prefix; verbose=true, platform=p)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(l, prefix; verbose=true, platform=p, isolate=true)

            # Check LibraryProduct objects with explicit directory paths
            ld = LibraryProduct("libfoo", :libfoo)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(ld, prefix; verbose=true, platform=p)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(ld, prefix; verbose=true, platform=p, isolate=true)
        end
    end

    # Ensure that the test suite thinks that these libraries are foreign
    # so that it doesn't try to `dlopen()` them:
    foreign_platform = if platform == Platform("aarch64", "linux")
        # Arbitrary architecture that is not dlopen()'able
        Platform("powerpc64le", "linux")
    else
        # If we're not Platform("aarch64", "linux"), then say the libraries are
        Platform("aarch64", "linux")
    end

    # Test for valid library name permutations
    for ext in ["so", "so.1", "so.1.2", "so.1.2.3"]
        temp_prefix() do prefix
            l_path = joinpath(last(libdirs(prefix, foreign_platform)), "libfoo.$ext")
            l = LibraryProduct("libfoo", :libfoo)
            mkdir(dirname(l_path))
            touch(l_path)
            @test @test_logs (:info, r"^Found a valid") (:info, r"matches our search criteria") satisfied(l, prefix; verbose=true, platform=foreign_platform)
        end
    end

    # Test for invalid library name permutations
    for ext in ["1.so", "so.1.2.3a", "so.1.a"]
        temp_prefix() do prefix
            l_path = joinpath(last(libdirs(prefix, foreign_platform)), "libfoo.$ext")
            l = LibraryProduct("libfoo", :libfoo)
            mkdir(dirname(l_path))
            touch(l_path)
            if ext == "1.so"
                @test_logs (:info, r"^Found a valid") (:info, r"^Could not locate") begin
                    @test !satisfied(l, prefix; verbose=true, platform=foreign_platform)
                end
            else
                @test_logs (:info, r"^Could not locate") begin
                    @test !satisfied(l, prefix; verbose=true, platform=foreign_platform)
                end
            end
        end
    end

    # Test for proper repr behavior
    temp_prefix() do prefix
        l = LibraryProduct("libfoo", :libfoo)
        @test repr(l) == "LibraryProduct($(repr(["libfoo"])), :libfoo)"
        l = LibraryProduct(["libfoo", "libfoo2"], :libfoo)
        @test repr(l) == "LibraryProduct($(repr(["libfoo", "libfoo2"])), :libfoo)"

        e = ExecutableProduct("fooifier", :fooifier)
        @test repr(e) == "ExecutableProduct([\"fooifier\"], :fooifier)"
        e = ExecutableProduct("fooifier", :fooifier, "bin/qux")
        @test repr(e) == "ExecutableProduct([\"fooifier\"], :fooifier, \"bin/qux\")"

        f = FileProduct(joinpath("etc", "fooifier"), :foo_conf)
        @test repr(f) == "FileProduct([$(repr(joinpath("etc", "fooifier")))], :foo_conf)"

        f = FileProduct(joinpath(prefix, "etc", "foo.conf"), :foo_conf)
        @test repr(f) == "FileProduct([$(repr(joinpath(prefix, "etc", "foo.conf")))], :foo_conf)"
    end

    # Test that FileProduct's can have `${target}` within their paths:
    temp_prefix() do prefix
        multilib_dir = joinpath(prefix, "foo", triplet(platform))
        mkpath(multilib_dir)
        touch(joinpath(multilib_dir, "bar"))

        for path in ("foo/\$target/bar", "foo/\${target}/bar")
            f = FileProduct(path, :bar)
            @test  @test_logs (:info, r"^FileProduct .* found at") satisfied(f, prefix; verbose=true, platform=platform)
        end
    end
end

@testset "Dlopen flags" begin
    lp = LibraryProduct("libfoo2", :libfoo2; dlopen_flags=[:RTLD_GLOBAL, :RTLD_NOLOAD])
    @test lp.dlopen_flags == [:RTLD_GLOBAL, :RTLD_NOLOAD]
    fp = FrameworkProduct("libfoo2", :libfoo2; dlopen_flags=[:RTLD_GLOBAL, :RTLD_NOLOAD])
    @test fp.libraryproduct.dlopen_flags == [:RTLD_GLOBAL, :RTLD_NOLOAD]
    for p in (lp, fp)
        flag_str = dlopen_flags_str(p)
        @test flag_str == "RTLD_GLOBAL | RTLD_NOLOAD"
        @test Libdl.eval(Meta.parse(flag_str)) == (Libdl.RTLD_NOLOAD | Libdl.RTLD_GLOBAL)
    end
    lp = LibraryProduct("libfoo2", :libfoo2; dont_dlopen=true)
    @test dlopen_flags_str(lp) == "nothing"
end