text:
"""
  msm_optimize!(sMMProblem::MSMProblem; verbose::Bool = true)

Function to launch an optimization. To be used after the following functions
have been called: (i) set_empirical_moments! (ii) set_priors!
(iii) set_simulate_empirical_moments! (iv) construct_objective_function!
"""
function msm_optimize!(sMMProblem::MSMProblem; verbose::Bool = true)

  # Initialize a BlackBoxOptim problem
  # this modifies sMMProblem.bbSetup
  #-----------------------------------
  set_global_optimizer!(sMMProblem)

  # If the global optimizer is using BlackBoxOptim
  #-----------------------------------------------
  if is_bb_optimizer(sMMProblem.options.globalOptimizer) == true

    # Store best fitness and best candidates
    listBestFitness = []
    listBestCandidates = []

    # Run the optimization with BlackBoxOptim
    #----------------------------------------
    sMMProblem.bbResults = bboptimize(sMMProblem.bbSetup)

    push!(listBestFitness, best_fitness(sMMProblem.bbResults))
    push!(listBestCandidates, best_candidate(sMMProblem.bbResults))

  # In the future, we may use other global minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.globalOptimizer = $(sMMProblem.options.globalOptimizer) is not supported.")

  end

  return listBestFitness, listBestCandidates

end

"""
  function msm_minimizer(sMMProblem::MSMProblem)

Function to get the parameter value minimizing the objective function
"""
function msm_minimizer(sMMProblem::MSMProblem)

  # If the global optimizer is using BlackBoxOptim
  #-----------------------------------------------
  if is_bb_optimizer(sMMProblem.options.globalOptimizer) == true

    best_candidate(sMMProblem.bbResults)

  # In the future, we may use other global minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.globalOptimizer = $(sMMProblem.options.globalOptimizer) is not supported.")

  end

end

"""
  function msm_minimum(sMMProblem::MSMProblem)

Function to get the minimum of the objective function
"""
function msm_minimum(sMMProblem::MSMProblem)

  # If the global optimizer is using BlackBoxOptim
  #-----------------------------------------------
  if is_bb_optimizer(sMMProblem.options.globalOptimizer) == true

    best_fitness(sMMProblem.bbResults)

  # In the future, we may use other global minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.globalOptimizer = $(sMMProblem.options.globalOptimizer) is not supported.")

  end

end

"""
  msm_refine_globalmin!(sMMProblem::MSMProblem; verbose::Bool = true)

Function to refine the global minimum using a local minimization routine.
To be used after the following functions have been called: (i) set_empirical_moments!
(ii) set_priors! (iii) set_simulate_empirical_moments! (iv) construct_objective_function!
(v) msm_optimize!
"""
function msm_refine_globalmin!(sMMProblem::MSMProblem; verbose::Bool = true)

  x0 = msm_minimizer(sMMProblem)

  # Let's use the result from the global maximizer as the starting value
  #---------------------------------------------------------------------
  if verbose == true
    info("Refining the global maximum using a local algorithm.")
    info("Using Fminbox = $(sMMProblem.options.minBox)")
    info("Starting value = $(x0)")
  end

  if is_local_optimizer(sMMProblem.options.localOptimizer) == true

    # If using Fminbox option is true
    #-------------------------------
    if sMMProblem.options.minBox == true

        lower = create_lower_bound(sMMProblem)
        upper = create_upper_bound(sMMProblem)

        sMMProblem.optimResults = optimize(sMMProblem.objective_function, lower, upper, x0, convert_to_fminbox(sMMProblem.options.localOptimizer), Optim.Options(iterations = sMMProblem.options.maxFuncEvals))

    else
        sMMProblem.optimResults = optimize(sMMProblem.objective_function, x0, convert_to_optim_algo(sMMProblem.options.localOptimizer), Optim.Options(iterations = sMMProblem.options.maxFuncEvals))
    end

  # In the future, we may use other local minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.localOptimizer = $(sMMProblem.options.localOptimizer) is not supported.")

  end

  return Optim.minimizer(sMMProblem.optimResults)

end

"""
  function msm_local_minimizer(sMMProblem::MSMProblem)

Function to get the parameter value minimizing the objective function (local)
"""
function msm_local_minimizer(sMMProblem::MSMProblem)

  # If the global optimizer is using BlackBoxOptim
  #-----------------------------------------------
  if is_optim_optimizer(sMMProblem.options.localOptimizer) == true

    Optim.minimizer(sMMProblem.optimResults)

  # In the future, we may use other global minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.localOptimizer = $(sMMProblem.options.localOptimizer) is not supported.")

  end

end

"""
  function msm_local_minimum(sMMProblem::MSMProblem)

Function to get the local minimum value of the objetive function
"""
function msm_local_minimum(sMMProblem::MSMProblem)

  # If the global optimizer is using BlackBoxOptim
  #-----------------------------------------------
  if is_optim_optimizer(sMMProblem.options.localOptimizer) == true

    Optim.minimum(sMMProblem.optimResults)

  # In the future, we may use other global minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.localOptimizer = $(sMMProblem.options.localOptimizer) is not supported.")

  end

end

"""
  function msm_multistart_minimizer(sMMProblem::MSMProblem)

Function to get the parameter value minimizing the objective function when
using the multistart algorithm
"""
function msm_multistart_minimizer(sMMProblem::MSMProblem)

  # Result given by msm_local_minimizer.
  msm_local_minimizer(sMMProblem::MSMProblem)

end

"""
  function msm_multistart_minimum(sMMProblem::MSMProblem)

Function to get the minimum value of the objetive function when
using the multistart algorithm
"""
function msm_multistart_minimum(sMMProblem::MSMProblem)

  # Result given by msm_local_minimum
  msm_local_minimum(sMMProblem::MSMProblem)

end

"""
  msm_multistart!(sMMProblem::MSMProblem; x0 = Array{Float64}(undef, 0,0), nums::Int64 = nworkers(), verbose::Bool = true)

Function to run several local minimization algorithms in parallel, with different
starting values. The minimum is calculated as the minimum of the local minima.
Changes sMMProblem.optimResults. This function also returns
a list containing Optim results.
"""
function msm_multistart!(sMMProblem::MSMProblem; x0 = Array{Float64}(undef, 0,0), nums::Int64 = nworkers(), verbose::Bool = true)

  # Safety checks
  #--------------
  if nums < nworkers()
    errors("nums < nworkers()")
  elseif nums > nworkers()
    info("nums > nworkers(). Some starting values will be ignored.")
  end

  # To store minization results
  #----------------------------
  results = []

  # Look for valid starting values (for which convergence is reached)
  #-------------------------------------------------------------------
  if x0 == Array{Float64}(undef, 0,0)
    myGrid = search_starting_values(sMMProblem, nums, verbose = verbose)
  # Using starting values provided by the user
  #-------------------------------------------
  else

    # Check that enough starting values were provided
    if size(x0, 1) < nworkers()
      info("$(size(x0, 1)) starting value(s) were provided.")
      error("The minimum number of starting value(s) to provide is $(nworkers()).")
    end

    myGrid = x0

  end

  # If the local optimizer is using Optim
  #--------------------------------------
  if is_optim_optimizer(sMMProblem.options.localOptimizer) == true

      # A. Starting tasks on available workers
      #---------------------------------------
      @sync for (workerIndex, w) in enumerate(workers())

        @async push!(results, @fetchfrom w wrap_msm_localmin(sMMProblem, myGrid[workerIndex,:], verbose = true))

      end

    # B. Looking for the minimum
    #----------------------------
    # Initialization
    minIndex = 0
    minValue = Inf
    minimizerValue = zeros(length(keys(sMMProblem.priors)))
    nbConvergenceReached = 0
    listOptimResults = []

    for (workerIndex, w) in enumerate(workers())

      push!(listOptimResults, results[workerIndex])

      try

        minimumValue = Optim.minimum(results[workerIndex])
        minimizer = Optim.minimizer(results[workerIndex])

        if minimumValue < minValue && Optim.converged(results[workerIndex]) == true

          minIndex = workerIndex
          minValue = minimumValue
          minimizerValue = minimizer
          nbConvergenceReached += 1

        end

      catch myError
        info("$(myError)")
      end

    end

    # D. If none of the optimization converged
    #-----------------------------------------
    if minIndex == 0
      info("None of the local optimizer algorithm converged.")
    else
      info("Convergence reached for $(nbConvergenceReached) worker(s).")
      info("Minimum value found with worker $(minIndex)")
      sMMProblem.optimResults = results[minIndex]
    end

  # In the future, we may use other local minimizer
  # routines. For the moment, let's return an error
  #-------------------------------------------------
  else

    error("sMMProblem.options.localOptimizer = $(sMMProblem.options.localOptimizer) is not supported.")

  end

  if verbose == true
    if nbConvergenceReached != 0
      info("Best value found with starting values = $(myGrid[minIndex,:]).")
      info("Best value = $(minValue).")
      info("Minimizer = $(minimizerValue)")
    end
  end

  return listOptimResults

end

"""
  msm_localmin(sMMProblem::MSMProblem, x0::Array{Float64,1}; verbose::Bool = true)

Function find a local minimum using a local minimization routine, with starting value x0.
To be used after the following functions have been called: (i) set_empirical_moments!
(ii) set_priors! (iii) set_simulate_empirical_moments! (iv) construct_objective_function!
"""
function msm_localmin(sMMProblem::MSMProblem, x0::Array{Float64,1}; verbose::Bool = true)

    # Let's use the result from the global maximizer as the starting value
    #---------------------------------------------------------------------
    if verbose == true
        info("Starting value = $(x0)")
        info("Using Fminbox = $(sMMProblem.options.minBox)")
    end

    if is_local_optimizer(sMMProblem.options.localOptimizer) == true

    # If using Fminbox option is true
    #-------------------------------
    if sMMProblem.options.minBox == true

        lower = create_lower_bound(sMMProblem)
        upper = create_upper_bound(sMMProblem)

        optimResults = optimize(sMMProblem.objective_function, lower, upper, x0, convert_to_fminbox(sMMProblem.options.localOptimizer), Optim.Options(iterations = sMMProblem.options.maxFuncEvals))

    else
        optimResults = optimize(sMMProblem.objective_function, x0, convert_to_optim_algo(sMMProblem.options.localOptimizer), Optim.Options(iterations = sMMProblem.options.maxFuncEvals))
    end

    # In the future, we may use other local minimizer
    # routines. For the moment, let's return an error
    #-------------------------------------------------
    else

    error("sMMProblem.options.localOptimizer = $(sMMProblem.options.localOptimizer) is not supported.")

    end

    return optimResults

end

"""

"""
function wrap_msm_localmin(sMMProblem::MSMProblem, x0::Array{Float64,1}; verbose::Bool = true)

  try
    msm_localmin(sMMProblem, x0, verbose = verbose)
  catch myError
    info("$(myError)")
    info("Error with msm_localmin")
  end

end

"""
  search_starting_values(sMMProblem::MSMProblem, numPoints::Int64; verbose::Bool = true)

Search for nums valid starting values. To be used after the following functions have been called:
(i) set_empirical_moments! (ii) set_priors! (iii) set_simulate_empirical_moments!
(iv) construct_objective_function!
"""
function search_starting_values(sMMProblem::MSMProblem, numPoints::Int64; verbose::Bool = true)

  # Safety Check
  #-------------
  if is_optim_optimizer(sMMProblem.options.localOptimizer) == false
    error("sMMProblem.options.localOptimizer = $(sMMProblem.options.localOptimizer) is not supported.")
  end

  if verbose == true
    info("Searching for $(numPoints) valid starting value(s)")
  end

  # Generate upper and lower bounds vector
  #--------------------------------------------------------------------------
  lower_bound = zeros(length(keys(sMMProblem.priors)))
  upper_bound = zeros(length(keys(sMMProblem.priors)))

  for (kIndex, k) in enumerate(keys(sMMProblem.priors))
    lower_bound[kIndex] = sMMProblem.priors[k][2]
    upper_bound[kIndex] = sMMProblem.priors[k][3]
  end

  #Each row is a new point and each column is a dimension of this points.
  #---------------------------------------------------------------------
  Validx0 = zeros(numPoints, length(lower_bound))
  distanceValue = zeros(numPoints) #to store the distance associated to each point
  nbValidx0Found = 0

  # Create many grids (stochastic draws) with many potential points
  #----------------------------------------------------------------
  if verbose == true
    info("Creating $(sMMProblem.options.maxTrialsStartingValues) potential starting value(s)")
    info("gridType = $(sMMProblem.options.gridType)")
  end

  # Generate many points for the grid
  #-----------------------------------------------------------------------------
  if sMMProblem.options.gridType == :LHC
    candidates_starting_values = latin_hypercube_sampling(generate_bbSearchRange(sMMProblem), Int(sMMProblem.options.maxTrialsStartingValues*numPoints))
  elseif sMMProblem.options.gridType == :Sobol
    candidates_starting_values = sobol_sampling(lower_bound, upper_bound, Int(sMMProblem.options.maxTrialsStartingValues*numPoints))
  else
    err("sMMProblem.options.gridType = $(sMMProblem.options.gridType) is not a valid sampling procedure.")
  end

  # Split the grid into chunks
  #-----------------------------------------------------------------------------
  listGrids = []
  i = 1;
  j = i + numPoints - 1;
  for k=1:sMMProblem.options.maxTrialsStartingValues
      push!(listGrids, candidates_starting_values[i:j,:])
      i = i + numPoints;
      j = i + numPoints - 1;
  end

  listGridsIndex = 0

  # Looping until numPoints valid points have been found
  #----------------------------------------------------------------------------
  while nbValidx0Found < numPoints

    results = []
    listGridsIndex += 1

    if listGridsIndex > sMMProblem.options.maxTrialsStartingValues
      error("Maximum number of attempts reached without success. maxTrialsStartingValues = $(sMMProblem.options.maxTrialsStartingValues)")
    end

    # Use available workers to simulate moments
    #------------------------------------------
    @sync for (workerIndex, w) in enumerate(workers())

      @async push!(results, @fetchfrom w sMMProblem.objective_function(listGrids[listGridsIndex][workerIndex, :]))

    end

    # Check for convergence
    #----------------------
    for (workerIndex, w) in enumerate(workers())

      # Set penalty value by default
      distanceValue[workerIndex] = sMMProblem.options.penaltyValue

      try

        distanceValue[workerIndex] = results[workerIndex]

        # discard inf distances, values equal to penaltyValue and values above the threshold
        if isinf(distanceValue[workerIndex]) == false && distanceValue[workerIndex] != sMMProblem.options.penaltyValue && distanceValue[workerIndex] < sMMProblem.options.thresholdStartingValue

          nbValidx0Found +=1

          if nbValidx0Found <= numPoints

            Validx0[nbValidx0Found,:] = listGrids[listGridsIndex][workerIndex, :]
            info("Valid starting value = $(Validx0[nbValidx0Found,:]), distance = $(distanceValue[workerIndex])")
          end

        end

      catch myError
        info("$(myError)")
      end

    end

  end

  # sorting starting values according to distance value (in ascending order)
  p = sortperm(distanceValue) #get the ascending order
  Validx0 = Validx0[p,:]  #re-order rows

  if verbose == true
    info("Found $(nbValidx0Found) valid starting value(s)")
  end

  # If requested, save (valid) starting values generated
  if sMMProblem.options.saveStartingValues == true
    if verbose == true
      info("Saving starting values to disk.")
    end

    tempfilename = "starting_values_"* sMMProblem.options.saveName * ".bson"
    bson(tempfilename, Dict(:Validx0=>Validx0))

    tempfilename = "starting_distances_"* sMMProblem.options.saveName * ".bson"
    bson(tempfilename, Dict(:distanceValue=>distanceValue))

  end

  return Validx0

end


text:
export OutdatedWarning

module OutdatedWarning
using Gumbo, AbstractTrees, Documenter

OLD_VERSION_CSS = replace("""
.outdated-warning-overlay {
  position: fixed;
  top: 0;
  left: 0;
  right: 0;
  box-shadow: 0 0 10px rgba(0, 0, 0, 0.3);
  z-index: 999;
  background-color: #ffaba7;
  color: rgba(0, 0, 0, 0.7);
  border-bottom: 3px solid #da0b00;
  padding: 10px 35px;
  text-align: center;
  font-size: 15px; }
  .outdated-warning-overlay .outdated-warning-closer {
    position: absolute;
    top: calc(50% - 10px);
    right: 18px;
    cursor: pointer;
    width: 12px; }
  .outdated-warning-overlay a {
    color: #2e63b8; }
    .outdated-warning-overlay a:hover {
      color: #363636; }
""", '\n' => "")

OUTDATED_VERSION_ATTR = isdefined(Documenter.Writers.HTMLWriter, :OUTDATED_VERSION_ATTR) ?
    Documenter.Writers.HTMLWriter.OUTDATED_VERSION_ATTR : "data-outdated-warner"

OLD_VERSION_WARNER = """
function maybeAddWarning () {
    const head = document.getElementsByTagName('head')[0];

    // Add a noindex meta tag (unless one exists) so that search engines don't index this version of the docs.
    if (document.body.querySelector('meta[name="robots"]') === null) {
        const meta = document.createElement('meta');
        meta.name = 'robots';
        meta.content = 'noindex';

        head.appendChild(meta);
    };

    // Add a stylesheet to avoid inline styling
    const style = document.createElement('style');
    style.type = 'text/css';
    style.appendChild(document.createTextNode('$(OLD_VERSION_CSS)'));
    head.appendChild(style);

    const div = document.createElement('div');
    div.classList.add('outdated-warning-overlay');
    const closer = document.createElement('div');
    closer.classList.add('outdated-warning-closer');

    // Icon by font-awesome (license: https://fontawesome.com/license, link: https://fontawesome.com/icons/times?style=solid)
    closer.innerHTML = '<svg aria-hidden="true" focusable="false" data-prefix="fas" data-icon="times" class="svg-inline--fa fa-times fa-w-11" role="img" xmlns="http://www.w3.org/2000/svg" viewBox="0 0 352 512"><path fill="currentColor" d="M242.72 256l100.07-100.07c12.28-12.28 12.28-32.19 0-44.48l-22.24-22.24c-12.28-12.28-32.19-12.28-44.48 0L176 189.28 75.93 89.21c-12.28-12.28-32.19-12.28-44.48 0L9.21 111.45c-12.28 12.28-12.28 32.19 0 44.48L109.28 256 9.21 356.07c-12.28 12.28-12.28 32.19 0 44.48l22.24 22.24c12.28 12.28 32.2 12.28 44.48 0L176 322.72l100.07 100.07c12.28 12.28 32.2 12.28 44.48 0l22.24-22.24c12.28-12.28 12.28-32.19 0-44.48L242.72 256z"></path></svg>';
    closer.addEventListener('click', function () {
        document.body.removeChild(div);
    });
    let href = '/stable';
    if (window.documenterBaseURL) {
        href = window.documenterBaseURL + '/../stable';
    }
    div.innerHTML = 'This documentation is not for the latest version.
 Go to the latest documentation.';
    div.appendChild(closer);
    document.body.appendChild(div);
};

if (document.readyState === 'loading') {
    document.addEventListener('DOMContentLoaded', maybeAddWarning);
} else {
    maybeAddWarning();
};
"""

function add_old_docs_notice(index, force = false)
    html = read(index, String)
    parsed = Gumbo.parsehtml(html)

    did_update = false
    for el in PreOrderDFS(parsed.root)
        if el isa HTMLElement && Gumbo.tag(el) == :head
            old_notice = get_notice(el)

            if old_notice === nothing
                push!(el.children, make_notice(el))
                did_update = true
            elseif force
                update_notice(old_notice)
                did_update = true
            end

            break
        end
    end

    if did_update
        open(index, "w") do io
            print(io, parsed)
        end
    end
    return did_update
end

function make_notice(parent)
    script = Gumbo.HTMLElement{:script}([], parent, Dict(OUTDATED_VERSION_ATTR => ""))
    content = Gumbo.HTMLText(script, OLD_VERSION_WARNER)
    push!(script.children, content)
    return script
end

function update_notice(script)
    content = Gumbo.HTMLText(script, OLD_VERSION_WARNER)
    empty!(script.children)
    push!(script.children, content)
    return script
end

function get_notice(html)
    for el in PreOrderDFS(html)
        if el isa HTMLElement && Gumbo.tag(el) == :script
            attrs = Gumbo.attrs(el)
            if haskey(attrs, OUTDATED_VERSION_ATTR)
                return el
            end
        end
    end
    return nothing
end

"""
    generate([io::IO = stdout,] root::String;force = false)

This function adds a warning (and `noindex` meta tag) to all versions of
the documentation in `root`.

`force` overwrites a previous injected warning message created by this function.

A typical use case is to run this on the `gh-pages` branch of a packge.
"""
generate(root::String; kwargs...) = generate(stdout, root; kwargs...)
function generate(io::IO, root::String;force = false)
    for dir in readdir(root)
        path = joinpath(root, dir)
        islink(path) && continue
        isdir(path) || continue
        index = joinpath(path, "index.html")
        isfile(index) || continue

        if endswith(path, "dev")
            println(io, "Skipping $(dir) since it's a dev version.")
            continue
        end

        print(io, "Processing $(dir): ")
        for (root, _, files) in walkdir(path)
            for file in files
                _, ext = splitext(file)
                if ext == ".html"
                    try
                        did_change = add_old_docs_notice(joinpath(root, file), force)
                        print(io, did_change ? "✓" : ".")
                    catch err
                        if err isa InterruptException
                            rethrow()
                        end
                        @debug "Fatally failed to add a outdated warning" exception = (err, catch_backtrace())
                        print(io, "!")
                    end
                end
            end
        end
        println(io)
    end
end

end

ext: jl


text:
module TestBench

using Reactive, ReactiveBasics, BenchmarkTools

const N = 10^3
const V = rand(N)
function testReactive()
    a = Reactive.Signal(0.0)
    b = Reactive.Signal(1.0)
    c = map(+, a, b)
    d = Reactive.foldp(+, 0.0, a)
    x = 0.0
    for i in 1:N
        push!(a, V[i])
        Reactive.run(1)
        x = Reactive.value(c) + Reactive.value(d)
    end
    x
end

function testReactiveBasics()
    a = ReactiveBasics.Signal(0.0)
    b = ReactiveBasics.Signal(1.0)
    c = map(+, a, b)
    d = ReactiveBasics.foldp(+, 0.0, a)
    x = 0.0
    for i in 1:N
        push!(a, V[i])
        x = ReactiveBasics.value(c) + ReactiveBasics.value(d)
    end
    x
end

@show @benchmark testReactive()
@show @benchmark testReactiveBasics()

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""NamespaceStatus is information about the current status of a Namespace.

    IoK8sApiCoreV1NamespaceStatus(;
        conditions=nothing,
        phase=nothing,
    )

    - conditions::Vector{IoK8sApiCoreV1NamespaceCondition} : Represents the latest available observations of a namespace's current state.
    - phase::String : Phase is the current lifecycle phase of the namespace. More info: https://kubernetes.io/docs/tasks/administer-cluster/namespaces/
"""
mutable struct IoK8sApiCoreV1NamespaceStatus <: SwaggerModel
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1NamespaceCondition} } # spec name: conditions
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase

    function IoK8sApiCoreV1NamespaceStatus(;conditions=nothing, phase=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NamespaceStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1NamespaceStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        o
    end
end # type IoK8sApiCoreV1NamespaceStatus

const _property_map_IoK8sApiCoreV1NamespaceStatus = Dict{Symbol,Symbol}(Symbol("conditions")=>Symbol("conditions"), Symbol("phase")=>Symbol("phase"))
const _property_types_IoK8sApiCoreV1NamespaceStatus = Dict{Symbol,String}(Symbol("conditions")=>"Vector{IoK8sApiCoreV1NamespaceCondition}", Symbol("phase")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1NamespaceStatus }) = collect(keys(_property_map_IoK8sApiCoreV1NamespaceStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1NamespaceStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NamespaceStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NamespaceStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NamespaceStatus[property_name]

function check_required(o::IoK8sApiCoreV1NamespaceStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NamespaceStatus }, name::Symbol, val)
end


text:
module CP_MC

include("../Misc.jl")
include("BackwardKernel.jl")
using .Misc, .BackwardKernel

using Distributions, SpecialFunctions

export Particle, smc_weight, particle_mcmc

mutable struct Particle{ℜ <: Real}
    S_R::Array{ℜ,1}
    S_IR::Array{ℜ,2}
    S_JR::Array{ℜ,2}
    S_KR::Array{ℜ,2}
    function Particle(I::Int,J::Int, K::Int, R::Int; γ::ℜ=0.1) where {ℜ<:Real}
        a::ℜ = I*J*K*γ
        return new{ℜ}(fill(a/R,R), fill(a/(R*I),I,R), fill(a/(R*J),J,R), fill(a/(R*K),K,R))
    end
    function Particle(S_R::Array{ℜ,1}, S_IR::Array{ℜ,2}, S_JR::Array{ℜ,2}, S_KR::Array{ℜ,2}) where {ℜ<:Real}
        return new{ℜ}(S_R, S_IR, S_JR, S_KR)
    end
end

function smc_weight(X::Array{ℜ,3}, R::Ƶ, N::Ƶ=1; μ::ℜ=3.5, γ::ℜ=0.1, adaptive::Bool=true) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ, K::Ƶ = size(X)
    T::Ƶ = Ƶ(sum(X))
    ESS::Array{ℜ} = zeros(ℜ,T)
    ESS_min::Array{ℜ} = adaptive ? fill(N/2.0,T) : fill(N+1.0,T)

    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    log_Zᵗ::ℜ = a*log(b) - (a+T)*log(b + 1.0) - sum(lgamma.(X .+ 1.0))
    P = [Particle(I,J,K,R; γ=γ) for n=1:N]
    P_temp = [Particle(I,J,K,R; γ=γ) for n=1:N]
    log_wᵗ::Array{ℜ}, Wᵗ::Array{ℜ}, log_Wᵗ::Array{ℜ}, cum_Wᵗ::Array{ℜ} = fill(log_Zᵗ,N), fill(1.0/N,N), fill(-log(N),N), zeros(ℜ,N)

    log_νᵗ::ℜ, log_q_r::Array{ℜ}, q_r::Array{ℜ}= 0.0, zeros(ℜ,R), zeros(ℜ,R)
    rᵗ::Ƶ, iᵗ::Ƶ, jᵗ::Ƶ, kᵗ::Ƶ, uᵗ::ℜ  = 0, 0, 0, 0, 0.0

    for (t,(iᵗ, jᵗ, kᵗ)) in enumerate(EventQueue(X))
        for (n,p) in enumerate(P)

            log_q_r .= log.(p.S_IR[iᵗ,:]) .+ log.(p.S_JR[jᵗ,:]) .+ log.(p.S_KR[kᵗ,:]) .- 2.0.*log.(p.S_R)
            log_νᵗ = logsumexp(log_q_r)
            log_q_r .-= log_νᵗ
            q_r .= exp.(log_q_r)

            rᵗ = rand(Categorical(q_r))
            p.S_R[rᵗ] += 1.0
            p.S_IR[iᵗ,rᵗ] += 1.0
            p.S_JR[jᵗ,rᵗ] += 1.0
            p.S_KR[kᵗ,rᵗ] += 1.0

            log_wᵗ[n] += log_νᵗ
        end

        log_Zᵗ = logsumexp(log_wᵗ)
        log_Wᵗ .= log_wᵗ .- log_Zᵗ
        log_Zᵗ -= log(N)
        Wᵗ .= exp.(log_Wᵗ)

        ESS[t] = 1.0/sum(Wᵗ .* Wᵗ)

        if ESS[t] < ESS_min[t]
            for n=1:N
                P_temp[n].S_R .= P[n].S_R
                P_temp[n].S_IR .= P[n].S_IR
                P_temp[n].S_JR .= P[n].S_JR
                P_temp[n].S_KR .= P[n].S_KR
            end

            cum_Wᵗ .= cumsum(Wᵗ)
            uᵗ = rand()/N
            p_id = 1

            for n=1:N # systematic resampling
                while cum_Wᵗ[p_id] < uᵗ
                    p_id += 1
                end
                P[n].S_R .= P_temp[p_id].S_R
                P[n].S_IR .= P_temp[p_id].S_IR
                P[n].S_JR .= P_temp[p_id].S_JR
                P[n].S_KR .= P_temp[p_id].S_KR
                uᵗ += 1.0/N
            end
            log_wᵗ .= log_Zᵗ
        end
    end
    p_id = rand(Categorical(Wᵗ))
    P[p_id].S_R .-= a/R
    P[p_id].S_IR .-= a/(R*I)
    P[p_id].S_JR .-= a/(R*J)
    P[p_id].S_KR .-= a/(R*K)
    return log_Zᵗ, ESS, P[p_id]
end

function particle_mcmc(X::Array{ℜ,3}, R::Ƶ, N::Ƶ=1; μ::ℜ=3.5, γ::ℜ=0.1, EPOCHS::Ƶ=1, adaptive=false) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ, K::Ƶ = size(X)
    a::ℜ, b::ℜ = I*J*K*γ, γ/μ

    log_PX::Array{ℜ}, stays::Ƶ = zeros(ℜ, EPOCHS), 0
    log_PX[1], ESS::Array{ℜ}, part::Particle{ℜ} = smc_weight(X, R, N; μ=μ, γ=γ, adaptive=adaptive)
    log_PX_prop::ℜ, ESS_prop::Array{ℜ}, part_prop::Particle{ℜ} = smc_weight(X, R, N; μ=μ, γ=γ, adaptive=adaptive)

    for eph=2:EPOCHS
        log_PX_prop, ESS_prop, part_prop = smc_weight(X, R, N; μ=μ, γ=γ, adaptive=adaptive)

        if log(rand()) < log_PX_prop - log_PX[eph-1]
            log_PX[eph] = log_PX_prop
            ESS .= ESS_prop
            part.S_R .= part_prop.S_R
            part.S_IR .= part_prop.S_IR
            part.S_JR .= part_prop.S_JR
            part.S_KR .= part_prop.S_KR
        else
            log_PX[eph] = log_PX[eph-1]
            stays += 1
        end
    end

    log_Z_τ::ℜ = logsumexp(log_PX) - log(EPOCHS)
    return log_Z_τ, log_PX, 1.0 - stays./EPOCHS, part
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
#####
##### Timescale for diffusion across one cell
#####

using Oceananigans.Grids: topology

function min_Δxyz(grid)
    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx, topo[1])
    Δy = minimum_grid_spacing(grid.Δy, topo[2])
    Δz = minimum_grid_spacing(grid.Δz, topo[3])

    return min(Δx, Δy, Δy)
end

cell_diffusion_timescale(model) = cell_diffusion_timescale(model.closure, model.diffusivities, model.grid)
cell_diffusion_timescale(::Nothing, diffusivities, grid) = Inf

maximum_numeric_diffusivity(κ::Number) = κ
maximum_numeric_diffusivity(κ::NamedTuple) = maximum(κ)
maximum_numeric_diffusivity(κ::NamedTuple{()}) = 0 # tracers=nothing means empty diffusivity tuples

# As the name suggests, we give up in the case of a function diffusivity
maximum_numeric_diffusivity(κ::Function) = 0

minimum_grid_spacing(Δx, TX)           = minimum(Δx)
minimum_grid_spacing(Δx, ::Type{Flat}) = Inf

function cell_diffusion_timescale(closure::IsotropicDiffusivity, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_κ = maximum_numeric_diffusivity(closure.κ)
    max_ν = maximum_numeric_diffusivity(closure.ν)
    return min(Δ^2 / max_ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::AnisotropicDiffusivity, diffusivities, grid)
    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx, topo[1])
    Δy = minimum_grid_spacing(grid.Δy, topo[2])
    Δz = minimum_grid_spacing(grid.Δz, topo[3])

    max_νx = maximum_numeric_diffusivity(closure.νx)
    max_νy = maximum_numeric_diffusivity(closure.νy)
    max_νz = maximum_numeric_diffusivity(closure.νz)

    max_κx = maximum_numeric_diffusivity(closure.κx)
    max_κy = maximum_numeric_diffusivity(closure.κy)
    max_κz = maximum_numeric_diffusivity(closure.κz)

    return min(Δx^2 / max_νx,
               Δy^2 / max_νy,
               Δz^2 / max_νz,
               Δx^2 / max_κx,
               Δy^2 / max_κy,
               Δz^2 / max_κz)
end

function cell_diffusion_timescale(closure::AnisotropicBiharmonicDiffusivity, diffusivities, grid)
    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx, topo[1])
    Δy = minimum_grid_spacing(grid.Δy, topo[2])
    Δz = minimum_grid_spacing(grid.Δz, topo[3])

    max_νx = maximum_numeric_diffusivity(closure.νx)
    max_νy = maximum_numeric_diffusivity(closure.νy)
    max_νz = maximum_numeric_diffusivity(closure.νz)

    max_κx = maximum_numeric_diffusivity(closure.κx)
    max_κy = maximum_numeric_diffusivity(closure.κy)
    max_κz = maximum_numeric_diffusivity(closure.κz)

    return min(Δx^4 / max_νx,
               Δy^4 / max_νy,
               Δz^4 / max_νz,
               Δx^4 / max_κx,
               Δy^4 / max_κy,
               Δz^4 / max_κz)
end

function cell_diffusion_timescale(closure::HorizontallyCurvilinearAnisotropicDiffusivity, diffusivities, grid)
    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx, topo[1])
    Δy = minimum_grid_spacing(grid.Δy, topo[2])
    Δz = minimum_grid_spacing(grid.Δz, topo[3])

    max_νh = maximum_numeric_diffusivity(closure.νh)
    max_νz = maximum_numeric_diffusivity(closure.νz)

    max_κh = maximum_numeric_diffusivity(closure.κh)
    max_κz = maximum_numeric_diffusivity(closure.κz)

    return min(Δx^2 / max_νh,
               Δy^2 / max_νh,
               Δz^2 / max_νz,
               Δx^2 / max_κh,
               Δy^2 / max_κh,
               Δz^2 / max_κz)
end

function cell_diffusion_timescale(closure::HorizontallyCurvilinearAnisotropicBiharmonicDiffusivity, diffusivities, grid)
    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx, topo[1])
    Δy = minimum_grid_spacing(grid.Δy, topo[2])
    Δz = minimum_grid_spacing(grid.Δz, topo[3])

    max_νh = maximum_numeric_diffusivity(closure.νh)
    max_νz = maximum_numeric_diffusivity(closure.νz)

    max_κh = maximum_numeric_diffusivity(closure.κh)
    max_κz = maximum_numeric_diffusivity(closure.κz)

    return min(Δx^4 / max_νh,
               Δy^4 / max_νh,
               Δz^4 / max_νz,
               Δx^4 / max_κh,
               Δy^4 / max_κh,
               Δz^4 / max_κz)
end

function cell_diffusion_timescale(closure::SmagorinskyLilly{FT, P, <:NamedTuple{()}},
                                  diffusivities, grid) where {FT, P}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::SmagorinskyLilly, diffusivities, grid)
    Δ = min_Δxyz(grid)
    min_Pr = minimum(closure.Pr)
    max_κ = maximum(closure.κ)
    max_νκ = maximum(diffusivities.νₑ.data.parent) * max(1, 1/min_Pr)
    return min(Δ^2 / max_νκ, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::AnisotropicMinimumDissipation{FT, PK, PN, <:NamedTuple{()}},
                                  diffusivities, grid) where {FT, PK, PN}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::AnisotropicMinimumDissipation, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    max_κ = max(Tuple(maximum(κₑ.data.parent) for κₑ in diffusivities.κₑ)...)
    return min(Δ^2 / max_ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::TwoDimensionalLeith, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

cell_diffusion_timescale(closure::Tuple, diffusivities, grid) =
    min(Tuple(cell_diffusion_timescale(c, diffusivities, grid) for c in closure)...)


text:
# Relational OverApproximation Utilities

mutable struct OverApproximation
    output::Union{Symbol, Real}
    output_range::Array{T, 1} where {T <: Real}
    ranges::Dict{Symbol, Array{T, 1}} where {T <: Real}
    nvars::Integer
    consts::Array{Symbol, 1}
    approx_eq::Array{Expr, 1}
    approx_ineq::Array{Expr, 1}
    fun_eq::Dict{Symbol, Any}
    N::Integer # number of points in a bound in a [0,1] interval
    ϵ::Real # construct overapprox with ϵ leeway from function
end
# default constructor
OverApproximation() = OverApproximation(:(null_output),
                                        Array{Float64, 1}(),
                                        Dict{Symbol, Array{Float64,1}}(),
                                        0,
                                        Array{Symbol, 1}(),
                                        Array{Expr, 1}(),
                                        Array{Expr, 1}(),
                                        Dict{Symbol, Any}(),
                                        3,
                                        1e-2)

Base.print(oA::OverApproximation) = print_overapproximate(oA)
Base.display(oA::OverApproximation) = print_overapproximate(oA)
function print_overapproximate(oA::OverApproximation)
    println("output = $(oA.output)")
    for eq in oA.approx_eq
        println(eq)
    end
    for ineq in oA.approx_ineq
        println(ineq)
    end
    println(oA.output_range)
end

"""
This function combines a list of OverApproximation objects and return one
    that contains all their equality and ineqaulity constraints.

    ** for output, output_range, nvars, const, fun_eq, N and ϵ, we keep the default values
"""
function add_overapproximate(list_oA::Array{OverApproximation, 1})

    out_oA = OverApproximation()
    out_oA.ranges = merge([oA.ranges for oA in list_oA]...)
    out_oA.approx_eq = vcat([oA.approx_eq for oA in list_oA]...)
    out_oA.approx_ineq = vcat([oA.approx_ineq for oA in list_oA]...)
    return out_oA
end


text:
@doc raw"""
    CompressibleEulerEquations2D

The compressible Euler equations for an ideal gas in two space dimensions.
"""
struct CompressibleEulerEquations2D <: AbstractCompressibleEulerEquations{2, 4}
  gamma::Float64
end

function CompressibleEulerEquations2D()
  gamma = parameter("gamma", 1.4)

  CompressibleEulerEquations2D(gamma)
end

get_name(::CompressibleEulerEquations2D) = "CompressibleEulerEquations2D"
varnames_cons(::CompressibleEulerEquations2D) = @SVector ["rho", "rho_v1", "rho_v2", "rho_e"]
varnames_prim(::CompressibleEulerEquations2D) = @SVector ["rho", "v1", "v2", "p"]

# Set initial conditions at physical location `x` for time `t`
function initial_conditions_density_pulse(x, t, equation::CompressibleEulerEquations2D)
  rho = 1 + exp(-(x[1]^2 + x[2]^2))/2
  v1 = 1
  v2 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 1
  rho_e = p/(equation.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return @SVector [rho, rho_v1, rho_v2, rho_e]
end

"""
    initial_conditions_density_wave(x, t, equation::CompressibleEulerEquations2D)

test case for stability of EC fluxes from paper: https://arxiv.org/pdf/2007.09026.pdf
domain [-1, 1]^2
mesh = 4x4, polydeg = 5
"""
function initial_conditions_density_wave(x, t, equation::CompressibleEulerEquations2D)
  v1 = 0.1
  v2 = 0.2
  rho = 1 + 0.98 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2)))
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 20
  rho_e = p / (equation.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return @SVector [rho, rho_v1, rho_v2, rho_e]
end

function initial_conditions_pressure_pulse(x, t, equation::CompressibleEulerEquations2D)
  rho = 1
  v1 = 1
  v2 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 1 + exp(-(x[1]^2 + x[2]^2))/2
  rho_e = p/(equation.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return @SVector [rho, rho_v1, rho_v2, rho_e]
end

function initial_conditions_density_pressure_pulse(x, t, equation::CompressibleEulerEquations2D)
  rho = 1 + exp(-(x[1]^2 + x[2]^2))/2
  v1 = 1
  v2 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  p = 1 + exp(-(x[1]^2 + x[2]^2))/2
  rho_e = p/(equation.gamma - 1) + 1/2 * rho * (v1^2 + v2^2)
  return @SVector [rho, rho_v1, rho_v2, rho_e]
end

function initial_conditions_constant(x, t, equation::CompressibleEulerEquations2D)
  rho = 1.0
  rho_v1 = 0.1
  rho_v2 = -0.2
  rho_e = 10.0
  return @SVector [rho, rho_v1, rho_v2, rho_e]
end

function initial_conditions_convergence_test(x, t, equation::CompressibleEulerEquations2D)
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  ini = c + A * sin(ω * (x[1] + x[2] - t))

  rho = ini
  rho_v1 = ini
  rho_v2 = ini
  rho_e = ini^2

  return @SVector [rho, rho_v1, rho_v2, rho_e]
end

function initial_conditions_isentropic_vortex(x, t, equation::CompressibleEulerEquations2D)
  # needs appropriate mesh size, e.g. [-10,-10]x[10,10]
  # make sure that the inicenter does not exit the domain, e.g. T=10.0
  # initial center of the vortex
  inicenter = SVector(0.0, 0.0)
  # size and strength of the vortex
  iniamplitude = 0.2
  # base flow
  rho = 1.0
  v1 = 1.0
  v2 = 1.0
  vel = SVector(v1, v2)
  p = 10.0
  vec = SVector(v1, v2)
  rt = p / rho                  # ideal gas equation
  cent = inicenter + vel*t      # advection of center
  cent = x - cent               # distance to centerpoint
  #cent=cross(iniaxis,cent)     # distance to axis, tangent vector, length r
  # cross product with iniaxis = [0,0,1]
  cent = SVector(-cent[2], cent[1])
  r2 = cent[1]^2 + cent[2]^2
  du = iniamplitude/(2*π)*exp(0.5*(1-r2)) # vel. perturbation
  dtemp = -(equation.gamma-1)/(2*equation.gamma*rt)*du^2            # isentrop
  rho = rho * (1+dtemp)^(1\(equation.gamma-1))
  vel = vel + du*cent
  v1, v2 = vel
  p = p * (1+dtemp)^(equation.gamma/(equation.gamma-1))
  prim = SVector(rho, v1, v2, p)
  return prim2cons(prim, equation)
end

function initial_conditions_weak_blast_wave(x, t, equation::CompressibleEulerEquations2D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up polar coordinates
  inicenter = SVector(0.0, 0.0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)
  phi = atan(y_norm, x_norm)
  sin_phi, cos_phi = sincos(phi)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1  = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  v2  = r > 0.5 ? 0.0 : 0.1882 * sin_phi
  p   = r > 0.5 ? 1.0 : 1.245

  return prim2cons(SVector(rho, v1, v2, p), equation)
end

function initial_conditions_blast_wave(x, t, equation::CompressibleEulerEquations2D)
  # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
  # Set up polar coordinates
  inicenter = SVector(0.0, 0.0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)
  phi = atan(y_norm, x_norm)
  sin_phi, cos_phi = sincos(phi)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1  = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  v2  = r > 0.5 ? 0.0 : 0.1882 * sin_phi
  p   = r > 0.5 ? 1.0E-3 : 1.245

  return prim2cons(SVector(rho, v1, v2, p), equation)
end

function initial_conditions_sedov_blast_wave(x, t, equation::CompressibleEulerEquations2D)
  # Set up polar coordinates
  inicenter = SVector(0.0, 0.0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  # r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 3 * (equation.gamma - 1) * E / (3 * pi * r0^2)
  p0_outer = 1.0e-5 # = true Sedov setup
  # p0_outer = 1.0e-3 # = more reasonable setup

  # Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  v2  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, v2, p), equation)
end

function initial_conditions_medium_sedov_blast_wave(x, t, equation::CompressibleEulerEquations2D)
  # Set up polar coordinates
  inicenter = SVector(0.0, 0.0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  # r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 3 * (equation.gamma - 1) * E / (3 * pi * r0^2)
  # p0_outer = 1.0e-5 # = true Sedov setup
  p0_outer = 1.0e-3 # = more reasonable setup

  # Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  v2  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, v2, p), equation)
end

function initial_conditions_khi(x, t, equation::CompressibleEulerEquations2D)
  # https://rsaa.anu.edu.au/research/established-projects/fyris/2-d-kelvin-helmholtz-test
  # change discontinuity to tanh
  # typical resolution 128^2, 256^2
  # domain size is [-0.5,0.5]^2
  dens0 = 1.0 # outside density
  dens1 = 2.0 # inside density
  velx0 = -0.5 # outside velocity
  velx1 = 0.5 # inside velocity
  slope = 50 # used for tanh instead of discontinuous initial condition
  # pressure equilibrium
  p     = 2.5
  #  y velocity v2 is only white noise
  v2  = 0.01*(rand(Float64,1)[1]-0.5)
  # density
  rho = dens0 + (dens1-dens0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))
  #  x velocity is also augmented with noise
  v1 = velx0 + (velx1-velx0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))+0.01*(rand(Float64,1)[1]-0.5)
  return prim2cons(SVector(rho, v1, v2, p), equation)
end

function initial_conditions_blob(x, t, equation::CompressibleEulerEquations2D)
  # blob test case, see Agertz et al. https://arxiv.org/pdf/astro-ph/0610051.pdf
  # other reference: https://arxiv.org/pdf/astro-ph/0610051.pdf
  # change discontinuity to tanh
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
  # resolution 128^2, 256^2
  # domain size is [-20.0,20.0]^2
  # gamma = 5/3 for this test case
  R = 1.0 # radius of the blob
  # background density
  dens0 = 1.0
  Chi = 10.0 # density contrast
  # reference time of characteristic growth of KH instability equal to 1.0
  tau_kh = 1.0
  tau_cr = tau_kh/1.6 # crushing time
  # determine background velocity
  velx0 = 2*R*sqrt(Chi)/tau_cr
  vely0 = 0.0
  Ma0 = 2.7 # background flow Mach number Ma=v/c
  c = velx0/Ma0 # sound speed
  # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
  p0 = c*c*dens0/equation.gamma
  # initial center of the blob
  inicenter = [-15,0]
  x_rel = x-inicenter
  r = sqrt(x_rel[1]^2 + x_rel[2]^2)
  # steepness of the tanh transition zone
  slope = 2
  # density blob
  dens = dens0 + (Chi-1) * 0.5*(1+(tanh(slope*(r+R)) - (tanh(slope*(r-R)) + 1)))
  # velocity blob is zero
  velx = velx0 - velx0 * 0.5*(1+(tanh(slope*(r+R)) - (tanh(slope*(r-R)) + 1)))
  return prim2cons(SVector(dens, velx, vely0, p0), equation)
end

function initial_conditions_jeans_instability(x, t, equation::CompressibleEulerEquations2D)
  # Jeans gravitational instability test case
  # see Derigs et al. https://arxiv.org/abs/1605.03572; Sec. 4.6
  # OBS! this uses cgs (centimeter, gram, second) units
  # periodic boundaries
  # domain size [0,L]^2 depends on the wave number chosen for the perturbation
  # OBS! Be very careful here L must be chosen such that problem is periodic
  # typical final time is T = 5
  # gamma = 5/3
  dens0  = 1.5e7 # g/cm^3
  pres0  = 1.5e7 # dyn/cm^2
  delta0 = 1e-3
  # set wave vector values for pertubation (units 1/cm)
  # see FLASH manual: https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf
  kx = 2.0*pi/0.5 # 2π/λ_x, λ_x = 0.5
  ky = 0.0   # 2π/λ_y, λ_y = 1e10
  k_dot_x = kx*x[1] + ky*x[2]
  # perturb density and pressure away from reference states ρ_0 and p_0
  dens = dens0*(1.0 + delta0*cos(k_dot_x))                # g/cm^3
  pres = pres0*(1.0 + equation.gamma*delta0*cos(k_dot_x)) # dyn/cm^2
  # flow starts as stationary
  velx = 0.0 # cm/s
  vely = 0.0 # cm/s
  return prim2cons(SVector(dens, velx, vely, pres), equation)
end

function initial_conditions_eoc_test_coupled_euler_gravity(x, t, equation::CompressibleEulerEquations2D)
  # OBS! this assumes that γ = 2 other manufactured source terms are incorrect
  if equation.gamma != 2.0
    error("adiabatic constant must be 2 for the coupling convergence test")
  end
  c = 2.0
  A = 0.1
  ini = c + A * sin(pi * (x[1] + x[2] - t))
  G = 1.0 # gravitational constant

  rho = ini
  v1 = 1.0
  v2 = 1.0
  p = ini^2 * G / pi # * 2 / ndims, but ndims==2 here

  return prim2cons(SVector(rho, v1, v2, p), equation)
end

function initial_conditions_sedov_self_gravity(x, t, equation::CompressibleEulerEquations2D)
  # Set up polar coordinates
  r = sqrt(x[1]^2 + x[2]^2)

  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000
  r0 = 0.125 # = 4.0 * smallest dx (for domain length=8 and max-ref=8)
  E = 1.0
  p_inner   = (equation.gamma - 1) * E / (pi * r0^2)
  p_ambient = 1e-5 # = true Sedov setup

  # Calculate primitive variables
  # use a logistic function to tranfer density value smoothly
  L  = 1.0    # maximum of function
  x0 = 1.0    # center point of function
  k  = -150.0 # sharpness of transfer
  logistic_function_rho = L/(1.0 + exp(-k*(r - x0)))
  rho_ambient = 1e-5
  rho = max(logistic_function_rho, rho_ambient) # clip background density to not be so tiny

  # velocities are zero
  v1 = 0.0
  v2 = 0.0

  # use a logistic function to tranfer pressure value smoothly
  logistic_function_p = p_inner/(1.0 + exp(-k*(r - r0)))
  p = max(logistic_function_p, p_ambient)

  return prim2cons(SVector(rho, v1, v2, p), equation)
end

# Apply boundary conditions
function boundary_conditions_convergence_test(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equation::CompressibleEulerEquations2D)
  u_boundary = initial_conditions_convergence_test(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

function boundary_conditions_sedov_self_gravity(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equation::CompressibleEulerEquations2D)
  # velocities are zero, density/pressure are ambient values according to
  # initial_conditions_sedov_self_gravity
  rho = 1e-5
  v1 = 0.0
  v2 = 0.0
  p = 1e-5

  u_boundary = prim2cons(SVector(rho, v1, v2, p), equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

# Apply source terms
function source_terms_convergence_test(ut, u, x, element_id, t, n_nodes, equation::CompressibleEulerEquations2D)
  # Same settings as in `initial_conditions`
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  γ = equation.gamma

  for j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, element_id]
    x2 = x[2, i, j, element_id]
    si, co = sincos((x1 + x2 - t)*ω)
    tmp1 = co * A * ω
    tmp2 = si * A
    tmp3 = γ - 1
    tmp4 = (2*c - 1)*tmp3
    tmp5 = (2*tmp2*γ - 2*tmp2 + tmp4 + 1)*tmp1
    tmp6 = tmp2 + c

    ut[1, i, j, element_id] += tmp1
    ut[2, i, j, element_id] += tmp5
    ut[3, i, j, element_id] += tmp5
    ut[4, i, j, element_id] += 2*((tmp6 - 1)*tmp3 + tmp6*γ)*tmp1

    # Original terms (without performanc enhancements)
    # ut[1, i, j, element_id] += cos((x1 + x2 - t)*ω)*A*ω
    # ut[2, i, j, element_id] += (2*sin((x1 + x2 - t)*ω)*A*γ - 2*sin((x1 + x2 - t)*ω)*A +
    #                             2*c*γ - 2*c - γ + 2)*cos((x1 + x2 - t)*ω)*A*ω
    # ut[3, i, j, element_id] += (2*sin((x1 + x2 - t)*ω)*A*γ - 2*sin((x1 + x2 - t)*ω)*A +
    #                             2*c*γ - 2*c - γ + 2)*cos((x1 + x2 - t)*ω)*A*ω
    # ut[4, i, j, element_id] += 2*((c - 1 + sin((x1 + x2 - t)*ω)*A)*(γ - 1) +
    #                               (sin((x1 + x2 - t)*ω)*A + c)*γ)*cos((x1 + x2 - t)*ω)*A*ω
  end

  return nothing
end

function source_terms_eoc_test_coupled_euler_gravity(ut, u, x, element_id, t, n_nodes, equation::CompressibleEulerEquations2D)
  # Same settings as in `initial_conditions_eoc_test_coupled_euler_gravity`
  c = 2.0
  A = 0.1
  G = 1.0 # gravitational constant, must match coupling solver
  C_grav = -2.0 * G / pi # 2 == 4 / ndims

  for j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, element_id]
    x2 = x[2, i, j, element_id]
    si, co = sincos(pi * (x1 + x2 - t))
    rhox = A * pi * co
    rho  = c + A *  si

    ut[1, i, j, element_id] += rhox
    ut[2, i, j, element_id] += rhox
    ut[3, i, j, element_id] += rhox
    ut[4, i, j, element_id] += (1.0 - C_grav*rho)*rhox
  end

  return nothing
end

function source_terms_eoc_test_euler(ut, u, x, element_id, t, n_nodes, equation::CompressibleEulerEquations2D)
  # Same settings as in `initial_conditions_eoc_test_coupled_euler_gravity`
  c = 2.0
  A = 0.1
  G = 1.0
  C_grav = -2 * G / pi # 2 == 4 / ndims

  for j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, element_id]
    x2 = x[2, i, j, element_id]
    si, co = sincos(pi * (x1 + x2 - t))
    rhox = A * pi * co
    rho  = c + A *  si

    ut[1, i, j, element_id] += rhox
    ut[2, i, j, element_id] += rhox * (1 -     C_grav * rho)
    ut[3, i, j, element_id] += rhox * (1 -     C_grav * rho)
    ut[4, i, j, element_id] += rhox * (1 - 3 * C_grav * rho)
  end

  return nothing
end

# Empty source terms required for coupled Euler-gravity simulations
function source_terms_harmonic(ut, u, x, element_id, t, n_nodes, equation::CompressibleEulerEquations2D)
  # OBS! used for the Jeans instability as well as self-gravitating Sedov blast
  # TODO: make this cleaner and let each solver have a different source term name
  return nothing
end

# Calculate 1D flux for a single point
@inline function calcflux(u, orientation, equation::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u
  v1 = rho_v1/rho
  v2 = rho_v2/rho
  p = (equation.gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))
  if orientation == 1
    f1 = rho_v1
    f2 = rho_v1 * v1 + p
    f3 = rho_v1 * v2
    f4 = (rho_e + p) * v1
  else
    f1 = rho_v2
    f2 = rho_v2 * v1
    f3 = rho_v2 * v2 + p
    f4 = (rho_e + p) * v2
  end
  return SVector(f1, f2, f3, f4)
end

"""
    function flux_shima_etal(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

This flux is is a modification of the original kinetic energy preserving two-point flux by
Kuya, Totani and Kawai (2018)
  Kinetic energy and entropy preserving schemes for compressible flows
  by split convective forms
  [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058)
The modification is in the energy flux to guarantee pressure equilibrium and was developed by
Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020)
  Preventing spurious pressure oscillations in split convective form discretizations for
  compressible flows
"""
@inline function flux_shima_etal(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  p_ll  = (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))
  p_rr  = (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg  = 1/2 * ( v1_ll +  v1_rr)
  v2_avg  = 1/2 * ( v2_ll +  v2_rr)
  p_avg   = 1/2 * (  p_ll +   p_rr)
  kin_avg = 1/2 * (v1_ll*v1_rr + v2_ll*v2_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    pv1_avg = 1/2 * (p_ll*v1_rr + p_rr*v1_ll)
    f1 = rho_avg * v1_avg
    f2 = rho_avg * v1_avg * v1_avg + p_avg
    f3 = rho_avg * v1_avg * v2_avg
    f4 = p_avg*v1_avg/(equation.gamma-1) + rho_avg*v1_avg*kin_avg + pv1_avg
  else
    pv2_avg = 1/2 * (p_ll*v2_rr + p_rr*v2_ll)
    f1 = rho_avg * v2_avg
    f2 = rho_avg * v2_avg * v1_avg
    f3 = rho_avg * v2_avg * v2_avg + p_avg
    f4 = p_avg*v2_avg/(equation.gamma-1) + rho_avg*v2_avg*kin_avg + pv2_avg
  end

  return SVector(f1, f2, f3, f4)
end

"""
    flux_kennedy_gruber(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

Kinetic energy preserving two-point flux by Kennedy and Gruber (2008)
  Reduced aliasing formulations of the convective terms within the
  Navier-Stokes equations for a compressible fluid
[DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020)
"""
@inline function flux_kennedy_gruber(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr

  # Average each factor of products in flux
  rho_avg = 1/2 * (rho_ll + rho_rr)
  v1_avg = 1/2 * (v1_ll + v1_rr)
  v2_avg = 1/2 * (v2_ll + v2_rr)
  p_avg = 1/2 * ((equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2)) +
                 (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2)))
  e_avg = 1/2 * (rho_e_ll/rho_ll + rho_e_rr/rho_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho_avg * v1_avg
    f2 = rho_avg * v1_avg * v1_avg + p_avg
    f3 = rho_avg * v1_avg * v2_avg
    f4 = (rho_avg * e_avg + p_avg) * v1_avg
  else
    f1 = rho_avg * v2_avg
    f2 = rho_avg * v2_avg * v1_avg
    f3 = rho_avg * v2_avg * v2_avg + p_avg
    f4 = (rho_avg * e_avg + p_avg) * v2_avg
  end

  return SVector(f1, f2, f3, f4)
end

"""
    flux_chandrashekar(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

Entropy conserving two-point flux by Chandrashekar (2013)
  Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes
  for Compressible Euler and Navier-Stokes Equations
[DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a)
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  p_ll =  (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))
  p_rr =  (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))
  beta_ll = 0.5*rho_ll/p_ll
  beta_rr = 0.5*rho_rr/p_rr
  specific_kin_ll = 0.5*(v1_ll^2 + v2_ll^2)
  specific_kin_rr = 0.5*(v1_rr^2 + v2_rr^2)

  # Compute the necessary mean values
  rho_avg  = 0.5*(rho_ll+rho_rr)
  rho_mean = ln_mean(rho_ll,rho_rr)
  beta_mean = ln_mean(beta_ll,beta_rr)
  beta_avg = 0.5*(beta_ll+beta_rr)
  v1_avg = 0.5*(v1_ll+v1_rr)
  v2_avg = 0.5*(v2_ll+v2_rr)
  p_mean = 0.5*rho_avg/beta_avg
  velocity_square_avg = specific_kin_ll + specific_kin_rr

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho_mean * v1_avg
    f2 = f1 * v1_avg + p_mean
    f3 = f1 * v2_avg
    f4 = f1 * 0.5*(1/(equation.gamma-1)/beta_mean - velocity_square_avg)+f2*v1_avg + f3*v2_avg
  else
    f1 = rho_mean * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_mean
    f4 = f1 * 0.5*(1/(equation.gamma-1)/beta_mean - velocity_square_avg)+f2*v1_avg + f3*v2_avg
  end

  return SVector(f1, f2, f3, f4)
end

"""
    flux_ranocha(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

Entropy conserving and kinetic energy preserving two-point flux by Ranocha (2018)
  Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods
  for Hyperbolic Balance Laws
[PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743)
See also Ranocha (2020)
  Entropy Conserving and Kinetic Energy Preserving Numerical Methods for
  the Euler Equations Using Summation-by-Parts Operators
[Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42)
"""
@inline function flux_ranocha(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)
  # Unpack left and right state
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  p_ll =  (equation.gamma - 1) * (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2))
  p_rr =  (equation.gamma - 1) * (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2))

  # Compute the necessary mean values
  rho_mean   = ln_mean(rho_ll, rho_rr)
  rho_p_mean = ln_mean(rho_ll / p_ll, rho_rr / p_rr)
  v1_avg = 0.5 * (v1_ll + v1_rr)
  v2_avg = 0.5 * (v2_ll + v2_rr)
  p_avg  = 0.5 * (p_ll + p_rr)
  velocity_square_avg = 0.5 * (v1_ll*v1_rr + v2_ll*v2_rr)

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = rho_mean * v1_avg
    f2 = f1 * v1_avg + p_avg
    f3 = f1 * v2_avg
    f4 = f1 * ( velocity_square_avg + 1 / ((equation.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)
  else
    f1 = rho_mean * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_avg
    f4 = f1 * ( velocity_square_avg + 1 / ((equation.gamma-1) * rho_p_mean) ) + 0.5 * (p_ll*v2_rr + p_rr*v2_ll)
  end

  return SVector(f1, f2, f3, f4)
end

function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v_mag_ll = sqrt(v1_ll^2 + v2_ll^2)
  p_ll = (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * v_mag_ll^2)
  c_ll = sqrt(equation.gamma * p_ll / rho_ll)
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v_mag_rr = sqrt(v1_rr^2 + v2_rr^2)
  p_rr = (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * v_mag_rr^2)
  c_rr = sqrt(equation.gamma * p_rr / rho_rr)

  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
  f1 = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (rho_rr    - rho_ll)
  f2 = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (rho_v1_rr - rho_v1_ll)
  f3 = 1/2 * (f_ll[3] + f_rr[3]) - 1/2 * λ_max * (rho_v2_rr - rho_v2_ll)
  f4 = 1/2 * (f_ll[4] + f_rr[4]) - 1/2 * λ_max * (rho_e_rr  - rho_e_ll)

  return SVector(f1, f2, f3, f4)
end

function flux_hll(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)
  # Calculate primitive variables and speed of sound
  rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  p_ll = (equation.gamma - 1) * (rho_e_ll - 1/2 * rho_ll * (v1_ll^2 + v2_ll^2))

  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  p_rr = (equation.gamma - 1) * (rho_e_rr - 1/2 * rho_rr * (v1_rr^2 + v2_rr^2))

  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  if orientation == 1 # x-direction
    Ssl = v1_ll - sqrt(equation.gamma * p_ll / rho_ll)
    Ssr = v1_rr + sqrt(equation.gamma * p_rr / rho_rr)
  else # y-direction
    Ssl = v2_ll - sqrt(equation.gamma * p_ll / rho_ll)
    Ssr = v2_rr + sqrt(equation.gamma * p_rr / rho_rr)
  end

  if Ssl >= 0.0 && Ssr > 0.0
    f1 = f_ll[1]
    f2 = f_ll[2]
    f3 = f_ll[3]
    f4 = f_ll[4]
  elseif Ssr <= 0.0 && Ssl < 0.0
    f1 = f_rr[1]
    f2 = f_rr[2]
    f3 = f_rr[3]
    f4 = f_rr[4]
  else
    f1 = (Ssr*f_ll[1] - Ssl*f_rr[1] + Ssl*Ssr*(rho_rr[1]    - rho_ll[1]))/(Ssr - Ssl)
    f2 = (Ssr*f_ll[2] - Ssl*f_rr[2] + Ssl*Ssr*(rho_v1_rr[1] - rho_v1_ll[1]))/(Ssr - Ssl)
    f3 = (Ssr*f_ll[3] - Ssl*f_rr[3] + Ssl*Ssr*(rho_v2_rr[1] - rho_v2_ll[1]))/(Ssr - Ssl)
    f4 = (Ssr*f_ll[4] - Ssl*f_rr[4] + Ssl*Ssr*(rho_e_rr[1]  - rho_e_ll[1]))/(Ssr - Ssl)
  end

  return SVector(f1, f2, f3, f4)
end

# Determine maximum stable time step based on polynomial degree and CFL number
function calc_max_dt(u, element_id, invjacobian, cfl,
                     equation::CompressibleEulerEquations2D, dg)
  λ_max = 0.0
  for j in 1:nnodes(dg), i in 1:nnodes(dg)
    rho, rho_v1, rho_v2, rho_e = get_node_vars(u, dg, i, j, element_id)
    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    v_mag = sqrt(v1^2 + v2^2)
    p = (equation.gamma - 1) * (rho_e - 1/2 * rho * v_mag^2)
    c = sqrt(equation.gamma * p / rho)
    λ_max = max(λ_max, v_mag + c)
  end

  dt = cfl * 2 / (nnodes(dg) * invjacobian * λ_max)

  return dt
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equation::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  p = (equation.gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))

  return SVector(rho, v1, v2, p)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equation::CompressibleEulerEquations2D)
  rho, rho_v1, rho_v2, rho_e = u

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v_square = v1^2 + v2^2
  p = (equation.gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s = log(p) - equation.gamma*log(rho)
  rho_p = rho / p

  w1 = (equation.gamma - s) / (equation.gamma-1) - 0.5 * rho_p * v_square
  w2 = rho_p * v1
  w3 = rho_p * v2
  w4 = -rho_p

  return SVector(w1, w2, w3, w4)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equation::CompressibleEulerEquations2D)
  rho, v1, v2, p = prim
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  rho_e  = p/(equation.gamma-1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2)
  return SVector(rho, rho_v1, rho_v2, rho_e)
end

@inline function density(u, equation::CompressibleEulerEquations2D)
 rho = u[1]
 return rho
end

@inline function pressure(u, equation::CompressibleEulerEquations2D)
 rho, rho_v1, rho_v2, rho_e = u
 p = (equation.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho)
 return p
end

@inline function density_pressure(u, equation::CompressibleEulerEquations2D)
 rho, rho_v1, rho_v2, rho_e = u
 rho_times_p = (equation.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))
 return rho_times_p
end

# Calculates the entropy flux in direction "orientation" and the entropy variables for a state cons
# NOTE: This method seems to work currently (b82534e) but is never used anywhere. Thus it is
# commented here until someone uses it or writes a test for it.
# @inline function cons2entropyvars_and_flux(gamma::Float64, cons, orientation::Int)
#   entropy = MVector{4, Float64}(undef)
#   v = (cons[2] / cons[1] , cons[3] / cons[1])
#   v_square= v[1]*v[1]+v[2]*v[2]
#   p = (gamma - 1) * (cons[4] - 1/2 * (cons[2] * v[1] + cons[3] * v[2]))
#   rho_p = cons[1] / p
#   # thermodynamic entropy
#   s = log(p) - gamma*log(cons[1])
#   # mathematical entropy
#   S = - s*cons[1]/(gamma-1)
#   # entropy variables
#   entropy[1] = (gamma - s)/(gamma-1) - 0.5*rho_p*v_square
#   entropy[2] = rho_p*v[1]
#   entropy[3] = rho_p*v[2]
#   entropy[4] = -rho_p
#   # entropy flux
#   entropy_flux = S*v[orientation]
#   return entropy, entropy_flux
# end

# Calculate thermodynamic entropy for a conservative state `cons`
@inline function entropy_thermodynamic(cons, equation::CompressibleEulerEquations2D)
  # Pressure
  p = (equation.gamma - 1) * (cons[4] - 1/2 * (cons[2]^2 + cons[3]^2) / cons[1])

  # Thermodynamic entropy
  s = log(p) - equation.gamma*log(cons[1])

  return s
end

# Calculate mathematical entropy for a conservative state `cons`
@inline function entropy_math(cons, equation::CompressibleEulerEquations2D)
  # Mathematical entropy
  S = -entropy_thermodynamic(cons, equation) * cons[1] / (equation.gamma - 1)

  return S
end

# Default entropy is the mathematical entropy
@inline entropy(cons, equation::CompressibleEulerEquations2D) = entropy_math(cons, equation)

# Calculate total energy for a conservative state `cons`
@inline energy_total(cons, ::CompressibleEulerEquations2D) = cons[4]

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(cons, equation::CompressibleEulerEquations2D)
  return 0.5 * (cons[2]^2 + cons[3]^2)/cons[1]
end

# Calculate internal energy for a conservative state `cons`
@inline function energy_internal(cons, equation::CompressibleEulerEquations2D)
  return energy_total(cons, equation) - energy_kinetic(cons, equation)
end


text:
@json struct Member
    id::Int
    username::String
    name::String
    state::String
    avatar_url::String
    web_url::String
    expires_at::DateTime
    access_level::Int
end

@json struct Group
    id::Int
    name::String
    path::String
    description::String
end

endpoint(::GitLabAPI, ::typeof(is_member), org::AStr, id::Integer) =
    Endpoint(:GET, "/groups/$(encode(org))/members/$id"; allow_404=true)
postprocessor(::GitLabAPI, ::typeof(is_member)) = DoSomething(ismember)
into(::GitLabAPI, ::typeof(is_member)) = Bool

endpoint(::GitLabAPI, ::typeof(groups)) =
    Endpoint(:GET, "/groups")
into(::GitLabAPI, ::typeof(groups)) = Vector{Group}


text:
"""
Basic model types and methods
"""
module Models

using Slomo.Constants: G
using Slomo.Integrate: solve, integrate

export Model, DensityModel
export update, density, density2d, mass

struct NotImplemented <: Exception end

abstract type Model end

abstract type DensityModel <: Model end

"""
    update(model)

Create a new model from the specified parameters.
"""
function update(model::Model; parameters...)
    args = []
    for kw in propertynames(model)
        if kw in keys(parameters)
            push!(args, parameters[kw])
        else
            push!(args, getproperty(model, kw))
        end
    end
    return typeof(model)(args...)
end

update(model::Array{T} where T<:Model; parameters...) = begin
    [update(m; parameters...) for m in model]
end

"""
    has_analytic_profile(f, model::Model)

Stupid hack to see if there is an analytic profile defined for a subtype
or if it needs to be numerically integrated.
"""
function has_analytic_profile(f, model::Model)
    method = methods(f, (typeof(model), Vararg)).ms[1]
    return occursin(string(typeof(model)), string(method))
end

"""
    density(model::DensityModel, r)

Local volume density at r.
"""
function density(model::DensityModel, r)
    throw(NotImplemented("no defined density profile"))
end

"""
    density2d(model::DensityModel, R)

Local surface density at R.  If not defined for a subtype of `DensityModel`, then calculate
numerically as the Abel transform from the volume density profile.
"""
function density2d(model::DensityModel, R)
    integrand(r) = density(model, r) * r / √(r ^ 2 - R ^ 2)
    return 2 * integrate(integrand, R, rmax)
end

"""
    mass(model::DensityModel, r)

Enclosed mass within r.  If not defined for a subtype of `DensityModel`, then calculate
numerically as the integral of the volume density.
"""
function mass(model::DensityModel, r)
    integrand(x) = 4π * x ^ 2 * density(model, x)
    return integrate(integrand, eps() ^ (1/3), r)
end

mass(model::Array{T} where T<:DensityModel, r) = sum([mass(m, r) for m in model])

"""
    potential(model::DensityModel, r)

Gravitational potential at r.  Equal to the square of the circular velocity.  If not
defined for a subtype of `DensityModel`, then calculate from the enclosed mass.
"""
function potential(model::DensityModel, r)
    return -G * mass(model, r) ./ r
end

end


text:
#################################
# Alerts API                    #
#################################
# Get severe weather alerts from official Government Meteorological Agencies and leading global weather alert providers.

# Contains a centralized function to make calls to the API.
include("request.jl")

function get_alerts_for_location(apikey::String, locationKey::String, language::String, details::Bool)
    # Returns alert data for a specific location. Alert searches require a location key. Please use the Locations API to obtain the location key for your desired location. By default, a truncated version of the alert data is returned. The full object can be obtained by passing "details=true" into the url string.

    url = "http://dataservice.accuweather.com/alerts/v1/"

    # Add required parameters to the URL
    url = url * locationKey
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end
    if details == true url = url * "&details=true" end

    return make_API_call(url)
end


text:
# * Common operators

"""
    OneBodyHamiltonian

The one-body Hamiltonian, may include external fields. It is diagonal
in spin, i.e. it does not couple orbitals of opposite spin.
"""
struct OneBodyHamiltonian <: OneBodyOperator end

Base.show(io::IO, ::OneBodyHamiltonian) = write(io, "ĥ")
Base.show(io::IO, me::OrbitalMatrixElement{1,A,OneBodyHamiltonian,B}) where{A,B} =
    write(io, "(", join(string.(me.a), " "), "|", join(string.(me.b), " "), ")")

"""
    iszero(me::EnergyExpressions.OrbitalMatrixElement{1,<:SpinOrbital{<:Orbital},OneBodyHamiltonian,<:SpinOrbital{<:Orbital}})

The matrix element vanishes if the spin-orbitals do not have the same spin.
"""
Base.iszero(me::OrbitalMatrixElement{1,A,OneBodyHamiltonian,B}) where {A<:SpinOrbital{<:Orbital},B<:SpinOrbital{<:Orbital}} =
    me.a[1].m[2] != me.b[1].m[2]

"""
    FieldFreeOneBodyHamiltonian

The one-body Hamiltonian, with no external fields. It is diagonal in
the orbital symmetry.
"""
struct FieldFreeOneBodyHamiltonian <: OneBodyOperator end

Base.show(io::IO, ::FieldFreeOneBodyHamiltonian) = write(io, "ĥ₀")
Base.show(io::IO, me::OrbitalMatrixElement{1,A,FieldFreeOneBodyHamiltonian,B}) where{A,B} =
    write(io, "(", join(string.(me.a), " "), "|", join(string.(me.b), " "), ")")
Base.iszero(me::OrbitalMatrixElement{1,<:SpinOrbital,FieldFreeOneBodyHamiltonian,<:SpinOrbital}) =
    symmetry(me.a[1]) != symmetry(me.b[1])

"""
    CoulombInteraction

Two-body Hamiltonian, representing the mutual Coulombic repulsion
between two electrons. Is diagonal in spin, i.e. the spin of the
orbitals associated with the same coordinate must be the same.

# Examples

```jldoctest
julia> EnergyExpressions.OrbitalMatrixElement((:a,:b), CoulombInteraction(), (:c,:d))
[a b|c d]

julia> EnergyExpressions.OrbitalMatrixElement((:a,:b), CoulombInteraction(), (:b,:a))
G(a,b)
```
"""
struct CoulombInteraction{O} <: TwoBodyOperator
    o::O # This can be used to indicate e.g. approximate Coulomb interaction
end
CoulombInteraction() = CoulombInteraction(nothing)

Base.hash(g::CoulombInteraction, h::UInt) = hash(g.o, h)

const CoulombPotential{A,B} = ContractedOperator{1,2,1,A,<:CoulombInteraction,B}

Base.show(io::IO, ::CoulombInteraction) = write(io, "ĝ")

function Base.show(io::IO, me::OrbitalMatrixElement{2,A,<:CoulombInteraction,B}) where {A,B}
    if me.a == me.b # Direct interaction
        write(io, "F($(me.a[1]),$(me.a[2]))")
    elseif me.a[1] == me.b[2] && me.a[2] == me.b[1] # Exchange interaction
        write(io, "G($(me.a[1]),$(me.a[2]))")
    else # General case
        write(io, "[", join(string.(me.a), " "), "|", join(string.(me.b), " "), "]")
    end
end

Base.show(io::IO, co::CoulombPotential{A,B}) where {A,B} =
    write(io, "[$(co.a[1])|$(co.b[1])]")

"""
    iszero(me::EnergyExpressions.OrbitalMatrixElement{2,<:SpinOrbital{<:Orbital},CoulombInteraction,<:SpinOrbital{<:Orbital}})

The matrix element vanishes if the (non-relativistic) spin-orbitals
associated with the same coordinate do not have the same spin.
"""
Base.iszero(me::OrbitalMatrixElement{2,A,<:CoulombInteraction,B}) where {A<:SpinOrbital{<:Orbital},B<:SpinOrbital{<:Orbital}} =
    me.a[1].m[2] != me.b[1].m[2] || me.a[2].m[2] != me.b[2].m[2]

export OneBodyHamiltonian, FieldFreeOneBodyHamiltonian, CoulombInteraction, CoulombPotential


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
divide.jl defines methods for the /(x, y) function
"""

# --- Exports

import Base.:(/)

# --- Method definitions

# ------ Specializations involving an AbstractMultivector instance

# M::AbstractMultivector, B::One
# B::One, M::AbstractMultivector
/(M::AbstractMultivector, B::One) = M
/(B::One, M::AbstractMultivector) = nothing  # TODO

# M::AbstractMultivector, B::Zero
# B::Zero, M::AbstractMultivector
/(B::Zero, M::AbstractMultivector) = B
/(M::AbstractMultivector, B::Zero) = Scalar{typeof(norm(M))}(Inf)

# ------ Specializations involving an AbstractBlade instance

/(B::AbstractBlade, C::AbstractBlade) = B * inv(C)

# B::AbstractBlade, C::One
# B::One, C::AbstractBlade
/(B::AbstractBlade, C::One) = B
/(B::One, C::AbstractBlade) = inv(B)

# B::Zero, C::AbstractBlade
/(B::Zero, C::AbstractBlade) = B

# B::AbstractBlade, x::Real
# x::Real, B::AbstractBlade
/(B::AbstractBlade, x::Real) = Blade(B, volume=volume(B) / x)
/(x::Real, B::AbstractBlade) = x * inv(B)

# ------ Specializations involving a Blade instance

# B::Blade, C::AbstractScalar
# B::AbstractScalar, C::Blade
/(B::Blade, C::AbstractScalar) = Blade(B, volume=volume(B) / value(C))
/(B::AbstractScalar, C::Blade) =
    mod(grade(C), 4) < 2 ?
        Blade(C, volume=value(B) / volume(C)) :
        Blade(C, volume=-value(B) / volume(C))

# B::Blade, C::One
# B::One, C::Blade
/(B::Blade, C::One) = B
/(B::One, C::Blade) = inv(C::Blade)

# B::Blade, C::Zero
# B::Zero, C::Blade
/(B::Blade, C::Zero) = Blade(B, volume=sign(B) * Inf)
/(B::Zero, C::Blade) = B

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
function /(B::Pseudoscalar, C::Pseudoscalar)
    assert_dim_equal(B, C)
    Scalar{typeof(value(B))}(value(B) / value(C))
end

# B::Pseudoscalar, C::AbstractScalar
# B::AbstractScalar, C::Pseudoscalar
/(B::Pseudoscalar, C::AbstractScalar) = B / value(C)
/(B::AbstractScalar, C::Pseudoscalar) = value(B) / C

# B::Pseudoscalar, C::One
# B::One, C::Pseudoscalar
/(B::Pseudoscalar, C::One) = B
/(B::One, C::Pseudoscalar) = inv(C)

# B::Pseudoscalar, C::Zero
# B::Zero, C::Pseudoscalar
/(B::Pseudoscalar, C::Zero) = Pseudoscalar(B, value=sign(B) * Inf)
/(B::Zero, C::Pseudoscalar) = B

# B::Pseudoscalar, x::Real
# x::Real, B::Pseudoscalar
/(B::Pseudoscalar, x::Real) = Pseudoscalar(B, value=value(B) / x)
/(x::Real, B::Pseudoscalar) =
    mod(dim(B), 4) < 2 ?
        Pseudoscalar(B, value=x / value(B)) :
        Pseudoscalar(B, value=-x / value(B))

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::One
# B::One, C::AbstractScalar
/(B::AbstractScalar, C::One) = B
/(B::One, C::AbstractScalar) = inv(C)

# B::AbstractScalar, C::Zero
# B::Zero, C::AbstractScalar
/(B::AbstractScalar, C::Zero) =
    sign(B) > 0 ?
        Scalar{typeof(norm(B))}(Inf) :
        Scalar{typeof(norm(B))}(-Inf)

/(B::Zero, C::AbstractScalar) = B

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
/(B::AbstractScalar, x::Real) = Scalar{typeof(value(B))}(value(B) / x)
/(x::Real, B::AbstractScalar) = Scalar{typeof(value(B))}(x / value(B))

# B::AbstractScalar, v::Vector
# v::Vector, B::AbstractScalar
/(v::Vector{<:Real}, B::AbstractScalar) = Blade(v) / value(B)

# ------ Specializations involving a Scalar instance

# B::Scalar, C::Scalar
/(B::Scalar, C::Scalar) = Scalar{typeof(value(B))}(value(B) / value(C))

# ------ Specializations involving a One instance

# B::One, C::One
/(B::One, C::One) = B

# B::One, C::Zero
# B::Zero, C::One
/(B::One, C::Zero) = inv(C)
/(B::Zero, C::One) = B

# ------ Specializations involving a Zero instance

# B::Zero, C::Zero
/(B::Zero, C::Zero) = Scalar{typeof(value(B))}(NaN)


text:
using LaTeXStrings
import MathTeXEngine:
    generate_tex_elements, inkwidth, inkheight, bottominkbound, TeXChar, HLine

"""
    texalign(halign, valign, bottom_pt, top_pt, font_size)
Helper function to align LaTeX text properly. Returns
`translate_x` and `translate_y` which consists of the amount
to be shifted depending on the type of alignment chosen
and the bounding box of the text.
"""
function texalign(halign, valign, bottom_pt::Point, top_pt::Point, font_size)
    textw = top_pt[1] - bottom_pt[1]
    translate_x, translate_y = 0, 0

    if halign === :left
        translate_x = 0
    elseif halign === :right
        translate_x = -textw * font_size
    elseif halign === :center
        translate_x = -textw / 2 * font_size
    end

    if valign === :baseline
        translate_y = 0
    elseif valign === :bottom
        translate_y = 0
    elseif valign === :top
        translate_y = -top_pt[2] * font_size
    elseif valign === :middle
        translate_y = -(bottom_pt[2] + top_pt[2]) * font_size / 2
    end
    return translate_x, translate_y
end

"""
    rawlatexboundingbox(lstr::LaTeXString, font_size=1)
Helper function that returns the coordinate points
of the bounding box containing the specific LaTeX text.
"""
function rawlatexboundingbox(lstr::LaTeXString)
    sentence = generate_tex_elements(lstr)
    els = filter(x -> x[1] isa TeXChar, sentence)
    texchars = [x[1] for x in els]
    pos_x = [x[2][1] for x in els]
    pos_y = [x[2][2] for x in els]
    scales = [x[3] for x in els]
    textw = maximum(inkwidth.(texchars) .* scales .+ pos_x)
    bottom = bottominkbound.(texchars)
    bottom_pt = Point(0, maximum(-bottom .* scales .- pos_y))
    top_pt = Point(textw, -maximum(inkheight.(texchars) .* scales .+ pos_y .+ bottom))

    return (bottom_pt, top_pt)
end

"""
    latexboundingbox(lstr::LaTeXString, font_size=get_fontsize(); halign=:left, valign=:right)
Returns the bounding box containing the latex text with
`(Lower Left Point, Upper Right Point)`.
Use `box(latex_bb(testext)...,:stroke)` to draw the bounding box.
"""
function latexboundingbox(
    lstr::LaTeXString, font_size=get_fontsize(); halign=:left, valign=:right
)
    bottom_pt, top_pt = rawlatexboundingbox(lstr)

    translate_x, translate_y = texalign(halign, valign, bottom_pt, top_pt, font_size)
    translate_pt = Point(translate_x, translate_y)

    return (bottom_pt * font_size + translate_pt, top_pt * font_size + translate_pt)
end

"""
    latextextsize(lstr::LaTeXString)
Returns the width and height of a latex string.
"""
function latextextsize(lstr::LaTeXString)
    bottom_pt, top_pt = latexboundingbox(lstr)
    textw = top_pt[1] - bottom_pt[1]
    texth = -(top_pt[2] - bottom_pt[2])
    return textw, texth
end

"""
    text(lstr::LaTeXString, pt::Point; valign=:baseline,
    halign=:left, rotationfixed = false, angle=0, kwargs...)

Draws LaTeX string using `MathTexEngine.jl`. Hence, uses
ModernCMU as font family. When `rotationfixed = true`,
the text will rotate around its own axis, instead of
rotating around `pt`.
"""
function text(
    lstr::LaTeXString,
    pt::Point;
    valign=:baseline,
    halign=:left,
    angle=0::Real,
    rotationfixed=false,
    kwargs...,
)

    # Function from MathTexEngine
    sentence = generate_tex_elements(lstr)

    # Get current font size.
    font_size = get_fontsize()

    textw, texth = latextextsize(lstr)
    bottom_pt, top_pt = rawlatexboundingbox(lstr)

    translate_x, translate_y = texalign(halign, valign, bottom_pt, top_pt, font_size)
    # translate_x, translate_y = texalign(halign, valign, textw, texth, font_size)

    # Writes text using ModernCMU font.
    for text in sentence
        @layer begin
            translate(pt)
            if !rotationfixed
                rotate(angle)
                translate(translate_x, translate_y)
            else
                l_pt, r_pt = latexboundingbox(lstr, halign = halign, valign = valign)
                translate((l_pt + r_pt)/2)
                rotate(angle)
                translate(Point(translate_x,translate_y)-(l_pt + r_pt)/2)
            end
            if text[1] isa TeXChar
                fontface(text[1].font.family_name)
                fontsize(font_size * text[3])
                Luxor.text(string(text[1].char), Point(text[2]...) * font_size * (1, -1))

            elseif text[1] isa HLine
                pointstart = Point(text[2]...) * font_size * (1, -1)
                pointend = pointstart + Point(text[1].width, 0) * font_size
                setline(0.5)
                line(pointstart, pointend, :stroke)
            end
        end
    end
end


text:
"""
   plot_decision_boundary_on_scatter(X, Y, Ŷ)

Plot a 2-dimensions scatter using X and Y and add the decision boundary computed in Ŷ

"""
function plot_decision_boundary_on_scatter(X, Y, Ŷ, title)

    # Create plot for training data (X matrix)
    training_plot = scatter(
        X[:,2],
        X[:,3],
        group = Y,
        palette = training_palette,
        ticks = false,
        legend = false,
        framestyle = :box,
        title = title,
    )

    # Prepare data to show the decision boundary
    # @TODO is there any more convenient way to do it?
    step = 0.01
    x1 = range(minimum(X[:,2]), maximum(X[:,2]), step=step)
    x2 = range(minimum(X[:,3]), maximum(X[:,3]), step=step)

    # Add the decision boundary at Ŷ = 0.5
    contour!(
        training_plot,
        x1,
        x2,
        Ŷ,
        levels = [0.5],
        color = :black,
    )

end


text:
module Random

using Base.LibRandom

export srand,
       rand, rand!,
       randn, randn!,
       randbool, randbool!,
       AbstractRNG, RNG, MersenneTwister

abstract AbstractRNG

type MersenneTwister <: AbstractRNG
    state::DSFMT_state
    seed::Union(Uint32,Vector{Uint32})

    function MersenneTwister(seed::Vector{Uint32})
        state = DSFMT_state()
        dsfmt_init_by_array(state, seed)
        return new(state, seed)
    end

    MersenneTwister(seed=0) = MersenneTwister(make_seed(seed))
end

function srand(r::MersenneTwister, seed)
    r.seed = seed
    dsfmt_init_gen_rand(r.state, seed)
    return r
end

## initialization

function __init__()

@unix_only begin
    try
        srand("/dev/urandom")
    catch
        println(STDERR, "Entropy pool not available to seed RNG; using ad-hoc entropy sources.")
        seed = reinterpret(Uint64, time())
        seed = bitmix(seed, uint64(getpid()))
        try
            seed = bitmix(seed, parseint(Uint64, readall(`ifconfig` |> `sha1sum`)[1:40], 16))
        catch
            # ignore
        end
        srand(seed)
    end
end

@windows_only begin
    a = zeros(Uint32, 2)
    win32_SystemFunction036!(a)
    srand(a)
end
end

## srand()

function srand(seed::Vector{Uint32})
    global RANDOM_SEED = seed
    dsfmt_gv_init_by_array(seed)
end
srand(n::Integer) = srand(make_seed(n))

function make_seed(n::Integer)
    n < 0 && throw(DomainError())
    seed = Uint32[]
    while true
        push!(seed, n & 0xffffffff)
        n >>= 32
        if n == 0
            return seed
        end
    end
end

function srand(filename::String, n::Integer)
    open(filename) do io
        a = Array(Uint32, int(n))
        read!(io, a)
        srand(a)
    end
end
srand(filename::String) = srand(filename, 4)

## random floating point values

rand(::Type{Float64}) = dsfmt_gv_genrand_close_open()
rand() = dsfmt_gv_genrand_close_open()

rand(::Type{Float32}) = float32(rand())
rand(::Type{Float16}) = float16(rand())

rand{T<:Real}(::Type{Complex{T}}) = complex(rand(T),rand(T))

rand(r::MersenneTwister) = dsfmt_genrand_close_open(r.state)

## random integers

dsfmt_randui32() = dsfmt_gv_genrand_uint32()
dsfmt_randui64() = uint64(dsfmt_randui32()) | (uint64(dsfmt_randui32())<<32)

rand(::Type{Uint8})   = uint8(rand(Uint32))
rand(::Type{Uint16})  = uint16(rand(Uint32))
rand(::Type{Uint32})  = dsfmt_randui32()
rand(::Type{Uint64})  = dsfmt_randui64()
rand(::Type{Uint128}) = uint128(rand(Uint64))<<64 | rand(Uint64)

rand(::Type{Int8})    = int8(rand(Uint8))
rand(::Type{Int16})   = int16(rand(Uint16))
rand(::Type{Int32})   = int32(rand(Uint32))
rand(::Type{Int64})   = int64(rand(Uint64))
rand(::Type{Int128})  = int128(rand(Uint128))

# Arrays of random numbers

rand!(A::Array{Float64}) = dsfmt_gv_fill_array_close_open!(A)
rand(::Type{Float64}, dims::Dims) = rand!(Array(Float64, dims))
rand(::Type{Float64}, dims::Int...) = rand(Float64, dims)

rand(dims::Dims) = rand(Float64, dims)
rand(dims::Int...) = rand(Float64, dims)

rand!(r::MersenneTwister, A::Array{Float64}) = dsfmt_fill_array_close_open!(r.state, A)
rand(r::AbstractRNG, dims::Dims) = rand!(r, Array(Float64, dims))
rand(r::AbstractRNG, dims::Int...) = rand(r, dims)

function rand!{T}(A::Array{T})
    for i=1:length(A)
        A[i] = rand(T)
    end
    A
end
rand(T::Type, dims::Dims) = rand!(Array(T, dims))
rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try a more specific type")
rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims)

# Generate random integer within a range

immutable RandIntGen{T<:Integer, U<:Unsigned}
    a::T   # first element of the range
    k::U   # range length
    u::U   # maximum multiple of k within the domain of U

    RandIntGen(a::T, k::U) = new(a, k, div(typemax(U),k)*k)
end

RandIntGen{T<:Unsigned}(r::UnitRange{T}) = RandIntGen{T,T}(first(r), convert(T, length(r)))

# specialized versions
for (T, U) in [(Uint8, Uint32), (Uint16, Uint32), (Int8, Uint32), (Int16, Uint32),
               (Int32, Uint32), (Int64, Uint64), (Int128, Uint128),
               (Bool, Uint32), (Char, Uint32)]

    @eval RandIntGen(r::UnitRange{$T}) = RandIntGen{$T, $U}(first(r), convert($U, length(r)))
end

function rand{T<:Integer,U<:Unsigned}(g::RandIntGen{T,U})
    x = rand(U)
    while x >= g.u
        x = rand(U)
    end
    convert(T, g.a + rem(x, g.k))
end

rand{T<:Union(Signed,Unsigned,Bool,Char)}(r::UnitRange{T}) = rand(RandIntGen(r))
rand{T<:Real}(r::Range{T}) = convert(T, first(r) + rand(0:(length(r)-1)) * step(r))

function rand!(g::RandIntGen, A::AbstractArray)
    for i = 1 : length(A)
        @inbounds A[i] = rand(g)
    end
    return A
end

rand!{T<:Union(Signed,Unsigned,Bool,Char)}(r::UnitRange{T}, A::AbstractArray) = rand!(RandIntGen(r), A)

function rand!{T<:Real}(r::Range{T}, A::AbstractArray)
    g = RandIntGen(0:(length(r)-1))
    f = first(r)
    s = step(r)
    if s == 1
        for i = 1 : length(A)
            @inbounds A[i] = convert(T, f + rand(g))
        end
    else
        for i = 1 : length(A)
            @inbounds A[i] = convert(T, f + rand(g) * s)
        end
    end
    return A
end

rand{T<:Real}(r::Range{T}, dims::Dims) = rand!(r, Array(T, dims))
rand(r::Range, dims::Int...) = rand(r, dims)

## random Bools

rand!(B::BitArray) = Base.bitarray_rand_fill!(B)

randbool(dims::Dims) = rand!(BitArray(dims))
randbool(dims::Int...) = rand!(BitArray(dims))

randbool() = ((dsfmt_randui32() & 1) == 1)
randbool!(B::BitArray) = rand!(B)

## randn() - Normally distributed random numbers using Ziggurat algorithm

# The Ziggurat Method for generating random variables - Marsaglia and Tsang
# Paper and reference code: http://www.jstatsoft.org/v05/i08/

randn() = randmtzig_gv_randn()
randn(rng::MersenneTwister) = randmtzig_randn(rng.state)
randn!(A::Array{Float64}) = randmtzig_gv_fill_randn!(A)
randn!(rng::MersenneTwister, A::Array{Float64}) = randmtzig_fill_randn!(rng.state, A)
randn(dims::Dims) = randn!(Array(Float64, dims))
randn(dims::Int...) = randn!(Array(Float64, dims...))

## random UUID generation

immutable UUID
    value::Uint128
end

function uuid4()
    u = rand(Uint128)
    u &= 0xffffffffffff0fff3fffffffffffffff
    u |= 0x00000000000040008000000000000000
    UUID(u)
end

function Base.convert(::Type{Vector{Uint8}}, u::UUID)
    u = u.value
    a = Array(Uint8,36)
    for i = [36:-1:25; 23:-1:20; 18:-1:15; 13:-1:10; 8:-1:1]
        d = u & 0xf
        a[i] = '0'+d+39*(d>9)
        u >>= 4
    end
    a[[24,19,14,9]] = '-'
    return a
end

Base.show(io::IO, u::UUID) = write(io,convert(Vector{Uint8},u))
Base.repr(u::UUID) = ASCIIString(convert(Vector{Uint8},u))

end # module

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using RLESUtils,RunUtils,IFTTTUtils

const COMP = gethostname()
const DATA = "098filt_10K"

template_gp(i) =
"""
using GrammarExpts,ACASX_GP_Tree,RLESUtils,IFTTTUtils
config=configure(ACASX_GP_Tree,"normal","nvn_libcas098smallfilt_10K")
acasx_gp_tree(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare_Tree/ACASX_GP_Tree/seed$i"), config...)
sendifttt(;value1="gp_tree,$i", value2=$COMP, value3=$DATA)
"""

template_ge(i) =
"""
using GrammarExpts,ACASX_GE_Tree,RLESUtils,IFTTTUtils
config=configure(ACASX_GE_Tree,"normal","nvn_libcas098smallfilt_10K")
acasx_ge_tree(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare_Tree/ACASX_GE_Tree/seed$i"), config...)
sendifttt(;value1="ge_tree,$i", value2=$COMP, value3=$DATA)
"""

template_mc(i) =
"""
using GrammarExpts,ACASX_MC_Tree,RLESUtils,IFTTTUtils
config=configure(ACASX_MC_Tree,"normal","nvn_libcas098smallfilt_10K")
acasx_mc_tree(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare_Tree/ACASX_MC_Tree/seed$i"), config...)
sendifttt(;value1="mc_tree,$i", value2=$COMP, value3=$DATA)
"""

template_mcts(i) =
"""
using GrammarExpts,ACASX_MCTS_Tree,RLESUtils,IFTTTUtils
config=configure(ACASX_MCTS_Tree,"normal","nvn_libcas098smallfilt_10K")
acasx_mcts(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare_Tree/ACASX_MCTS_Tree/seed$i"), config...)
sendifttt(;value1="mcts_tree,$i", value2=$COMP, value3=$DATA)
"""

notifydone() = sendifttt(;value1="$COMP done")

A = JuliaSource[]
append!(A, [JuliaSource(template_ge(i)) for i=11:13])
append!(A, [JuliaSource(template_gp(i)) for i=11:13])
append!(A, [JuliaSource(template_mc(i)) for i=11:13])
append!(A, [JuliaSource(template_mcts(i)) for i=1:3])
append!(A, [JuliaSource(template_mcts(i)) for i=11:13])


text:
"""
    check_vcf_sites("infile.vcf.gz", "datafile.txt"; c=1)

Check that variant sites (e.g. 9:43985033) listed in a tab-delimited file
are found in a vcf file. List the sites not found in the
vcf file.  Write a new file that contains only the sites found
in the vcf file.

Options:

c=n    Specify the column number n that contains the sites [1].

"""
function check_vcf_sites(infile::String, datafile::String; c=1)

    print("Reading all sites from $datafile ... ")
    filein = readdlm(datafile, '\t')
    filesites = Set( readdlm(datafile, '\t')[:,c] )
    println("(e.g. ", rand(filesites), ")" , " ... done.")

    print("Reading all sites from $infile ... ")
    vcfsites = Set( split( read(`bcftools query -f '%CHROM:%POS,' $infile`, String), ',' ) )
    println("done.")

    fileonly = setdiff(filesites, vcfsites)
    isect = intersect(filesites, vcfsites)

    ni = length(isect)
    n = length(fileonly)

    println("Found $ni sites in vcf file.")
    println("Found $n sites not in vcf file:")
    println(fileonly)

    gsites = filein[ findall(in(isect), filein[:,c]), : ]

    writedlm("$datafile.filtered", gsites, '\t')
    println("Wrote filtered file to $datafile.filtered.")

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
# To determine time of signal transmission, this function must be solved

# Returns the difference between the distance light travels in signal reception time (t) - signal transmission time (tt)
# and the geometric distance between the two satelites
function transmissionTimeFunction(tt, t, r, X, Y, Z)
    # Speed of light
    c = 299792.458 # [km / s]

    rsat = @SVector [X(tt), Y(tt), Z(tt)]
    func = c*(t - tt) - norm(r .- rsat)

    return func
end

# Determin transmission time through itteration
#  https://gssc.esa.int/navipedia/index.php/Emission_Time_Computation
function transmissionTimeFunction(t,r,X,Y,Z)
    # Speed of light
    c = 299792.458 # [km / s]

    # Convergence criteria
    iMax = 100
    tol  = 1e-12

    rsat    = @SVector [X(t), Y(t), Z(t)]
    stop    = false
    iter    = 0
    tt      = 0.0
    while stop == false && iter < iMax
        iter += 1

        Δt = sqrt((r[1] - rsat[1])^2 + (r[2] - rsat[2])^2 + (r[3] - rsat[3])^2) / c
        tt = t - Δt

        rdiff = sqrt((rsat[1] - X(tt))^2 + (rsat[2] - Y(tt))^2 + (rsat[3] - Z(tt))^2)
        if rdiff < tol
            stop = true
        else
            rsat = @SVector [X(tt), Y(tt), Z(tt)]
        end
    end

    if iter == iMax
        throw(ErrorException("Transmission time function did not converge."))
    end

    return tt
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
## shows some data analysis using DataFrames
using Econometrics, CSV, DataFrames, GLM
cd(@__DIR__)
nerlove = DataFrame(CSV.File("../Data/nerlove.csv"))
f = @formula(log(cost) ~ 1 + log(output) + log(labor) + log(fuel) + log(capital))

##
println("using the GLM package")
@show lm(f, nerlove)

##
println("using ols from this package")
println("note: this uses White's st. errors, GLM uses plain vanilla")
ols(f, nerlove)
nothing

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics