text:
function FieldCircularOrbit(arg0::FieldOrbit)
    return FieldCircularOrbit((FieldOrbit,), arg0)
end

function FieldCircularOrbit(arg0::FieldPVCoordinates, arg1::Frame, arg2::FieldAbsoluteDate, arg3::RealFieldElement)
    return FieldCircularOrbit((FieldPVCoordinates, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3)
end

function FieldCircularOrbit(arg0::RealFieldElement, arg1::RealFieldElement, arg10::RealFieldElement, arg11::RealFieldElement, arg12::PositionAngle, arg13::Frame, arg14::FieldAbsoluteDate, arg15::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::RealFieldElement, arg7::RealFieldElement, arg8::RealFieldElement, arg9::RealFieldElement)
    return FieldCircularOrbit((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, PositionAngle, Frame, FieldAbsoluteDate, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldCircularOrbit(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::PositionAngle, arg7::Frame, arg8::FieldAbsoluteDate, arg9::RealFieldElement)
    return FieldCircularOrbit((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, PositionAngle, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldCircularOrbit(arg0::TimeStampedFieldPVCoordinates, arg1::Frame, arg2::RealFieldElement)
    return FieldCircularOrbit((TimeStampedFieldPVCoordinates, Frame, RealFieldElement), arg0, arg1, arg2)
end

function add_kepler_contribution(obj::FieldCircularOrbit, arg0::PositionAngle, arg1::RealFieldElement, arg2::Vector{RealFieldElement})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, RealFieldElement, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function eccentric_to_mean(::Type{FieldCircularOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldCircularOrbit, "eccentricToMean", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function eccentric_to_true(::Type{FieldCircularOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldCircularOrbit, "eccentricToTrue", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function get_a(obj::FieldCircularOrbit)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_a_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getADot", RealFieldElement, ())
end

function get_alpha(obj::FieldCircularOrbit, arg0::PositionAngle)
    return jcall(obj, "getAlpha", RealFieldElement, (PositionAngle,), arg0)
end

function get_alpha_dot(obj::FieldCircularOrbit, arg0::PositionAngle)
    return jcall(obj, "getAlphaDot", RealFieldElement, (PositionAngle,), arg0)
end

function get_alpha_e(obj::FieldCircularOrbit)
    return jcall(obj, "getAlphaE", RealFieldElement, ())
end

function get_alpha_e_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getAlphaEDot", RealFieldElement, ())
end

function get_alpha_m(obj::FieldCircularOrbit)
    return jcall(obj, "getAlphaM", RealFieldElement, ())
end

function get_alpha_m_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getAlphaMDot", RealFieldElement, ())
end

function get_alpha_v(obj::FieldCircularOrbit)
    return jcall(obj, "getAlphaV", RealFieldElement, ())
end

function get_alpha_v_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getAlphaVDot", RealFieldElement, ())
end

function get_circular_ex(obj::FieldCircularOrbit)
    return jcall(obj, "getCircularEx", RealFieldElement, ())
end

function get_circular_ex_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getCircularExDot", RealFieldElement, ())
end

function get_circular_ey(obj::FieldCircularOrbit)
    return jcall(obj, "getCircularEy", RealFieldElement, ())
end

function get_circular_ey_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getCircularEyDot", RealFieldElement, ())
end

function get_e(obj::FieldCircularOrbit)
    return jcall(obj, "getE", RealFieldElement, ())
end

function get_e_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getEDot", RealFieldElement, ())
end

function get_equinoctial_ex(obj::FieldCircularOrbit)
    return jcall(obj, "getEquinoctialEx", RealFieldElement, ())
end

function get_equinoctial_ex_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getEquinoctialExDot", RealFieldElement, ())
end

function get_equinoctial_ey(obj::FieldCircularOrbit)
    return jcall(obj, "getEquinoctialEy", RealFieldElement, ())
end

function get_equinoctial_ey_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getEquinoctialEyDot", RealFieldElement, ())
end

function get_hx(obj::FieldCircularOrbit)
    return jcall(obj, "getHx", RealFieldElement, ())
end

function get_hx_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getHxDot", RealFieldElement, ())
end

function get_hy(obj::FieldCircularOrbit)
    return jcall(obj, "getHy", RealFieldElement, ())
end

function get_hy_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getHyDot", RealFieldElement, ())
end

function get_i(obj::FieldCircularOrbit)
    return jcall(obj, "getI", RealFieldElement, ())
end

function get_i_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getIDot", RealFieldElement, ())
end

function get_le(obj::FieldCircularOrbit)
    return jcall(obj, "getLE", RealFieldElement, ())
end

function get_le_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getLEDot", RealFieldElement, ())
end

function get_lm(obj::FieldCircularOrbit)
    return jcall(obj, "getLM", RealFieldElement, ())
end

function get_lm_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getLMDot", RealFieldElement, ())
end

function get_lv(obj::FieldCircularOrbit)
    return jcall(obj, "getLv", RealFieldElement, ())
end

function get_lv_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getLvDot", RealFieldElement, ())
end

function get_right_ascension_of_ascending_node(obj::FieldCircularOrbit)
    return jcall(obj, "getRightAscensionOfAscendingNode", RealFieldElement, ())
end

function get_right_ascension_of_ascending_node_dot(obj::FieldCircularOrbit)
    return jcall(obj, "getRightAscensionOfAscendingNodeDot", RealFieldElement, ())
end

function get_type(obj::FieldCircularOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function has_derivatives(obj::FieldCircularOrbit)
    return jcall(obj, "hasDerivatives", jboolean, ())
end

function interpolate(obj::FieldCircularOrbit, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldCircularOrbit, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function mean_to_eccentric(::Type{FieldCircularOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldCircularOrbit, "meanToEccentric", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function shifted_by(obj::FieldCircularOrbit, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldCircularOrbit, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldCircularOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldCircularOrbit, (jdouble,), arg0)
end

function to_orbit(obj::FieldCircularOrbit)
    return jcall(obj, "toOrbit", CircularOrbit, ())
end

function to_string(obj::FieldCircularOrbit)
    return jcall(obj, "toString", JString, ())
end

function true_to_eccentric(::Type{FieldCircularOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldCircularOrbit, "trueToEccentric", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end


text:
# plots a reach set

# The next line is needed for a world age issue in Plots:
# https://github.com/JuliaPlots/Plots.jl/issues/1047
# This line must also come before 'using Plots' to prevent an error:
# ERROR: MethodError: no method matching convert(::Type{AssertionError}, ::String)
import GR

using Reachability, Plots

import Reachability.@timing

function plot_reach(result::AbstractSolution, name::String="plot")
    if result.options[:mode] != "reach"
        # ignore plotting command
        return
    end
    info("Plotting...")
    # GR plotting backend does not support subindices like x₁ in savefig
    is_gr = Plots.backend() == Plots.GRBackend()
    @timing begin
        plot(result, use_subindices=!is_gr)
        savefig("$name.png")
    end
end

ext: jl

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
function array_string_view(data)
    isa(data, Missing) && return data

    if length(data) == 1
        return @sprintf "[%f]" data[begin]
    elseif length(data) == 2
        return @sprintf "[%f, ..., %f]" data[begin] data[end]
    elseif length(data) == 3
        return @sprintf "[%f, ..., %f, %f]" data[begin] data[end-1] data[end]
    elseif length(data) ≥ 4
        return @sprintf "[%f, ..., %f, %f, %f]" data[begin] data[end-2] data[end-1] data[end]
    end

    return data # data is empty
end

function header_view(is_suptype::Bool, name)
    is_suptype && begin
        return @sprintf "%4s%s:\n\n" "" name
    end
    return @sprintf "%s:\n--------------------------------------\n" name
end

# TODO: Add a print_field method, that handles offset

# TODO: Add a pretty_print method to round floats

function Base.show(io::IO, m::Model)
    println(io, "  Model: $(name(m))")
    formula_here = (isa(formula(m), Missing) ? formula(m) : "loaded")
    println(io, "  -------------------")
    println(io, "    objective:         $(objective(m))")
    println(io, "    gradient:          $(gradient(m))")
    println(io, "    initial iterate:   $(array_string_view(initial_iterate(m)))")
    println(io, "    dimension:         $(dimension(m))")
    println(io, "    directory:         $(directory(m))")
    println(io, "    objective formula: $(formula_here)")
    flush(io)
    return nothing
end

function Base.show(io::IO, o::DriverOptions)
    println(io, "    Options:")
    println(io, "      samples:        $(samples(o))")
    println(io, "      Δ_max:          $(Δ_max(o))")
    println(io, "      δ_tol:          $(δ_tol(o))")
    println(io, "      ϵ_tol:          $(ϵ_tol(o))")
    println(io, "      max_iterations: $(max_iterations(o))")
    flush(io)
    return nothing
end

function Base.show(io::IO, d::Driver)
    println(io, "  Driver:")
    println(io, "  -------------------")
    println(io, "    S_update:  $(S_update(d))")
    println(io, "    QN_update: $(QN_update(d))")
    println(io, "    pflag:     $(pflag(d))")
    show(io, options(d))
    flush(io)
    return nothing
end

function Base.show(io::IO, w::Weaver)
    println(io, "    Weaver:")
    println(io, "      f_vals:   $(array_string_view(f_vals(w)))")
    println(io, "      ∇f_norms: $(array_string_view(∇f_norms(w)))")
    println(io, "      Δ_vals:   $(array_string_view(Δ_vals(w)))")
    println(io, "      p_norms:  $(array_string_view(p_norms(w)))")
    println(io, "      ρ_vals:   $(array_string_view(ρ_vals(w)))")
    flush(io)
    return nothing
end

function Base.show(io::IO, p::BlockOptProfile)
    println(io, "    Profile:")
    println(io, "      trs_counter: $(evaluations(trs_counter(p)))")
    println(io, "      trs_timer:   $(Δt(trs_timer(p)))")
    println(io, "      ghs_counter: $(evaluations(ghs_counter(p)))")
    println(io, "      ghs_timer:   $(Δt(ghs_timer(p)))")
    flush(io)
    return nothing
end

function Base.show(io::IO, t::BlockOptTrace)
    println(io, "  Trace:")
    println(io, "  -------------------")
    show(io, weaver(t))
    show(io, profile(t))
    flush(io)
    return nothing
end

function Base.show(io::IO, b::BlockOptBackend)
    show(io, model(b))
    show(io, driver(b))
    flush(io)
    return nothing
end

function Base.show(io::IO, s::Simulation)
    if evaluations(trs_counter(s)) == 0
        # initial simulation display
        show(io, model(trace(s)))
        show(io, driver(trace(s)))
    else
        # display after simulation complete
        pass = (∇fₖ_norm(s) ≤ ϵ_tol(backend(s)) ? true : false)
        print(io, "$(pass ? "SUCCESS" : "FAIL")")
        print(io, " $(∇fₖ_norm(s)) $(pass ? "≤" : "≰") $(ϵ_tol(backend(s)))")
        println(io, " in $(evaluations(trs_counter(s))) steps")
        println(io, "--------------------------------------")
        println(io, "  Minimum f:      $(minimum(f_vals(s)))")
        println(io, "  Minimum ||∇f||: $(minimum(∇f_norms(s)))")
        println(io, "  Minimum Δ:      $(minimum(Δ_vals(s)))")
        println(io, "  Minimum Step:   $(minimum(p_norms(s)))")
        println(io, "")
        show(io, backend(s))
        show(io, trace(s))
    end
    println(io, "")
    flush(io)
    return nothing
end

# function show_driver(io::IO, d::Driver)
#     println(io,"\n%s %s Driver: n = %d, w = %d\n" d.name d.QN d.n d.w
#     println(io,"-----------------------------------------------------------\n"
#     println(io,"  Hessian Samples: %d\n" 2*d.w
#     println(io,"  Terminal Conditions: ||∇fₖ|| < %1.2e or %d iterations\n\n" d.ϵ d.max_iter
#     return nothing
# end


text:
"""
```
type optimization_result{T}
    minimizer::Vector{T}
    minimum::T
    converged::Bool
    iterations::Int

```
Container type for various optimization outputs
"""
type optimization_result{T}
    minimizer::Vector{T}
    minimum::T
    converged::Bool
    iterations::Int
end

"""
```
optimize!(m::AbstractModel, data::Matrix;
          method::Symbol       = :csminwel,
          xtol::Real           = 1e-32,  # default from Optim.jl
          ftol::Float64        = 1e-14,  # Default from csminwel
          grtol::Real          = 1e-8,   # default from Optim.jl
          iterations::Int      = 1000,
          store_trace::Bool    = false,
          show_trace::Bool     = false,
          extended_trace::Bool = false,
          verbose::Symbol      = :none)
```

Wrapper function to send a model to csminwel (or another optimization routine).
"""
function optimize!(m::AbstractModel,
                   data::Matrix;
                   method::Symbol       = :csminwel,
                   xtol::Real           = 1e-32,  # default from Optim.jl
                   ftol::Float64        = 1e-14,  # Default from csminwel
                   grtol::Real          = 1e-8,   # default from Optim.jl
                   iterations::Int      = 1000,
                   store_trace::Bool    = false,
                   show_trace::Bool     = false,
                   extended_trace::Bool = false,
                   mle::Bool            = false, # default from estimate.jl
                   step_size::Float64   = .01,
                   verbose::Symbol      = :none)

    ########################################################################################
    ### Step 1: Setup
    ########################################################################################

    # For now, only csminwel should be used
    optimizer = if method == :csminwel
        csminwel
    elseif method == :simulated_annealing
        simulated_annealing
    elseif method == :nelder_mead
        nelder_mead
    elseif method == :combined_optimizer
        combined_optimizer
    elseif method == :lbfgs
        lbfgs
    else
        error("Method ", method, " is not supported.")
    end

    # Inputs to optimization
    H0             = 1e-4 * eye(n_parameters_free(m))
    para_free_inds = find([!θ.fixed for θ in m.parameters])
    x_model        = transform_to_real_line(m.parameters)
    x_opt          = x_model[para_free_inds]

    ########################################################################################
    ### Step 2: Initialize f_opt
    ########################################################################################

    function f_opt(x_opt)
        try
            x_model[para_free_inds] = x_opt
            transform_to_model_space!(m,x_model)
        catch
            return Inf
        end
        if mle
            out = -likelihood(m, data; catch_errors = true)
        else
            out = -posterior(m, data; catch_errors=true)
        end
        out = !isnan(out) ? out : Inf
        return out
    end

    ########################################################################################
    ### Step 3: Optimizer-specific setup, call optimizer
    ########################################################################################

    # variables used across several optimizers
    rng = m.rng
    temperature = get_setting(m, :simulated_annealing_temperature)
    max_cycles = get_setting(m, :combined_optimizer_max_cycles)
    block_frac = get_setting(m, :simulated_annealing_block_proportion)
    H_ = nothing

    function neighbor_dsge!(x, x_proposal)
        # This function computes a proposal "next step" during simulated annealing.
        # Inputs:
        # - `x`: current position (of non-fixed states)
        # - `x_proposal`: proposed next position (of non-fixed states).
        #                 (passed in for pre-allocation purposes)
        # Outputs:
        # - `x_proposal`

        T = eltype(x)
        npara = length(x)
        subset_inds = []
        while length(subset_inds) == 0
            subset_inds = randsubseq(para_free_inds,block_frac)
        end

        # Convert x_proposal to model space and expand to full parameter vector
        x_all = T[p.value for p in m.parameters]  # to get fixed values
        x_all[para_free_inds] = x                 # this is from real line

        x_all_model = transform_to_model_space(m.parameters, x_all)
        x_proposal_all = copy(x_all_model)

        success = false
        while !success
            # take a step in model space
            for i in subset_inds
                prior_var = moments(m.parameters[i])[2]#moments(get(m.parameters[i].prior))[2]
                proposal_in_bounds = false
                proposal = x_all_model[i]
                lower = m.parameters[i].valuebounds[1]
                upper = m.parameters[i].valuebounds[2]
                # draw a new parameter value, and redraw if out of bounds
                while !proposal_in_bounds
                    r = rand([-1 1]) * rand()
                    proposal = x_all_model[i] + (r * step_size * prior_var)
                    if lower < proposal < upper
                        proposal_in_bounds = true
                    end
                end
                @inbounds x_proposal_all[i] = proposal
            end

            # check that model can be solved
            try
                update!(m, x_proposal_all)
                solve(m)
                x_proposal_all = transform_to_real_line(m.parameters, x_proposal_all)
                success = true
            catch ex
                if !(typeof(ex) in [DomainError, ParamBoundsError, GensysError])
                    rethrow(ex)
                end
            end
        end

        x_proposal[1:end] = x_proposal_all[para_free_inds]

        return
    end

    rng = m.rng
    temperature = get_setting(m, :simulated_annealing_temperature)
    if method == :simulated_annealing
       opt_result = optimizer(f_opt, x_opt;
                        iterations = iterations, step_size = step_size,
                        store_trace = store_trace, show_trace = show_trace, extended_trace = extended_trace,
                        neighbor! = neighbor_dsge!, verbose = verbose, rng = rng, temperature = temperature)
       converged = opt_result.iteration_converged
       out = optimization_result(opt_result.minimizer, opt_result.minimum, converged, opt_result.iterations)

    elseif method == :nelder_mead
        opt_result = optimizer(f_opt, x_opt;
                               iterations = iterations,
                               store_trace = store_trace, show_trace = show_trace,
                               extended_trace = extended_trace, verbose = verbose, rng = rng)
       converged = opt_result.iteration_converged
       out = optimization_result(opt_result.minimizer, opt_result.minimum, converged, opt_result.iterations)

    elseif method == :csminwel
        opt_result, H_ = optimizer(f_opt, x_opt, H0;
                        xtol = xtol, ftol = ftol, grtol = grtol, iterations = iterations,
                        store_trace = store_trace, show_trace = show_trace, extended_trace = extended_trace,
                        verbose = verbose, rng = rng)
        converged = opt_result.g_converged || opt_result.f_converged #|| opt_result.x_converged
        out = optimization_result(opt_result.minimizer, opt_result.minimum, converged, opt_result.iterations)

    elseif method == :lbfgs
        opt_result = optimizer(f_opt, x_opt;
                        xtol = xtol, ftol = ftol, grtol = grtol, iterations = iterations,
                        store_trace = store_trace, show_trace = show_trace, extended_trace = extended_trace,
                        verbose = verbose, rng = rng)
        converged = opt_result.g_converged || opt_result.f_converged #|| opt_result.x_converged
        out = optimization_result(opt_result.minimizer, opt_result.minimum, converged, opt_result.iterations)

    elseif method == :combined_optimizer
        opt_result = optimizer(f_opt, x_opt;
                        xtol = xtol, ftol = ftol, grtol = grtol, iterations = iterations, step_size = step_size,
                        store_trace = store_trace, show_trace = show_trace, extended_trace = extended_trace,
                        neighbor! = neighbor_dsge!, verbose = verbose, rng = rng, temperature = temperature,
                               max_cycles = max_cycles)
        converged = opt_result.g_converged || opt_result.f_converged || opt_result.x_converged
        converged = opt_result.method == "Simulated Annealing" ? opt_result.iteration_converged : converged
        out = optimization_result(opt_result.minimizer, opt_result.minimum, converged, opt_result.iterations)
    end

    ########################################################################################
    ### Step 4: transform output, populate Hessian
    ########################################################################################

    x_model[para_free_inds] = out.minimizer
    transform_to_model_space!(m, x_model)

    # Match original dimensions
    out.minimizer = map(θ -> θ.value, m.parameters)

    H = zeros(n_parameters(m), n_parameters(m))
    if H_ != nothing

        # Fill in rows/cols of zeros corresponding to location of fixed parameters
        # For each row corresponding to a free parameter, fill in columns corresponding to
        # free parameters. Everything else is 0.
        for (row_free, row_full) in enumerate(para_free_inds)
            H[row_full,para_free_inds] = H_[row_free,:]
        end

    end

    return out, H
end


text:
# Autogenerated wrapper script for Libuuid_jll for powerpc64le-linux-gnu
export libuuid

JLLWrappers.@generate_wrapper_header("Libuuid")
JLLWrappers.@declare_library_product(libuuid, "libuuid.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libuuid,
        "lib/libuuid.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Langmuir turbulence example
#
# This example implements the Langmuir turbulence simulation reported in section
# 4 of
#
# [McWilliams, J. C. et al., "Langmuir Turbulence in the ocean," Journal of Fluid Mechanics (1997)](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/langmuir-turbulence-in-the-ocean/638FD0E368140E5972144348DB930A38).
#
# This example demonstrates:
#
#   * how to run large eddy simulations with surface wave effects
#     via the Craik-Leibovich approximation

using Oceananigans

# ## Model set-up
#
# To build the model, we specify the grid, Stokes drift, boundary conditions, and
# Coriolis parameter.
#
# ### Domain specification and Grid construction
#
# We create a grid with modest resolution. The grid extent is similar, but not
# exactly the same as that in McWilliams et al. (1997).

using Oceananigans.Grids

grid = RegularCartesianGrid(size=(64, 64, 96), extent=(128, 128, 96))

# ### The Stokes Drift profile
#
# The surface wave Stokes drift profile used in McWilliams et al. (1997)
# corresponds to a 'monochromatic' (that is, single-frequency) wave field with

const wavenumber = 2π / 60 # m⁻¹
nothing # hide

# and

const amplitude = 0.0 # m
nothing # hide

# The `const` declarations ensure that Stokes drift functions compile on the GPU.
# To run this example on the GPU, write `architecture = GPU()` in the constructor
# for `IncompressibleModel` below.
#
# The Stokes drift at the surface for a monochromatic, deep water wave is

using Oceananigans.Buoyancy: g_Earth

const Uˢ = amplitude^2 * wavenumber * sqrt(g_Earth * wavenumber) # m s⁻¹

# The Stokes drift profile is then,

uˢ(z) = Uˢ * exp(2wavenumber * z)
nothing # hide

# which we need for the initial condition.
#
# Note that `Oceananigans.jl` implements the Lagrangian-mean form of the Craik-Leibovich
# equations. This means that our model takes the *vertical derivative* as an input,
# rather than the Stokes drift profile itself.
#
# The vertical derivative of the Stokes drift is

∂z_uˢ(z, t) = 2wavenumber * Uˢ * exp(2wavenumber * z)
nothing # hide

# Adjusting this wind stress to 1e-4 or less allows the simulation to run

Qᵘ = -1e-3 # m² s⁻² approx. Ramudu 30cm/s ice drift (excluding relative motion)
nothing # hide

# and weak cooling with temperature flux
#Qᵇ = 2.307e-9 # m³ s⁻²
Qᵇ = 0
nothing # hide

# Oceananigans uses "positive upward" conventions for all fluxes. In consequence,
# a negative flux at the surface drives positive velocities, and a positive flux of
# buoyancy drives cooling.
#
# The initial condition and bottom boundary condition for buoyancy
# impose a linear stratification with buoyancy frequency

N² = 1.936e-5 # s⁻²
nothing # hide

# To summarize, we impose a surface flux on $u$,

using Oceananigans.BoundaryConditions

# Physical constants.
ρ₀ = 1027  # Density of seawater [kg/m³]
heat_cap = 4000  # Specific heat capacity of seawater at constant pressure [J/(kg·K)]

uᶠ = sqrt(-Qᵘ/ρ₀)

Laᵗ = sqrt(uᶠ/uˢ(0))

# We impose the wind stress Fu as a flux at the surface.
# To impose a flux boundary condition, the top flux imposed should be negative
# for a heating flux and positive for a cooling flux, thus the minus sign on Fθ.
Fu = Qᵘ
Qᵀ = -Qᵇ / (ρ₀*heat_cap)
Qˢ = 0

u_boundary_conditions = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))
nothing # hide

# and a surface flux and bottom linear gradient on buoyancy, $b$,

T_boundary_conditions = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵀ),
                                                       bottom = BoundaryCondition(Gradient, N²))
nothing # hide

S_boundary_conditions = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qˢ),
                                                       bottom = BoundaryCondition(Gradient, N²))
nothing # hide

# ### Coriolis parameter
#
# McWilliams et al. (1997) use

f = 1e-4 # s⁻¹
nothing # hide

FT = Float64

# which is typical for mid-latitudes on Earth.

# ## Model instantiation
#
# Finally, we are ready to build the model. We use the AnisotropicMinimumDissipation
# model for large eddy simulation. Because our Stokes drift does not vary in $x, y$,
# we use `UniformStokesDrift`, which expects Stokes drift functions of $z, t$ only.

using Oceananigans.Buoyancy: BuoyancyTracer
using Oceananigans.SurfaceWaves: UniformStokesDrift
using SeawaterPolynomials.TEOS10
using SeawaterPolynomials: haline_contraction, thermal_expansion

equation_of_state = TEOS10EquationOfState(FT)
EOS = TEOS10.EOS₁₀
β = haline_contraction(-1.5, 27.0, 0.0, equation_of_state)
α = thermal_expansion(-1.5, 27.0, 0.0, equation_of_state)

model = IncompressibleModel(        architecture = CPU(),
                                            grid = grid,
                                        buoyancy = SeawaterBuoyancy(gravitational_acceleration = g_Earth,
                                                    equation_of_state = TEOS10EquationOfState(FT),
                                                    constant_temperature = false, constant_salinity = false),
                                        coriolis = FPlane(f=f),
                                         closure = AnisotropicMinimumDissipation(),
                                   surface_waves = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                             boundary_conditions = (u=u_boundary_conditions,
                                                    T=T_boundary_conditions,
                                                    S=S_boundary_conditions),
                            )

# ## Initial conditions
#
# We make use of random noise concentrated in the upper 4 meters
# for buoyancy and velocity initial conditions,

Ξ(z) = randn() * exp(z / 4)
nothing # hide

# Create a buoyancy profile with a subsurface maximum
# First create Heaviside and Interbal functions
function heaviside(t)
   0.5 * (sign(t) + 1)
end
function interval(t, a, b)
   heaviside(t-a) - heaviside(t-b)
end

# Define initial temperature profile to appoximate Ramudu - first a constant profile
# Add a Gaussian maximum between 15m and 25m depth (peaking at 20m)
# Add a second Gaussian maximum between 40 and 50 m depth (triple amplitude of other peak)
# The second Gaussian is skewed to have a sharp temperature change at its
# upper boundary

η(z) = 10.0 +
    + sqrt(abs(Qᵀ)) * 1e-1 * Ξ(z)
nothing # hide

#Define temperature and convert from Celsius to Kelvin for EOS
Tᵢ(x, y, z) = η(z)
nothing # hide

# Define initial salinity profile (PSW only)
# constant above 40m, increase quadratically to 33 below
#salinity(z) = 27.5 + (0.02/30.)*heaviside(-z-10)*(-z-10) +
#    + 4*heaviside(-z-50)*tanh(-(z+50)/25) +
#    + sqrt(abs(Qˢ)) * 1e-1 * Ξ(z)
#nothing # hide

# Define initial salinity profile
# constant above 40m, increase quadratically to 33 below
salinity(z) = 34 +
    + sqrt(abs(Qˢ)) * 1e-1 * Ξ(z)
nothing # hide

Sᵢ(x, y, z) = salinity(z)
nothing # hide

# The velocity initial condition is zero *Eulerian* velocity. This means that we
# must add the Stokes drift profile to the $u$ velocity field. We also add noise scaled
# by the friction velocity to $u$ and $w$.

uᵢ(x, y, z) = uˢ(z) + sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

wᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

set!(model, u=uᵢ, w=wᵢ, T=Tᵢ, S=Sᵢ)

# ## Setting up the simulation
#
# We use the `TimeStepWizard` for adaptive time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 0.2,

wizard = TimeStepWizard(cfl=0.2, Δt=5.0, max_change=1.1, max_Δt=10.0)
nothing # hide

# ### Nice progress messaging
#
# We define a function that prints a helpful message with
# maximum absolute value of $u, v, w$ and the current wall clock time.

using Oceananigans.Diagnostics, Printf

umax = FieldMaximum(abs, model.velocities.u)
vmax = FieldMaximum(abs, model.velocities.v)
wmax = FieldMaximum(abs, model.velocities.w)

wall_clock = time_ns()

function print_progress(simulation)
    model = simulation.model

    ## Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   model.clock.iteration,
                   prettytime(model.clock.time),
                   prettytime(wizard.Δt),
                   umax(), vmax(), wmax(),
                   prettytime(1e-9 * (time_ns() - wall_clock))
                  )

    @info msg

    return nothing
end

# Now we create the simulation,

using Oceananigans.Utils: hour # correpsonds to "1 hour", in units of seconds

simulation = Simulation(model, progress_frequency = 100,
                                               Δt = wizard,
                                        stop_time = 4hour,
                                         progress = print_progress)

# ## Output
#
# We set up an output writer for the simulation that saves all velocity fields,
# tracer fields, and the subgrid turbulent diffusivity every 2 minutes.

using Oceananigans.OutputWriters
using Oceananigans.Utils: minute

field_outputs = FieldOutputs(merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,)))

simulation.output_writers[:fields] = JLD2OutputWriter(model, field_outputs,
                                                      interval = 2minute,
                                                        prefix = "langmuir_turbulence",
                                                         force = true)
nothing # hide

# ## Running the simulation
#
# This part is easy,

run!(simulation)

# # Making a neat movie
#
# We look at the results by plotting vertical slices of $u$ and $w$, and a horizontal
# slice of $w$ to look for Langmuir cells.

k = searchsortedfirst(grid.zF[:], -8)
nothing # hide

# Making the coordinate arrays takes a few lines of code,

xw, yw, zw = nodes(model.velocities.w)
xu, yu, zu = nodes(model.velocities.u)
xT, yT, zT = nodes(model.tracers.T)
xS, yS, zS = nodes(model.tracers.S)

xw, yw, zw = xw[:], yw[:], zw[:]
xu, yu, zu = xu[:], yu[:], zu[:]
xT, yT, zT = xT[:], yT[:], zT[:]
xS, yS, zS = xS[:], yS[:], zS[:]
nothing # hide

# Next, we open the JLD2 file, and extract the iterations we ended up saving at,

using JLD2, Plots

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

# This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim)
    levels = range(-clim, stop=clim, length=10)

    cmax = maximum(abs, c)

    if clim < cmax # add levels on either end
        levels = vcat([-cmax], range(-clim, stop=clim, length=10), [cmax])
    end

    return levels
end
nothing # hide

# Finally, we're ready to animate.

@info "Making an animation from the saved data..."

anim = @animate for (i, iter) in enumerate(iterations)

    @info "Drawing frame $i from iteration $iter \n"

    ## Load 3D fields from file, omitting halo regions
    w = file["timeseries/w/$iter"][2:end-1, 2:end-1, 2:end-1]
    u = file["timeseries/u/$iter"][2:end-1, 2:end-1, 2:end-1]
    T = file["timeseries/T/$iter"][2:end-1, 2:end-1, 2:end-1]
    S = file["timeseries/S/$iter"][2:end-1, 2:end-1, 2:end-1]

    ## Extract slices
    wxy = w[:, :, k]
    wxz = w[:, 1, :]
    uxz = u[:, 1, :]
    Tz = T[1, 1, :]
    Sz = S[1, 1, :]

    wlim = 0.02
    ulim = 0.05
    Tlim = 0.002
    Slim = 1.5
    wlevels = nice_divergent_levels(w, wlim)
    ulevels = nice_divergent_levels(w, ulim)

    wxy_plot = contourf(xw, yw, wxy';
                              color = :balance,
                        aspectratio = :equal,
                              clims = (-wlim, wlim),
                             levels = wlevels,
                              xlims = (0, grid.Lx),
                              ylims = (0, grid.Ly),
                             xlabel = "x (m)",
                             ylabel = "y (m)")

    wxz_plot = contourf(xw, zw, wxz';
                              color = :balance,
                        aspectratio = :equal,
                              clims = (-wlim, wlim),
                             levels = wlevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    uxz_plot = contourf(xu, zu, uxz';
                              color = :balance,
                        aspectratio = :equal,
                              clims = (-ulim, ulim),
                             levels = ulevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    Tz_plot = plot(Tz, zT;
                            xlims = (-1.6, 0.2),
                            ylims = (-grid.Lz, 0),
                             xlabel = "Pot. Temp. (C)",
                             ylabel = "z (m)",
                             legend = false)

    Sz_plot = plot(Sz, zS;
                            xlims = (27, 33),
                            ylims = (-grid.Lz, 0),
                            xlabel = "S (g/kg)",
                            ylabel = "z (m)",
                            legend = false)

    l = @layout [a b c; d e]
    plot(wxy_plot, wxz_plot, uxz_plot, Tz_plot, Sz_plot, layout=l, size=(1300, 800),
         title = ["w(x, y, z=-8, t) (m/s)" "w(x, y=0, z, t) (m/s)" "u(x, y = 0, z, t) (m/s)" "T(x=0, y = 0, z, t) (m/s)" "S(x=0, y = 0, z, t) (m/s)"])

    iter == iterations[end] && close(file)
end

mp4(anim, "subsurface_heat_langmuir_HR.mp4", fps = 5) # hide


text:
function cr3bpEomControlMTKSys(μ::AbstractFloat)
    @parameters atx aty atz
    @variables t x(t) y(t) z(t) vx(t) vy(t) vz(t)
    D = Differential(t)

    # Compute r1 and r2
    r1 = sqrt((x+μ)^2 + y^2 + z^2)
    r2 = sqrt((x+μ-1)^2 + y^2 + z^2)

    # Compute g
    gx = x - (1 - μ)*(x + μ)/r1^3 - μ*(x + μ - 1)/r2^3
    gy = y - (1 - μ)*y/r1^3 - μ*y/r2^3
    gz = -(1 - μ)*z/r1^3 - μ*z/r2^3

    # Compute h
    hx = 2*vy
    hy = -2*vx

    eqs = [D(x) ~ vx,
           D(y) ~ vy,
           D(z) ~ vz,
           D(vx) ~ gx + hx + atx,
           D(vy) ~ gy + hy + aty,
           D(vz) ~ gz + atz]

    return ODESystem(eqs)
end

function cr3bpEomNoControlMTKSys(μ::AbstractFloat)
    _cr3bpEomWControl = cr3bpEomControlMTKSys(μ)
    (atx, aty, atz) = parameters(_cr3bpEomWControl)
    eqs = [eq.lhs ~ simplify(
        substitute(substitute(substitute(
            eq.rhs, atx => 0), aty => 0), atz => 0)) for
            eq in equations(_cr3bpEomWControl)]
    return ODESystem(eqs)
end

function cr3bpEOMIndirect(μ::AbstractFloat,    TU::AbstractFloat,
                          LU::AbstractFloat,   MU::AbstractFloat,
                          tMax::AbstractFloat, c::AbstractFloat)
    @parameters u
    @variables t x(t) y(t) z(t) vx(t) vy(t) vz(t) m(t)
    @variables λx(t) λy(t) λz(t) λvx(t) λvy(t) λvz(t) λm(t) H(t)
    D    = Differential(t)

    # Compute r1 and r2
    r1 = sqrt((x+μ)^2 + y^2 + z^2)
    r2 = sqrt((x+μ-1)^2 + y^2 + z^2)

    # Compute g
    gx = x - (1 - μ)*(x + μ)/r1^3 - μ*(x + μ - 1)/r2^3
    gy = y - (1 - μ)*y/r1^3 - μ*y/r2^3
    gz = -(1 - μ)*z/r1^3 - μ*z/r2^3

    # Compute h
    hx = 2*vy
    hy = -2*vx

    # Compute G
    G11     = 1 - (1 - μ)/r1^3 + 3*(1 - μ)*(x + μ)^2/r1^5 -
                μ/r2^3 + 3*μ*(x + μ - 1)^2/r2^5;
    G22     = 1 - (1 - μ)/r1^3 + 3*(1 - μ)*y^2/r1^5 -
                μ/r2^3 + 3*μ*y^2/r2^5;
    G33     = -(1 - μ)/r1^3 + 3*(1 - μ)*z^2/r1^5 -
                μ/r2^3 + 3*μ*z^2/r2^5;
    G12     = 3*(1 - μ)*(x + μ)*y/r1^5 +
                3*μ*(x + μ - 1)*y/r2^5;
    G13     = 3*(1 - μ)*(x + μ)*z/r1^5 +
                3*μ*(x + μ - 1)*z/r2^5;
    G23     = 3*(1 - μ)*y*z/r1^5 +
                3*μ*y*z/r2^5;

    # Compute scaled max thrust and exhaust velocity
    tMaxSc = tMax * TU * TU / (MU*LU*1000.0)
    cSc = c * TU / (LU*1000.0)

    # Compute thrust direction
    λv = sqrt(λvx^2 + λvy^2 + λvz^2)
    αx = -λvx / λv
    αy = -λvy / λv
    αz = -λvz / λv

    # Compute thrust acceleration
    atx = u*tMaxSc*αx / m
    aty = u*tMaxSc*αy / m
    atz = u*tMaxSc*αz / m

    # Dynamics
    dx   = vx
    dy   = vy
    dz   = vz
    ddx  = gx + hx + atx
    ddy  = gy + hy + aty
    ddz  = gz + atz
    dm   = -u*tMaxSc / cSc
    dλx  = -G11*λvx - G12*λvy - G13*λvz
    dλy  = -G12*λvx - G22*λvy - G23*λvz
    dλz  = -G13*λvx - G23*λvy - G33*λvz
    dλvx = -λx + 2.0*λvy
    dλvy = -λy - 2.0*λvx
    dλvz = -λz
    dλm  = -λv*u*tMaxSc / (m^2)

    # Equations
    eqs = [D(x)   ~ dx,
           D(y)   ~ dy,
           D(z)   ~ dz,
           D(vx)  ~ ddx,
           D(vy)  ~ ddy,
           D(vz)  ~ ddz,
           D(m)   ~ dm,
           D(λx)  ~ dλx,
           D(λy)  ~ dλy,
           D(λz)  ~ dλz,
           D(λvx) ~ dλvx,
           D(λvy) ~ dλvy,
           D(λvz) ~ dλvz,
           D(λm)  ~ dλm]

    return ODESystem(eqs, t)
end


text:
using Pkg
Pkg.add("PkgTemplates")
using PkgTemplates

templ = Template(;
        user="sambuddhac",
        authors=["Sambuddha Chakrabarti", "Hosna Khajeh", "Tom Nudell", "Mohammad Reza Hesamzadeh", "Ross Baldick"],
        dir="/Users/sc87/code/JuliaJuMPTestOpt",
        julia=v"1.0.0",
        host="github.com",
        plugins=[
            ProjectFile(;
                version=v"0.1.0"),
            License(; name="MIT", path=nothing, destination="LICENSE"),
            Readme(;
                file="/Users/sc87/.julia/packages/PkgTemplates/60BEI/templates/README.md",
                destination="README.md",
                inline_badges=true,
            ),
            Tests(; file="/Users/sc87/.julia/packages/PkgTemplates/60BEI/templates/test/runtests.jl", project=false),
            Documenter{TravisCI}(;
                make_jl="/Users/sc87/.julia/packages/PkgTemplates/60BEI/templates/docs/make.jl",
                index_md="/Users/sc87/.julia/packages/PkgTemplates/60BEI/templates/docs/src/index.md",
                assets=String[],
                logo=Logo(),
                #canonical_url=make_canonical(TravisCI),
                makedocs_kwargs=Dict{Symbol, Any}(),
            ),
            TravisCI(),
            Codecov(),
            Coveralls(),
            AppVeyor(),
            CirrusCI(),
            DroneCI(),
            GitLabCI(),
            Logo(; light=nothing, dark=nothing),
            ColPracBadge(),
        ],
    )

templ("TestOptJulia")

max_stars_repo_name: sambuddhac/TestOptJulia

max_issues_repo_name: sambuddhac/TestOptJulia

max_forks_repo_name: sambuddhac/TestOptJulia


text:
"""
    project!(fun, grid, i, cⁱ)

Project the function `fun` onto the finite element of `grid` with index
`i`, writing the results in the vector of expansion coefficients `cⁱ`.
"""
function project!(fun, grid::Grid, i::Integer,
                  cⁱ::AbstractVector)
    n = order(grid)
    cⁱ[1] = grid.N[i,1]*grid.W[i,1]*fun(grid.X[i,1])
    for m = 2:n-1
        cⁱ[m] = √(grid.W[i,m])*fun(grid.X[i,m])
    end
    cⁱ[end] = grid.N[i,end]*grid.W[i,end]*fun(grid.X[i,end])
end

"""
    project(fun, basis)

Project the function `fun` onto the FEDVR `basis`.
"""
function project(fun, basis::Basis)
    grid = basis.grid
    C = zero(grid.X)
    for i = elems(grid)
        project!(fun, grid, i, view(C, i, :))
    end
    # The bridge functions overlap the function both in element i and
    # element i+1; thus the contributions need to be added into one
    # expansion coefficient.
    C[2:elcount(grid),1] += C[1:elcount(grid)-1,end]

    boundary_sel(grid, [C[:,1:end-1]'[:]..., C[end]])
end

export project


text:
function mark_cross!(r::Robot)
    a = 0
    b = 0
    A = 0
    B = 0
    k = 0
    while !isborder(r, Nord)
        move!(r, Nord)
        a = a + 1
    end
    while !isborder(r, Sud)
        A = A + 1
        move!(r, Sud)
    end
    a = A - a
    while !isborder(r, Ost)
        b = b + 1
        move!(r, Ost)
    end
    putmarker!(r)
    while !isborder(r, West)
        B = B + 1
        move!(r, West)
        putmarker!(r)
    end
    b = B - b
    B = B - 1
    k = B
    while A > 0
        if (k = 0)
            break
        end
        move!(r, Nord)
        putmarker!(r)
        while B > 0
            move!(r, Ost)
            putmarker!(r)
        end
        B = k
        while B > 0
            move!(r, West)
            B = B - 1
        end
        k = k - 1
        B = k
        A = A - 1
    end
end


text:
## Adapted from sundialsTB/kinsol/examples_ser/mkinTest_nds.m

## %mkinTest_dns - KINSOL example problem (serial, dense)
## %   Simple test problem for the Dense linear solver in KINSOL
## %   This example solves the system
## %       y(1)^2 + y(2)^2 = 1
## %       y(2) = y(1)^2
## %

## % Radu Serban <radu@llnl.gov>

using Sundials

## function to be optimized
function sysfn(y_nv, fy_nv, a_in)
    y = convert(Vector, y_nv)
    fy = convert(Vector, fy_nv)
    fy[1] = y[1]^2 + y[2]^2 - 1.0
    fy[2] = y[2] - y[1]^2
    return Sundials.KIN_SUCCESS
end

## Initialize problem
const neq = 2
kmem = Sundials.KINCreate()
Sundials.@checkflag Sundials.KINSetFuncNormTol(kmem, 1.0e-5)
Sundials.@checkflag Sundials.KINSetScaledStepTol(kmem, 1.0e-4)
Sundials.@checkflag Sundials.KINSetMaxSetupCalls(kmem, 1)
y = ones(neq)
Sundials.@checkflag Sundials.KINInit(kmem, sysfn, y)
Sundials.@checkflag Sundials.KINDense(kmem, neq)
## Solve problem
scale = ones(neq)
Sundials.@checkflag Sundials.KINSol(kmem, y, Sundials.KIN_LINESEARCH,
                                    scale, scale)

println("Solution: ", y)
residual = ones(2)
sysfn(y, residual, [1, 2])
println("Residual: ", residual)


text:
# example of a unit test in Julia
# runs a test for a certain module
# to run the tests, first open the package manager (`]` in REPL),
# activate the project if not done so and then enter `test`

# wrap all your tests and subgroups in a `@testset` block
@testset "BeesAlgorithm" begin
    using STMOZOO.BeesAlgorithm  # load YOUR module
    using STMOZOO.BeesAlgorithm: initialize_population, compute_objective, compute_fitness, employed_bee_phase, scouting_phase, onlooker_bee_phase, foodsource_info_prob, create_newsolution

    S = 24
    bounds_lower = [-100,-100,-100,-100]
    bounds_upper = [100,100,100,100]
    D = 4
    limit = D * (S/2)
    T = 500
    n = Int8(S/2)

    @testset "initialize_population" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        @test length(population) == n
        @test all([bounds_lower <= population[i] <= bounds_upper for i in 1:length(population)])

    end

    @testset "compute_objective" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        objectives = compute_objective(population,sphere)
        @test length(objectives) == length(population)

    end

    @testset "compute_fitness" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        objectives = compute_objective(population,sphere)
        fitness = compute_fitness(objectives)
        @test length(fitness) == length(objectives) == length(population)

    end

    @testset "foodsource_info_prob" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        objectives = compute_objective(population,sphere)
        fitness = compute_fitness(objectives)
        probab = foodsource_info_prob(fitness)
        @test length(probab) == length(fitness)
        @test all(0 .≤ probab .≤ 1)

    end

    @testset "create_newsolution" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        solution = [4, 0, 1, 4]
        new_solution = create_newsolution(solution, population, bounds_lower, bounds_upper)
        @test length(new_solution) == length(solution)
        @test sum(new_solution .==  solution) == D-1

    end

    @testset "employed_bee_phase" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        trial = zeros(size(population)[1])
        population_new_evolved, fitness_new_evolved, objective_new_evolved, newtrial = employed_bee_phase(population, bounds_lower, bounds_upper, trial,n ,sphere)

        @test length(population_new_evolved) == length(population)
        @test length(newtrial) == length(trial)
    end

    @testset "onlooker_bee_phase" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        trial = zeros(size(population)[1])
        population_new_evolved, fitness_new_evolved, objective_new_evolved, newtrial = onlooker_bee_phase(population, bounds_lower, bounds_upper, trial,n ,sphere)

        @test length(population_new_evolved) == length(population)
        @test length(newtrial) == length(trial)
    end

    @testset "scouting_phase" begin

        population = initialize_population(D, bounds_lower, bounds_upper, n)
        trial = ones(size(population)[1])
        objective= compute_objective(population,sphere)
        fitness = compute_fitness(objective)
        population_new_evolved, fitness_new_evolved, objective_new_evolved, newtrial  = scouting_phase(population, bounds_lower, bounds_upper,D, trial, fitness, objective, 0, sphere)

        @test length(objective) == length(fitness)
        @test length(population_new_evolved) == length(population)
        @test length(newtrial) == length(trial)
    end

    @testset "ArtificialBeeColonization" begin

        optimal_solution,populations,best_fitness_tracker = ArtificialBeeColonization(D, bounds_lower, bounds_upper, S, T, limit, sphere)

        @test length(optimal_solution) == length(bounds_lower) == length(bounds_upper) == D
        @test optimal_solution isa Vector
        @test optimal_solution >= bounds_lower
        @test optimal_solution <= bounds_upper

    end

end

max_stars_repo_name: kirstvh/STMOZOO

max_issues_repo_name: kirstvh/STMOZOO

max_forks_repo_name: kirstvh/STMOZOO


text:
 @testset "inventory.jl" begin
    inv = Inventory(3, 10)
    @test state(inv) == [inv.onhand] == [10]
    @test length(state(inv)) == state_size(inv) == 1
    @test InventoryModels.inventory_position(inv) == 10
    dest = []
    InventoryModels.pull!(inv, 7, dest)
    @test only(inv.pull_orders) == (dest => 7)
    InventoryModels.dispatch!(inv)
    @test InventoryModels.reward!(inv) == -3*3
    @test state(inv) == [3]
    @test dest == [7, inv]
    InventoryModels.pull!(inv, 7, dest)
    @test only(inv.pull_orders) == (dest => 7)
    InventoryModels.dispatch!(inv)
    @test InventoryModels.reward!(inv) == 0
    @test state(inv) == [0]
    @test dest == [7, inv, 3, inv]
    reset!(inv)
    @test state(inv) == [inv.onhand] == [10]
    @test InventoryModels.inventory_position(inv) == 10
    dest2 = []
    InventoryModels.pull!(inv, 10, dest)
    InventoryModels.pull!(inv, 5, dest2)
    @test length(inv.pull_orders) == 2
    InventoryModels.dispatch!(inv)
    @test InventoryModels.reward!(inv) == 0
    @test state(inv) == [0]
    @test dest == [7, inv, 3, inv, 2/3*10, inv]
    @test dest2 == [2/3*5, inv]

    inv = Inventory(3, Uniform(10,20), capacity = 10)
    @test state(inv) == [10]
    push!(inv, 10, dest)
    @test state(inv) == [10]
    inv.capacity = 20
    reset!(inv)
    @test 10 < inv.onhand < 20

end


text:
# using LinearAlgebra, StatsModels, Random, Distributions, MixedModelsBLB, DataFrames

# # Run blb
# feformula   = @formula(y ~ 1 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 +
#                             x11 + x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 +
#                             x21 + x22 + x23 + x24 + x25 + x26 + x27 + x28 + x29 + x30 +
#                             x31 + x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 + x40 +
#                             x41 + x42 + x43 + x44 + x45 + x46 + x47 + x48 + x49 + x50 +
#                             x51 + x52 + x53 + x54 + x55 + x56 + x57 + x58 + x59 + x60 +
#                             x61 + x62 + x63 + x64 + x65 + x66 + x67 + x68 + x69 + x70 +
#                             x71 + x72 + x73 + x74 + x75 + x76 + x77 + x78 + x79 + x80 +
#                             x81 + x82 + x83 + x84 + x85 + x86 + x87 + x88 + x89 + x90 +
#                             x91 + x92 + x93 + x94 + x95 + x96 + x97 + x98 + x99)
# reformula   = @formula(y ~ 1 + z1)

# # datatable = JuliaDB.loadtable("test/data/testdata-N-20000-REPS-10.csv")
# datatable = DataFrame(CSV.read("test/data/testdata-N-20000-REPS-10.csv"))
# blb_full_data(
#         MersenneTwister(1),
#         datatable;
#         feformula = feformula,
#         reformula = reformula,
#         id_name = "id", cat_names = Array{String,1}(),
#         subset_size = 1000, n_subsets = 1, n_boots = 10,
#         solver = Ipopt.IpoptSolver(print_level=5, max_iter=100, mehrotra_algorithm = "yes", warm_start_init_point = "yes"),
#         # solver = Ipopt.IpoptSolver(
#         #   print_level = 5, derivative_test = "second-order", derivative_test_print_all = "yes", check_derivatives_for_naninf = "yes"
#         # ),
#         verbose = true, use_threads = false, use_groupby = true, nonparametric_boot = true
# )

# # # Generate Data
# # # Global parameters
# # N = 20000 # number of individuals
# # reps = 10 # number of observations from each individual
# # # id = repeat(1:N, inner = reps)
# # p = 100 # number of fixed effects
# # q = 2 # number of random effects
# # use_threads = false

# # # true parameter values
# # βtrue  = ones(100) #[0.1; 6; -3; 1]
# # σ²true = 1
# # σtrue  = sqrt(σ²true)
# # Σtrue  = [1.0 1; 1 4.0] # Matrix(Diagonal([2.0; 1.2; 1.0]))
# # Ltrue  = Matrix(cholesky(Symmetric(Σtrue)).L)

# # obsvec = Vector{blblmmObs{Float64}}(undef, N)
# # # initialize arrays
# # X = Matrix{Float64}(undef, reps, p)
# # X[:, 1] = ones(reps)
# # Z = Matrix{Float64}(undef, reps, q)
# # Z[:, 1] = ones(reps)
# # storage_q = Vector{Float64}(undef, q)
# # re_storage = Vector{Float64}(undef, q)
# # y = Vector{Float64}(undef, reps)
# # fenames = vcat("Intercept", "x" .* string.([1:1:(p-1);]))
# # renames = ["Intercept", "z1"]

# # # Generate data
# # for i in 1:N
# #     randn!(y) # y = standard normal error
# #     # first column intercept, remaining entries iid std normal
# #     @views randn!(X[:, 2:p]) #Distributions.rand!(Normal(), X[:, 2:p])
# #     BLAS.gemv!('N', 1., X, βtrue, σtrue, y) # y = Xβ + σtrue * standard normal error
# #     randn!(storage_q)
# #     BLAS.gemv!('N', 1., Ltrue, storage_q, 0., re_storage)
# #     # first column intercept, remaining entries iid std normal
# #     @views randn!(Z[:, 2:q]) #Distributions.rand!(Normal(), Z[:, 2:q])
# #     BLAS.gemv!('N', 1., Z, re_storage, 1., y) # y = Xβ + Zα + error
# #     # y = X * βtrue .+ Z * (Ltrue * randn(q)) .+ σtrue * randn(ns[i])
# #     # form a blblmmObs instance
# #     obsvec[i] = blblmmObs(copy(y), copy(X), copy(Z))
# # end
# # # form a LmmModel instance
# # lmm = blblmmModel(obsvec, fenames, renames, N, use_threads)

# # # Save the data to a CSV file.
# # open(string("test/data/testdata-N-", N, "-REPS-", reps, ".csv"), "w") do io
# #     p = size(lmm.data[1].X, 2)
# #     q = size(lmm.data[1].Z, 2)
# #     # print header
# #     print(io, "id,y,")
# #     for j in 1:(p-1)
# #         print(io, "x" * string(j) * ",")
# #     end
# #     for j in 1:(q-1)
# #         print(io, "z" * string(j) * (j < q-1 ? "," : "\n"))
# #     end
# #     # print data
# #     for i in eachindex(lmm.data)
# #         obs = lmm.data[i]
# #         for j in 1:length(obs.y)
# #             # id
# #             print(io, i, ",")
# #             # y
# #             print(io, obs.y[j], ",")
# #             # X data
# #             for k in 2:p
# #                 print(io, obs.X[j, k], ",")
# #             end
# #             # Z data
# #             for k in 2:q-1
# #                 print(io, obs.Z[j, k], ",")
# #             end
# #             print(io, obs.Z[j, q], "\n")
# #         end
# #     end
# # end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/relax.jl
# Defines routines used construct the relaxed subproblem.
#############################################################################

"""
$(FUNCTIONNAME)

Applies the safe cut checks detailed in Khajavirad, 2018 [Khajavirad, Aida,
and Nikolaos V. Sahinidis. "A hybrid LP/NLP paradigm for global optimization
relaxations." Mathematical Programming Computation 10.3 (2018): 383-421] to
ensure that only numerically safe affine relaxations are added. Checks that
i) ``|b| <= safe b`, ii) `safe_l <= abs(ai) <= safe u`, and iii) violates
`safe_l <= abs(ai/aj) <= safe_u`.
"""
function is_safe_cut!(m::Optimizer, f::SAF)

    safe_l = m._parameters.cut_safe_l
    safe_u = m._parameters.cut_safe_u
    safe_b = m._parameters.cut_safe_b

    # violates |b| <= safe_b
    (abs(f.constant) > safe_b) && return false

    term_count = length(f.terms)
    for i = 1:term_count

        ai = (@inbounds f.terms[i]).coefficient
        if ai !== 0.0

            # violates safe_l <= abs(ai) <= safe_u
            ai_abs = abs(ai)
            !(safe_l <= abs(ai) <= safe_u) && return false

            # violates safe_l <= abs(ai/aj) <= safe_u
            for j = i:term_count
                aj = (@inbounds f.terms[j]).coefficient
                if aj !== 0.0
                    !(safe_l <= abs(ai/aj) <= safe_u) && return false
                end
            end
        end
    end

    return true
end

"""
$(FUNCTIONNAME)

Relaxs the constraint by adding an affine constraint to the model.
"""
function relax! end

"""
$(FUNCTIONNAME)

Default routine for relaxing quadratic constraint `func < 0.0` on node `n`.
Takes affine bounds of convex part at point `x0` and secant line bounds on
concave parts.
"""
function affine_relax_quadratic!(func::SQF, buffer::Dict{Int,Float64}, saf::SAF,
                                 n::NodeBB, sol_to_branch_map::Vector{Int},
                                 x::Vector{Float64})

    lower_bounds = n.lower_variable_bounds
    upper_bounds = n.upper_variable_bounds
    quadratic_constant = func.constant

    # Affine terms only contribute coefficients, so the respective
    # values do not contribute to the cut. Since all quadratic terms
    # are considered to be branch variables we exclude any potential
    # need to retrieve variable bounds from locations other than
    # the node.
    for term in func.quadratic_terms

        a = term.coefficient
        idx1 = term.variable_index_1.value
        idx2 = term.variable_index_2.value
        sol_idx1 = sol_to_branch_map[idx1]
        sol_idx2 = sol_to_branch_map[idx2]
        x0_1 = x[sol_idx1]
        xL_1 = lower_bounds[sol_idx1]
        xU_1 = upper_bounds[sol_idx1]

        if idx1 === idx2

            if a > 0.0
                buffer[idx1] += a*x0_1
                quadratic_constant -= 0.5*a*x0_1*x0_1

            else
                if !isinf(xL_1) && !isinf(xU_1)
                    buffer[idx1] += 0.5*a*(xL_1 + xU_1)
                    quadratic_constant -= 0.5*a*xL_1*xU_1
                else
                    return false
                end
            end

        else
            x0_2 = x[sol_idx2]
            xL_2 = lower_bounds[sol_idx2]
            xU_2 = upper_bounds[sol_idx2]

            if a > 0.0
                if (!isinf(xL_1) && !isinf(xL_2)) &&
                   ((xU_1 - xL_1)*x0_2 + (xU_2 - xL_2)*x0_1 <= xU_1*xU_2 - xL_1*xL_2)
                    buffer[idx1] += a*xL_2
                    buffer[idx2] += a*xL_1
                    quadratic_constant -= a*xL_1*xL_2

                elseif !isinf(xU_1) && !isinf(xU_2)
                    buffer[idx1] += a*xU_2
                    buffer[idx2] += a*xU_1
                    quadratic_constant -= a*xU_1*xU_2

                else
                    return false

                end
            else
                if (!isinf(xU_1) && !isinf(xL_2)) &&
                   ((xU_1 - xL_1)*x0_2 - (xU_2 - xL_2)*x0_1 <= xU_1*xL_2 - xL_1*xU_2)

                    buffer[idx1] += a*xL_2
                    buffer[idx2] += a*xU_1
                    quadratic_constant -= a*xU_1*xL_2

                elseif !isinf(xL_1) && !isinf(xU_2)
                    buffer[idx1] += a*xU_2
                    buffer[idx2] += a*xL_1
                    quadratic_constant -= a*xL_1*xU_2

                else
                    return false
                end
            end
        end
    end

    for term in func.affine_terms
        a0 = term.coefficient
        idx = term.variable_index.value
        buffer[idx] += a0
    end

    count = 1
    for (key, value) in buffer
        saf.terms[count] = SAT(value, VI(key))
        buffer[key] = 0.0
        count += 1
    end
    saf.constant = quadratic_constant

    return true
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::Optimizer, f::BufferedQuadraticIneq, indx::Int, check_safe::Bool)

    constraint_tol = m._parameters.absolute_constraint_feas_tolerance
    finite_cut_generated = affine_relax_quadratic!(f.func, f.buffer, f.saf, m._current_node, m._sol_to_branch_map, m._current_xref)
    if finite_cut_generated
        if !check_safe || is_safe_cut!(m, f.saf)
            lt = LT(-f.saf.constant + constraint_tol)
            f.saf.constant = 0.0
            ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
            push!(m._buffered_quadratic_ineq_ci, ci)
        end
    end
    #m.relaxed_to_problem_map[ci] = indx

    return nothing
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::Optimizer, f::BufferedQuadraticEq, indx::Int, check_safe::Bool)

    constraint_tol = m._parameters.absolute_constraint_feas_tolerance
    finite_cut_generated = affine_relax_quadratic!(f.func, f.buffer, f.saf, m._current_node, m._sol_to_branch_map, m._current_xref)
    if finite_cut_generated
        if !check_safe || is_safe_cut!(m, f.saf)
            lt = LT(-f.saf.constant + constraint_tol)
            f.saf.constant = 0.0
            ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
            push!(m._buffered_quadratic_eq_ci, ci)
        end
    end
    #m.relaxed_to_problem_map[ci] = indx

    finite_cut_generated = affine_relax_quadratic!(f.minus_func, f.buffer, f.saf, m._current_node, m._sol_to_branch_map, m._current_xref)
    if finite_cut_generated
        if !check_safe || is_safe_cut!(m, f.saf)
            lt = LT(-f.saf.constant + constraint_tol)
            f.saf.constant = 0.0
            ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
            push!(m._buffered_quadratic_eq_ci, ci)
        end
    end
    #m.relaxed_to_problem_map[ci] = indx

    return nothing
end

"""
$(FUNCTIONNAME)
"""
function affine_relax_nonlinear!(f::BufferedNonlinearFunction{MC{N,T}}, evaluator::Evaluator,
                                 use_cvx::Bool, new_pass::Bool, is_constraint::Bool) where {N,T<:RelaxTag}

    if new_pass
        forward_pass!(evaluator, f)
    end
    x = evaluator.x
    finite_cut = true

    expr = f.expr
    grad_sparsity = expr.grad_sparsity
    if expr.isnumber[1]
        f.saf.constant = expr.numberstorage[1]
        for i = 1:N
            vval = @inbounds grad_sparsity[i]
            f.saf.terms[i] = SAT(0.0, VI(vval))
        end

    else
        setvalue = expr.setstorage[1]
        finite_cut &= !(isempty(setvalue) || isnan(setvalue))

        if finite_cut
            value = f.expr.setstorage[1]
            f.saf.constant = use_cvx ? value.cv : -value.cc
            for i = 1:N
                vval = @inbounds grad_sparsity[i]
                if use_cvx
                    coef = @inbounds value.cv_grad[i]
                else
                    coef = @inbounds -value.cc_grad[i]
                end
                f.saf.terms[i] = SAT(coef, VI(vval))
                xv = @inbounds x[vval]
                f.saf.constant = sub_round(f.saf.constant , mul_round(coef, xv, RoundUp), RoundDown)
            end
            if is_constraint
                bnd_used =  use_cvx ? -f.upper_bound : f.lower_bound
                f.saf.constant = add_round(f.saf.constant, bnd_used, RoundDown)
            end
        end
    end

    return finite_cut
end

"""
$(TYPEDSIGNATURES)
"""
function check_set_affine_nl!(m::Optimizer, f::BufferedNonlinearFunction{MC{N,T}}, finite_cut_generated::Bool, check_safe::Bool) where {N,T<:RelaxTag}

    constraint_tol = m._parameters.absolute_constraint_feas_tolerance
    if finite_cut_generated
        if !check_safe || is_safe_cut!(m, f.saf)
            lt = LT(-f.saf.constant + constraint_tol)
            f.saf.constant = 0.0
            ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
            push!(m._buffered_nonlinear_ci, ci)
        end
    end

    return nothing
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::Optimizer, f::BufferedNonlinearFunction{MC{N,T}}, indx::Int, check_safe::Bool) where {N,T<:RelaxTag}
    evaluator = m._working_problem._relaxed_evaluator

    finite_cut_generated = affine_relax_nonlinear!(f, evaluator, true, true, true)
    check_set_affine_nl!(m, f, finite_cut_generated, check_safe)

    finite_cut_generated = affine_relax_nonlinear!(f, evaluator, false, false, true)
    check_set_affine_nl!(m, f, finite_cut_generated, check_safe)

    return nothing
end

"""
$(TYPEDSIGNATURES)
"""
function bound_objective(t::ExtensionType, m::Optimizer)

    n = m._current_node
    sb_map = m._sol_to_branch_map
    wp = m._working_problem
    obj_type = wp._objective_type

    if obj_type === NONLINEAR

        # assumes current node has already been loaded into evaluator
        objective_lo = lower_interval_bound(m, wp._objective_nl, n)

    elseif obj_type === SINGLE_VARIABLE
        obj_indx = @inbounds sb_map[wp._objective_sv.variable.value]
        objective_lo = @inbounds n.lower_variable_bounds[obj_indx]

    elseif obj_type === SCALAR_AFFINE
        objective_lo = lower_interval_bound(wp._objective_saf_parsed, n)

    elseif obj_type === SCALAR_QUADRATIC
        objective_lo = lower_interval_bound(wp._objective_sqf, n)

    end

    return objective_lo
end
bound_objective(m::Optimizer) = bound_objective(m.ext_type, m)

"""
$(TYPEDSIGNATURES)
"""
function relax_objective_nonlinear!(m::Optimizer, wp::ParsedProblem, check_safe::Bool)

    relaxed_optimizer = m.relaxed_optimizer
    relaxed_evaluator = wp._relaxed_evaluator
    buffered_nl = wp._objective_nl

    new_flag = m._new_eval_objective
    relaxed_evaluator.is_first_eval = new_flag
    finite_cut_generated = affine_relax_nonlinear!(buffered_nl, relaxed_evaluator, true, new_flag, false)
    relaxed_evaluator.is_first_eval = false

    if finite_cut_generated
        if !check_safe || is_safe_cut!(m, buffered_nl.saf)
            copyto!(wp._objective_saf.terms, buffered_nl.saf.terms)
            wp._objective_saf.constant = buffered_nl.saf.constant
            MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SAF}(), wp._objective_saf)
        end
    end

    return nothing
end

"""
$(TYPEDSIGNATURES)

Triggers an evaluation of the objective function and then updates
the affine relaxation of the objective function.
"""
function relax_objective!(t::ExtensionType, m::Optimizer, q::Int64)

    relaxed_optimizer = m.relaxed_optimizer
    m._working_problem._relaxed_evaluator

    # Add objective
    wp = m._working_problem
    obj_type = wp._objective_type
    check_safe = (q === 1) ? false : m._parameters.cut_safe_on

    if obj_type === SINGLE_VARIABLE
        MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SV}(), wp._objective_sv)

    elseif obj_type === SCALAR_AFFINE
        MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SAF}(), wp._objective_saf)

    elseif obj_type === SCALAR_QUADRATIC
        buffered_sqf = wp._objective_sqf
        finite_cut_generated = affine_relax_quadratic!(buffered_sqf.func, buffered_sqf.buffer, buffered_sqf.saf,
                                m._current_node, m._sol_to_branch_map, m._current_xref)
        if finite_cut_generated
            if !check_safe || is_safe_cut!(m, buffered_sqf.saf)
                copyto!(wp._objective_saf.terms, buffered_sqf.saf.terms)
                wp._objective_saf.constant = buffered_sqf.saf.constant
                MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SAF}(), wp._objective_saf)
            end
        end

    elseif obj_type === NONLINEAR
        relax_objective_nonlinear!(m, wp, check_safe)
    end

    m._new_eval_objective = false

    return nothing
end
relax_objective!(m::Optimizer, q::Int64) = relax_objective!(m.ext_type, m, q)

"""

Triggers an evaluation of the nonlinear objective function (if necessary)
and adds the corresponding `<a,x> <= b` objective cut.
"""
function objective_cut_nonlinear!(m::Optimizer, wp::ParsedProblem, UBD::Float64, check_safe::Bool)

    relaxed_optimizer = m.relaxed_optimizer
    relaxed_evaluator = wp._relaxed_evaluator
    buffered_nl = wp._objective_nl

    # if the objective cut is the first evaluation of the objective expression
    # then perform a a forward pass
    new_flag = m._new_eval_objective
    relaxed_evaluator.is_first_eval = new_flag
    finite_cut_generated = affine_relax_nonlinear!(buffered_nl, relaxed_evaluator, true, new_flag, false)

    constraint_tol = m._parameters.absolute_constraint_feas_tolerance
    if finite_cut_generated
        copyto!(wp._objective_saf.terms, buffered_nl.saf.terms)
        wp._objective_saf.constant = 0.0
        if !check_safe || is_safe_cut!(m,  buffered_nl.saf)
            # TODO: When we introduce numerically safe McCormick operators we'll need to replace
            # the UBD - buffered_nl.saf.constant with a correctly rounded version. For now,
            # a small factor is added to the UBD calculation initially which should be sufficient.
            ci_saf = MOI.add_constraint(m.relaxed_optimizer, wp._objective_saf, LT(UBD - buffered_nl.saf.constant + constraint_tol))
            push!(m._objective_cut_ci_saf, ci_saf)
        end
    end

    m._new_eval_objective = false

    return nothing
end

"""
$(FUNCTIONNAME)

Adds linear objective cut constraint to the `x.relaxed_optimizer`.
"""
function objective_cut!(m::Optimizer, check_safe::Bool)

    UBD = m._global_upper_bound
    constraint_tol = m._parameters.absolute_constraint_feas_tolerance
    if m._parameters.objective_cut_on && m._global_upper_bound < Inf

        wp = m._working_problem
        obj_type = wp._objective_type

        if obj_type === SINGLE_VARIABLE
            if !isinf(UBD) && (m._objective_cut_ci_sv.value === -1)
                m._objective_cut_ci_sv = CI{SV,LT}(wp._objective_sv.variable.value)
                MOI.set(m.relaxed_optimizer, MOI.ConstraintSet(), m._objective_cut_ci_sv, LT(UBD))
            else
                MOI.set(m.relaxed_optimizer, MOI.ConstraintSet(), m._objective_cut_ci_sv, LT(UBD))
            end

        elseif obj_type === SCALAR_AFFINE
            formulated_constant = wp._objective_saf.constant
            wp._objective_saf.constant = 0.0
            if check_safe && is_safe_cut!(m, wp._objective_saf)
                ci_saf = MOI.add_constraint(m.relaxed_optimizer, wp._objective_saf, LT(UBD - wp._objective_saf.constant + constraint_tol))
                push!(m._objective_cut_ci_saf, ci_saf)
            end
            wp._objective_saf.constant = formulated_constant

        elseif obj_type === SCALAR_QUADRATIC
            buffered_sqf = wp._objective_sqf
            finite_cut_generated =  affine_relax_quadratic!(buffered_sqf.func, buffered_sqf.buffer,
                                                            buffered_sqf.saf, m._current_node, m._sol_to_branch_map,
                                                            m._current_xref)

            if finite_cut_generated
                if !check_safe || is_safe_cut!(m, buffered_sqf.saf)
                    copyto!(wp._objective_saf.terms, buffered_sqf.saf.terms)
                    wp._objective_saf.constant = 0.0
                    ci_saf = MOI.add_constraint(m.relaxed_optimizer, wp._objective_saf, LT(UBD - buffered_sqf.saf.constant + constraint_tol))
                    push!(m._objective_cut_ci_saf, ci_saf)
                end
            end

        elseif obj_type === NONLINEAR
            objective_cut_nonlinear!(m, wp, UBD, check_safe)
        end

        m._new_eval_objective = false
    end

    return nothing
end

"""
$(TYPEDSIGNATURES)

A routine that adds relaxations for all nonlinear constraints and quadratic constraints
corresponding to the current node to the relaxed problem. This adds an objective cut
(if specified by `objective_cut_on`) and then sets the `_new_eval_constraint` flag
to false indicating that an initial evaluation of the constraints has occurred. If
the `objective_cut_on` flag is `true` then the `_new_eval_objective` flag is also
set to `false` indicating that the objective expression was evaluated.
"""
function relax_all_constraints!(t::ExtensionType, m::Optimizer, q::Int64)

    check_safe = (q === 1) ? false : m._parameters.cut_safe_on
    m._working_problem._relaxed_evaluator.is_first_eval = m._new_eval_constraint

    sqf_leq_list = m._working_problem._sqf_leq
    for i = 1:m._working_problem._sqf_leq_count
        sqf_leq = @inbounds sqf_leq_list[i]
        relax!(m, sqf_leq, i, check_safe)
    end

    sqf_eq_list = m._working_problem._sqf_eq
    for i = 1:m._working_problem._sqf_eq_count
        sqf_eq = @inbounds sqf_eq_list[i]
        relax!(m, sqf_eq, i, check_safe)
    end

    nl_list = m._working_problem._nonlinear_constr
    for i = 1:m._working_problem._nonlinear_count
        nl = @inbounds nl_list[i]
        relax!(m, nl, i, check_safe)
    end

    m._new_eval_constraint = false

    objective_cut!(m, check_safe)

    return nothing
end
relax_constraints!(t::ExtensionType, m::Optimizer, q::Int64) = relax_all_constraints!(t, m, q)
relax_constraints!(m::Optimizer, q::Int64) = relax_constraints!(m.ext_type, m, q)

"""
$(FUNCTIONNAME)

Deletes all nonlinear constraints added to the relaxed optimizer.
"""
function delete_nl_constraints!(m::Optimizer)

    # delete affine relaxations added from quadratic inequality
    for ci in m._buffered_quadratic_ineq_ci
        MOI.delete(m.relaxed_optimizer, ci)
    end
    empty!(m._buffered_quadratic_ineq_ci)

    # delete affine relaxations added from quadratic equality
    for ci in m._buffered_quadratic_eq_ci
        MOI.delete(m.relaxed_optimizer, ci)
    end
    empty!(m._buffered_quadratic_eq_ci)

    # delete affine relaxations added from nonlinear inequality
    for ci in m._buffered_nonlinear_ci
        MOI.delete(m.relaxed_optimizer, ci)
    end
    empty!(m._buffered_nonlinear_ci)

    return nothing
end

"""
$(FUNCTIONNAME)

Deletes all scalar-affine objective cuts added to the relaxed optimizer.
"""
function delete_objective_cuts!(m::Optimizer)

    for ci in m._objective_cut_ci_saf
        MOI.delete(m.relaxed_optimizer, ci)
    end
    empty!(m._objective_cut_ci_saf)

    return nothing
end

"""
$(FUNCTIONNAME)

"""
function set_first_relax_point!(m::Optimizer)

    m._working_problem._relaxed_evaluator.is_first_eval = true
    m._new_eval_constraint = true
    m._new_eval_objective = true

    n = m._current_node
    @__dot__ m._current_xref = 0.5*(n.upper_variable_bounds + n.lower_variable_bounds)
    unsafe_check_fill!(isnan, m._current_xref, 0.0, length(m._current_xref))

    return nothing
end


text:
__precompile__()
module SeisXcorrelation

using SeisIO, Noise, Printf, Dates, FFTW, JLD2, MPI

export seisxcorrelation

"""

    seisxcorrelation(maxtimelag::Real, finame::String, foname::String; IsAllComponent::Bool=false)

Compute cross-correlation function and save data in jld2 file with SeisData format.

# Arguments
- `finame::String,`    : Input file name e.g. network = "BPnetwork"
- `maxtimelag::Real`    : Maximum lag time e.g. maxtimelag = 100 [s]
- `IsAllComponent::Bool`   : If true, compute 3 X 3 components cross-correlation

# Output
- `foname.jld2`                 : contains SeisData structure with a hierarchical structure (CC function, metadata)

"""
function seisxcorrelation(maxtimelag::Real, finame::String, foname::String; IsAllComponents::Bool=false) #Add whatever necessary arguments

    MPI.Init()
    # establish the MPI communicator and obtain rank
    comm = MPI.COMM_WORLD
    size = MPI.Comm_size(comm)
    rank = MPI.Comm_rank(comm)

    #Read validation data from ["../EXAMPLE/codevalidation"]

    #Make Station pairs (after validating CC part, read from actual dataset and use statino list in SeisData structure)
    #Take into account IsAllComponents == true or not

    mpiitrcount = 0
    baton = Array{Int32, 1}([0]) # for Relay Dumping algorithm

    for i = 1:size(data["info/timestamplist"])

        #parallelize one processor for one time stamp
        processID = stid - (size * mpiitrcount)

        # if this mpiitrcount is final round or not
        length(stlist) - size >= size * mpiitrcount ? anchor_rank = size-1 : anchor_rank = mod(length(stlist), size)-1

        if rank == processID-1

            #roop each station pair

            for j = 1:length("""stationpairlist`""")

                # compute crosscorrelation using Noise.jl

                # store in SeisChannel structure
                # time, fs, misc["stationpair"], misc["cc1"]

                # save dataset

                # Save data one by one
                if size == 1
                    save_SeisData2JLD2(foname, varname, S)

                else
                    if rank == 0
                        save_SeisData2JLD2(foname, varname, S)

                        if anchor_rank != 0
                            MPI.Send(baton, rank+1, 11, comm)
                            MPI.Recv!(baton, anchor_rank, 12, comm)
                        end

                    elseif rank == anchor_rank
                        MPI.Recv!(baton, rank-1, 11, comm)
                        save_SeisData2JLD2(foname, varname, S)
                        MPI.Send(baton, 0, 12, comm)

                    else
                        MPI.Recv!(baton, rank-1, 11, comm)
                        save_SeisData2JLD2(foname, varname, S)
                        MPI.Send(baton, rank+1, 11, comm)
                    end
                end
            end

            mpiitrcount += 1

        end
    end

    if rank == 0 println("Crosscorrelation and Saving data is successfully done.\njob ended at "*string(now())) end

    MPI.Finalize()

    return 0
end

end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
@testset "completions" begin
    using REPL.REPLCompletions: completions

    comps(line; mod = Main, force = false, lineNumber = 1, column = 1, text = "") =
        Atom.basecompletionadapter(line, mod, force, lineNumber, column, text)[1]

    @testset "module completion" begin
        ## basic
        let line = "@", cs = comps(line, mod = Atom)
            # test no error occurs in completion processing
            @test length(cs) == length(completions(line, lastindex(line), Atom)[1])
            # test detecting all the macros in Atom module
            @test filter(cs) do c
                c[:type] == "snippet" &&
                c[:icon] == "icon-mention" &&
                c[:rightLabel] == "Atom"
            end |> length == filter(names(Atom, all=true, imported=true)) do n
                startswith(string(n), '@')
            end |> length
        end

        ## advanced
        let cs = comps("match(code", mod = Atom)
            filter(cs) do c
                c[:text] == "codeblock_regex" &&
                c[:rightLabel] == "Atom"
            end |> !isempty
        end

        ## don't error for `ModuleCompletion` of modules where `@doc` isn't defined
        let line = "Core.Compiler.type"
            @test length(comps(line)) == length(completions(line, lastindex(line))[1])
        end
    end

    # method completion
    @testset "method completion" begin
        line = "push!("

        ## basic
        let cs = comps(line)
            # test no error occurs in completion processing
            @test length(cs) == length(completions(line, lastindex(line))[1])
            # test detecting all the `push!` methods available
            @test filter(cs) do c
                c[:type] == "method"
            end |> length == length(methods(push!))
        end

        ## features
        # - dynamic method addition
        # - shows module where the method defined in right label
        # - shows the infered return type in left label
        let cs = comps(line)
            @test filter(cs) do c
                c[:type] == "method" &&
                c[:rightLabel] == "Atom" &&
                c[:leftLabel] == "String"
            end |> isempty
        end

        @eval Atom begin
            import Base: push!
            push!(::Undefined) = "i'm such a dynamic push!"
        end

        let cs = comps(line)
            @test filter(cs) do c
                c[:type] == "method" &&
                c[:rightLabel] == "Atom" &&
                c[:leftLabel] == "String" &&
                c[:text] == "push!(::Atom.Undefined)"
            end |> !isempty
        end

        ## advanced
        line *= "Undefined(), "
        let cs = comps(line) # in Main: Undefined is not defined -- should show all push! methods
            @test length(cs) == length(methods(push!))
        end

        let c = comps(line, mod = Atom)[1] # in Atom: show push!(::Undefined) **first**
            @test c[:type] == "method"
            @test c[:rightLabel] == "Atom"
            @test c[:leftLabel] == "String"
            @test c[:text] == "push!(::Atom.Undefined)"
        end
    end

    @eval Main dict = Dict(:a => 1, :b => 2)

    @testset "property completion" begin
        @test map(comps("dict.")) do c
            c[:type] == "property" &&
            c[:rightLabel] == "Main" &&
            c[:text] ∈ string.(propertynames(dict))
        end |> all
    end

    @testset "field completion" begin
        @test map(comps("split(\"\", \"\")[1].")) do c
            c[:type] == "property" &&
            c[:rightLabel] == "SubString{String}" &&
            c[:leftLabel] ∈ string.(ntuple(i -> fieldtype(SubString{String}, i), fieldcount(SubString{String}))) &&
            c[:text] ∈ string.(fieldnames(SubString{String}))
        end |> all
    end

    @testset "dictionary completion" begin
        @test map(comps("dict[")) do c
            c[:type] == "property" &&
            c[:icon] == "icon-key" &&
            c[:rightLabel] == "Main" &&
            c[:leftLabel] == string(Int) &&
            c[:text] ∈ sprint.(show, keys(dict))
        end |> all
    end

    @testset "keyword completion" begin
        @test filter(comps("begin")) do c
            c[:type] == "keyword" &&
            c[:rightLabel] |> isempty &&
            c[:description] |> !isempty
        end |> !isempty
    end

    @testset "path completion" begin
        let line = "\"", cs = comps(line)
            @test length(cs) ==
                  length(completions(line, lastindex(line))[1]) ==
                  length(readdir(@__DIR__))
            @test filter(cs) do c
                c[:type] == "path" &&
                c[:icon] == "icon-file" &&
                c[:rightLabel] |> isempty
            end |> !isempty
        end
    end

    # completion suppressing
    @testset "suppressing: $(troublemaker)" for troublemaker ∈ [" ", "(", "[", "\$"]
        @test length(comps(troublemaker, force = false)) == 0 # surpressed
        @test length(comps(troublemaker, force = true)) > Atom.MAX_COMPLETIONS # invoked
    end
end

@testset "local completions" begin
    # local completions ordered by proximity
    let str = """
        function foo(x)
            aaa = 3
            a

            foo = x
            a
            abc = 3
            a
        end
        """
        @test Atom.basecompletionadapter("a", Main, true, 3, 6, str)[1][1][:text] == "aaa"
        @test Atom.basecompletionadapter("a", Main, true, 6, 6, str)[1][1][:text] == "aaa"
        @test Atom.basecompletionadapter("a", Main, true, 8, 6, str)[1][1][:text] == "abc"
    end

    # don't error on fallback case
    @test_nowarn @test Atom.localcompletions("", 1, 1) == []
end


text:
import TOML
using OrderedCollections
using Conda
using PyCall
import YAML
using Pkg.Artifacts
import GitHub
import GitHub: gh_get_json, DEFAULT_API
using HTTP
import JSON
using ghr_jll

include("generator/generator.jl")

sources = TOML.parsefile(joinpath(@__DIR__, "Sources.toml"))

build_dir = joinpath(@__DIR__, "build")

artifact_file = joinpath(@__DIR__, "..", "Artifacts.toml")

generate(ARGS, sources, build_dir, artifact_file)


text:
function shufflesim(decks::DecklistArray; shuffledecks = shuffle, parameters = HearthstoneSimParameters(),
    sets = sets, reps, batchsize)
   total_ouput = SimOutput(sets,  reps * batchsize,  length(decks))

   index = 1
   for i = 1:reps
       simres = hssim(batchsize, shuffledecks(decks), sets,deepcopy(parameters))
       index = copy_into_output!(total_ouput, simres; startindex = index)
   end

   total_ouput
end

function shuffle(x::DecklistArray)
    ans = deepcopy(x)
    sets = keys(x.commons)
    n = length(x.commons[first(sets)])
    p = randperm(n)
    for r in 1:4
        for set in sets
            permute!(ans[r][set],p)
        end
    end
    ans
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
include("step_script.jl")

# Simulate physics for 1/240 seconds (the default timestep)
function simstep(s::mjSim)
    # Create local simulator d (data), and m (model) variables
    d = s.d
    m = s.m

    if s.robot !== nothing
        # Check for gamepad input
        gamepad(s, GLFW.JOYSTICK_1)

        # Execute next step in command script
        step_script(s::mjSim, s.robot)

        # Step the controller forward by dt
        run!(s.robot)

        # Apply updated joint angles to sim
        d.ctrl .= unsafe_wrap(Array{Float64,1}, pointer(s.robot.state.joint_angles), 12)

        # If Pupper controller, subtract the l1 joint angles from the l2 joint angles
        # to fake the kinematics of the parallel linkage
        if true  # TODO: verify that the controller is a Pupper quadruped controller
            d.ctrl[[3,6,9,12]] .= d.ctrl[[3,6,9,12]] - d.ctrl[[2,5,8,11]]
        end
    end

    if s.paused
        if s.pert[].active > 0
            mjv_applyPerturbPose(m, d, s.pert, 1)  # move mocap and dynamic bodies
            mj_forward(m, d)
        end
    else
        #slow motion factor: 10x
        factor = s.slowmotion ? 10 : 1

        # advance effective simulation time by 1/refreshrate
        startsimtm = d.time
        starttm = time()
        refreshtm = 1.0/(factor*s.refreshrate)
        updates = refreshtm / m.opt.timestep

        steps = round(Int, round(s.framecount+updates)-s.framecount)
        s.framecount += updates

        for i=1:steps
            # clear old perturbations, apply new
            d.xfrc_applied .= 0.0
            if s.pert[].select > 0
                mjv_applyPerturbPose(m, d, s.pert, 0) # move mocap bodies only
                mjv_applyPerturbForce(m, d, s.pert)
            end

            mj_step(m, d)

            # break on reset
            (d.time < startsimtm) && break
        end
    end
end


text:
function build_cache!(cache::TimeStatusChange, op_model::OperationModel)
    build_cache!(cache, op_model.canonical)
end

function build_cache!(cache::TimeStatusChange, canonical::CanonicalModel)
    parameter = get_value(canonical, cache.ref)
    value_array = JuMP.Containers.DenseAxisArray{Dict{Symbol, Float64}}(undef, axes(parameter)...)

    for name in parameter.axes[1]
        status = PJ.value(parameter[name])
        value_array[name] = Dict(:count => 999.0, :status => status)
    end

    setfield!(cache, :value, value_array)

    return

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test

@testset "density_sample" begin

    param1 = [Float64(1.),-0.1,0.5]
    param4 = ones(Float64, 3)
    ds1 = @inferred DensitySample(param1, Float32(-3.3868156), 1)
    ds2 = similar(ds1)
    ds4 = @inferred DensitySample(param4, Float32(-4.2568156), 1)

    @testset "DensitySample" begin
        @test typeof(ds1)  <: AbstractDensitySample
        @test typeof(ds1)  <: DensitySample{Float64,Float32,Int}

        @test length(ds1) == 3
        @test nparams(ds1) == 3

        @test typeof(ds2)  <: DensitySample{Float64,Float32,Int}
        @test all(BAT.isoob.(ds2.params))
        @test isnan(ds2.log_value)
        @test iszero(ds2.weight)

        @test ds1 != ds4
        @test ds1 != ds2
        ds3 = @inferred DensitySample(param1, Float32(-3.3868156), 1)
        @test ds1 == ds3

        copy!(ds2, ds1)
        @test ds2 == ds3
    end

    @testset "DensitySampleVector" begin
        dsv1 = @inferred DensitySampleVector{Float64,Float32,Int}(3)
        @test typeof(dsv1) <: BAT.BATDataVector{DensitySample{Float64,Float32,Int,Vector{Float64}}}
        @test typeof(dsv1) <: DensitySampleVector{Float64,Float32,Int,ElasticArrays.ElasticArray{Float64,2,1},Vector{Float32},Vector{Int}}

        @test size(dsv1) == (0,)
        push!(dsv1, ds1)
        @test size(dsv1) == (1,)
        @test dsv1[1] == ds1
        @test IndexStyle(dsv1) == IndexLinear()

        push!(dsv1, ds2)
        dsv2 = @inferred DensitySampleVector{Float64, Float64, Float32}(3)
        push!(dsv2, ds2)
        push!(dsv2, ds4)
        append!(dsv1, dsv2)
        @test dsv1[4] == ds4
        @test dsv1[2] == ds1
    end
end


text:
sum(x, y) = x + y

@show sum(1, 1)

@show @code_lowered sum(1, 1)

@show @code_llvm sum(1, 1)

@show @code_native sum(1, 1)

max_stars_repo_name: dev-juliacn/JuliaTDD

max_issues_repo_name: dev-juliacn/JuliaTDD

max_forks_repo_name: dev-juliacn/JuliaTDD


text:
"""
    GeoAxis(args...; kwargs...) → ax
Create a new axis instance `ax` that is a modified `Axis` of the Makie.jl ecosystem.
All Makie.jl plotting functions work directly on `GeoAxis`, e.g., `scatter!(ax, x, y)`.

`GeoAxis` is appropriate for geospatial plotting because it automatically transforms all
plotted data given a user-defined map projection. See keyword arguments below and examples
in the online documentation. Longitude and latitude values in GeoMakie.jl are always
assumed to be **in degrees**.

In the call signature, `args...` is a standard figure location, e.g., `fig[1,1]` as given in
`Axis`. The keyword arguments decide the geospatial projection:

* `source = "+proj=longlat +datum=WGS84", dest = "+proj=eqearth"`: These two keywords
  configure the map projection to be used for the given field using Proj4.jl.
  See also online the section [Changing central longitude](@ref) for data that may not
  span the (expected by default) longitude range from -180 to 180.
* `transformation = Proj4.Transformation(source, dest, always_xy=true)`: Instead of
  `source, dest` you can directly use the Proj4.jl package to define the projection.
  lines!(ax, GeoMakie.coastlines(); coastkwargs...)
* `lonticks = -180:60:180, latticks = -90:30:90` ticks for the longitude and latitude
  dimensions. The grid lines of the axis are also spanning these tick values.
* `hidespines = true` Hide the axis spines (rectangle surrounding the axis).
* `coastlines = false` draw coastlines
* `coastline_attributes = (;)` named tuple that gets passed to the `lines` call drawing the coastline

## Example
```julia
using GeoMakie, GLMakie

lons = -180:180
lats = -90:90
field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

# Plot coastlines
coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, overdraw = true, coastkwargs...)
translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements

# Surface example
fig = Figure()
ax = GeoAxis(fig[1,1])
surface!(ax, lons, lats, field)
display(fig)

# Scatter example
slons = rand(lons, 2000)
slats = rand(lats, 2000)
sfield = [exp(cosd(l)) + 3(y/90) for (l,y) in zip(slons, slats)]

fig = Figure()
ax = GeoAxis(fig[1,1])
el = scatter!(slons, slats; color = sfield)
display(fig)
```
"""
function GeoAxis(args...;
        source = "+proj=longlat +datum=WGS84", dest = "+proj=eqearth",
        transformation = Proj4.Transformation(source, dest, always_xy=true),
        lonticks = -180:60:180,
        latticks = -90:30:90,
        hidespines = true,
        coastlines = false,
        coastline_attributes = (;),
        aspect = DataAspect(),
        kw...
    )
    ptrans = Makie.PointTrans{2}() do p
        transformation(p) ./ 10_0000f0
    end

    xticks = first.(Makie.apply_transform(ptrans, Point2f0.(lonticks, latticks[1])))
    yticks = last.(Makie.apply_transform(ptrans, Point2f0.(lonticks[1], latticks)))
    # Generate Axis instance
    ax = Axis(args...;
        aspect = aspect,
        xticks = (xticks, string.(lonticks, 'ᵒ')),
        yticks = (yticks, string.(latticks, 'ᵒ')),
        limits = (extrema(lonticks), extrema(latticks)),
        kw...)

    # Set axis transformation
    ax.scene.transformation.transform_func[] = ptrans

    # Plot coastlines
    if coastlines
        coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, coastline_attributes...)
        translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements
    end
    # Draw tick lines
    ax.xgridvisible=false; ax.ygridvisible=false
    # TODO: How to get "default" grid line style from the theme?
    for lon in lonticks
        coords = [Point2f(lon, l) for l in range(latticks[1], latticks[end]; length = 100)]
        gridplot = lines!(ax, coords; color = :gray20, linewidth = 0.5)
        translate!(gridplot, 0, 0, 100) # ensure they are on top of other plotted elements
    end

    for lat in latticks
        coords = [Point2f(l, lat) for l in range(lonticks[1], lonticks[end]; length = 100)]
        gridplot = lines!(ax, coords; color = :gray20, linewidth = 0.5)
        translate!(gridplot, 0, 0, 100) # ensure they are on top of other plotted elements
    end

    hidespines && hidespines!(ax)
    return ax
end


text:
module NumericalFluxes
using StaticArrays

"""
    rusanov!(F::MArray, nM, QM, QVM, auxM, QP, QVP, auxP, t, flux!, wavespeed,
             [computeQjump!])

Calculate the Rusanov (aka local Lax-Friedrichs) numerical flux given the plus
and minus side states/viscous states `QP`/`QVP` and `QM`/`QVM` using the physical
flux function `flux!` and `wavespeed` calculation.

The `flux!` has almost the same calling convention as `flux!` from
[`DGBalanceLaw`](@ref).

The function `wavespeed` should return the maximum wavespeed for a state and is
called as `wavespeed(nM, QM, auxM, t)` and `wavespeed(nM, QP, auxP, t)` where
`nM` is the outward unit normal for the minus side.

When present `computeQjump!(ΔQ, QM, auxM, QP, auxP)` will be called after so
that the user specify the value to use for `QM - QP`; this is useful for
correcting `Q` to include discontinuous reference states.

"""
function rusanov!(F::MArray{Tuple{nstate}}, nM,
                  QM, QVM, auxM,
                  QP, QVP, auxP,
                  t, flux!, wavespeed,
                  computeQjump! = nothing
                 ) where {nstate}
  λM = wavespeed(nM, QM, auxM, t)
  FM = similar(F, Size(3, nstate))
  flux!(FM, QM, QVM, auxM, t)

  λP = wavespeed(nM, QP, auxP, t)
  FP = similar(F, Size(3, nstate))
  flux!(FP, QP, QVP, auxP, t)

  λ  =  max(λM, λP)

  if computeQjump! === nothing
    @inbounds for s = 1:nstate
      F[s] = (nM[1] * (FM[1, s] + FP[1, s]) + nM[2] * (FM[2, s] + FP[2, s]) +
              nM[3] * (FM[3, s] + FP[3, s]) + λ * (QM[s] - QP[s])) / 2
    end
  else
    ΔQ = copy(QM)
    computeQjump!(ΔQ, QM, auxM, QP, auxP)
    @inbounds for s = 1:nstate
      F[s] = (nM[1] * (FM[1, s] + FP[1, s]) + nM[2] * (FM[2, s] + FP[2, s]) +
              nM[3] * (FM[3, s] + FP[3, s]) + λ * ΔQ[s]) / 2
    end
  end
end

"""
    rusanov_boundary_flux!(F::MArray{Tuple{nstate}}, nM, QM, QVM, auxM, QP, QVP,
                           auxP, bctype, t, flux!, bcstate!, wavespeed,
                           computeQjump! = nothing) where {nstate}

The function `bcstate!` is used to calculate the plus side state for the
boundary condition `bctype`. The calling convention is:
```
bcstate!(QP, QVP, auxP, nM, QM, QVM, auxM, bctype, t)
```
where `QP`, `QVP`, and `auxP` are the plus side state, viscous state, and
auxiliary state to be filled from the given data; other arguments should not be
modified.
"""
function rusanov_boundary_flux!(F::MArray{Tuple{nstate}}, nM,
                                QM, QVM, auxM,
                                QP, QVP, auxP,
                                bctype, t,
                                flux!, bcstate!,
                                wavespeed,
                                computeQjump! = nothing
                               ) where {nstate}
  bcstate!(QP, QVP, auxP, nM, QM, QVM, auxM, bctype, t)
  rusanov!(F, nM, QM, QVM, auxM, QP, QVP, auxP, t, flux!, wavespeed,
           computeQjump!)
end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
import Theta: parity_char, remainder_char

@testset "Characteristics" begin
    @test parity_char([[0,0],[0,1]]) == 0
    @test parity_char([[1,1],[0,1]]) == 1
    @test remainder_char([[0,1,2],[3,1,5]]) == [[0,1,0],[1,1,1]]
    @test Set(theta_char(1)) == Set([[[0],[0]], [[0],[1]], [[1],[0]], [[1],[1]]])
    @test Set(even_theta_char(1)) == Set([[[0],[0]], [[0],[1]], [[1],[0]]])
    @test odd_theta_char(1) == [[[1],[1]]]
    @test check_azygetic([[[1,0,1,0], [1,0,1,0]], [[0,0,0,1], [1,0,0,0]], [[0,0,1,1], [1,0,1,1]]]) == true
end


text:
@defcomp impactextratropicalstorms begin
    regions = Index()

    extratropicalstormsdam = Variable(index=[time,regions])
    extratropicalstormsdead = Variable(index=[time,regions])

    extratropicalstormsbasedam = Parameter(index=[regions])
    extratropicalstormsdamel = Parameter(default = -0.514)
    extratropicalstormspar = Parameter(index=[regions])
    extratropicalstormsbasedead = Parameter(index=[regions])
    extratropicalstormsdeadel = Parameter(default = -0.501)
    extratropicalstormsnl = Parameter(default = 1)

    gdp90 = Parameter(index=[regions])
    pop90 = Parameter(index=[regions])

    population = Parameter(index=[time,regions])
    income = Parameter(index=[time,regions])

    acco2 = Parameter(index=[time])
    co2pre = Parameter()

    function run_timestep(p, v, d, t)

        for r in d.regions
            ypc = p.income[t, r] / p.population[t, r] * 1000.0
            ypc90 = p.gdp90[r] / p.pop90[r] * 1000.0

            v.extratropicalstormsdam[t, r] = p.extratropicalstormsbasedam[r] * p.income[t, r] * (ypc / ypc90)^p.extratropicalstormsdamel * ((1.0 + (p.extratropicalstormspar[r] * (p.acco2[t] / p.co2pre)))^p.extratropicalstormsnl - 1.0)
            v.extratropicalstormsdead[t, r] = 1000.0 * p.extratropicalstormsbasedead[r] * p.population[t, r] * (ypc / ypc90)^p.extratropicalstormsdeadel * ((1.0 + (p.extratropicalstormspar[r] * (p.acco2[t] / p.co2pre)))^p.extratropicalstormsnl - 1.0)
        end
    end
end


text: (x + 1) y

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia