text:
#simulation of subsidizing
using Optim, LinearAlgebra, Random, Distributions, Suppressor, Statistics

struct UnobservedType
    v::Array{Float64}
    c::Array{Float64, 2}
end

function UnobservedType(n)
    v = rand(Uniform(7, 15), n)
    cW = rand(Uniform(5, 10), n)
    cM = rand(Uniform(2, 5), n)
    return UnobservedType(v, [cW cM])
end

function calculate_ground_truth()
    n = 1000000
    types = UnobservedType(n)
    p0 = eqm_price(zeros(n, 2), types)
    function direct_effect(η)
        y = get_supplier_welfare(transpose(p0) .+ zeros(n, 2), types; η = η)
        return mean(y)
    end

    function total_effect(η)
        p1 = eqm_price(zeros(n, 2), types; η = η)
        y1 = get_supplier_welfare(transpose(p1) .+ zeros(n, 2), types; η = η)
        return mean(y1)
    end
    #println(eqm_price(zeros(n, 2), types; η = 0.1) - eqm_price(zeros(n, 2), types; η = 0.0))
    println((total_effect(0.2) - total_effect(-0.2))/0.4)
    println((direct_effect(0.2) - direct_effect(-0.2))/0.4)
    #τ = central_fdm(5, 1; factor=1e6)(total_effect, 0.0)
    #τD = central_fdm(5, 1; factor=1e6)(direct_effect, 0.0)

    #return (τ = τ, τD = τD)
end

function mc_simulation(S, n)
    results = zeros((S, 7))
    for i in 1:S
        estimates = run_simulation(n)
        results[i, :] .= collect(estimates)
    end
    return results
end

function run_simulation(n)
    types = UnobservedType(n)
    return calculate_effects(types)
end

function get_supplier_welfare(P, types; η= 0.0)
    qs = aggregate_supply(P, types; η=η)
    return qs[:, 1].*(P[:, 1] .+ η .- P[:, 2] .- types.c[:, 1])
end

# calculate general eqm effect using local experimemnt
function calculate_effects(types)
    n = length(types.v)
    U = 0.6 .*rand([-1,1], (n, 2))
    pstar = eqm_price(U, types)
    Z = aggregate_excess_demand(transpose(pstar) .+ U, types)
    qs = aggregate_supply(transpose(pstar) .+ U, types)
    qd = aggregate_demand(transpose(pstar) .+ U, types)
    ys = get_supplier_welfare(transpose(pstar) .+ U, types)
    ∂Zπ = zeros(2, 1)
    ∂Zp = zeros(2, 2)
    ∂ys = inv(transpose(U)*U)*(transpose(U)*ys)

    for i in 1:2
        ∂Zp[:, i] .= inv(transpose(U)*U)*(transpose(U)*Z[:, i])
    end
    ∂Zπ[1] = -1.0 .* inv(transpose(U[:, 1])*U[:, 1])*(transpose(U[:, 1])*qs[:, 1])
    ∂Zπ[2] = inv(transpose(U[:, 1])*U[:, 1])*(transpose(U[:, 1])*Z[:, 2])
    ∂pstar = - inv(∂Zp)*(∂Zπ)
    τI = transpose(∂ys)*∂pstar
    dy = ∂ys[1] + τI[1]
    return (τ = dy, τD  = ∂ys[1], τI = τI[1], ∂pstarG = ∂pstar[1], ∂pstarH = ∂pstar[2],
            ∂ysG = ∂ys[1], ∂ysW = ∂ys[2])
end

function aggregate_demand(P, types)
    n = size(types.c)[1]
    qd = zeros((n, 2))
    for i in 1:n
        qd[i, :] .=  individual_excess_demand(P[i, :], types.c[i, :], types.v[i]).qd
    end
    return qd
end

function aggregate_supply(P, types; η=0.0)
    n = size(types.c)[1]
    qs = zeros((n, 2))
    for i in 1:n
        qs[i, :] .=  individual_excess_demand(P[i, :], types.c[i, :], types.v[i]; η= η).qs
    end
    return qs
end

function aggregate_excess_demand(P, types; η = 0.0)
    n = size(types.c)[1]
    Z = zeros((n, 2))
    for i in 1:n
        Z[i, :] .=  individual_excess_demand(P[i, :], types.c[i, :], types.v[i]; η = η).z
    end
    return Z
end

function individual_excess_demand(p, c, v; η = 0.0)
     #println(p)
     profits = [p[1] + η - p[2] - c[1], p[2] - c[2]]
     k = argmax(profits)
     l = argmin(profits)
     z = zeros(2)
     z[k] = -1.0 .* (profits[k] .>0)
     z[l] = 0.0
     supply = -1.0 .* z
     demand = [(p[1] < v).*1.0, supply[1]]
     z[2] = z[2] + demand[2]
     z[1] = z[1] + demand[1]
     return (z= z, qs = supply, qd = demand)
 end

function eqm_price(U, types::UnobservedType; η = 0.0)
    function eqm_condition(p)
        zp = aggregate_excess_demand(transpose(p) .+ U, types; η= η)
        return norm(mean(zp, dims=1))
    end
    r = optimize(eqm_condition, [8.0, 3.0]).minimizer
    #rintln(r)
    return r
end

max_stars_repo_name: evanmunro/market-equilibrium

max_issues_repo_name: evanmunro/market-equilibrium

max_forks_repo_name: evanmunro/market-equilibrium


text:
struct RiskFunctions{T <: DiseaseStateSequence}
  sparks::Union{Nothing, Function}
  susceptibility::Union{Nothing, Function}
  infectivity::Union{Nothing, Function}
  transmissibility::Union{Nothing, Function}
  latency::Union{Nothing, Function}
  removal::Union{Nothing, Function}

  RiskFunctions{SEIR}(ϵ, Ωs, κ, Ωt, Ωl, Ωr) = new{SEIR}(ϵ, Ωs, κ, Ωt, Ωl, Ωr)
  RiskFunctions{SEI}(ϵ, Ωs, κ, Ωt, Ωl) = new{SEI}(ϵ, Ωs, κ, Ωt, Ωl, nothing)
  RiskFunctions{SIR}(ϵ, Ωs, κ, Ωt, Ωr) = new{SIR}(ϵ, Ωs, κ, Ωt, nothing, Ωr)
  RiskFunctions{SI}(ϵ, Ωs, κ, Ωt) = new{SI}(ϵ, Ωs, κ, Ωt, nothing, nothing)
end

function Base.show(io::IO, x::RiskFunctions{T}) where T <: DiseaseStateSequence
  return print(io, "$T model risk functions")
end

ext: jl


text:
# Ray - 3D Plane Intersection
# Returns: hit, distance
# hit: a boolean indicating if an intersection occurred
# distance: the distance to intersection
# Edge case policies:
#    1. Ray is inside the plane: Returns (true, NaN)
#    2. Ray is parallel to the plan, but not inside: Returns (false, NaN)
#    3. Ray never hits plane: Returns (false, NaN)
"""
    function raytrace(ray::Ray, surface::Surface)

Determines if a `Ray` and a `Surface` intersect, and the distance to that intersection.

# Returns
* `Bool`: Indicates if the ray intersects the surface or not
* `Float64`: The distance between the ray's origin and the point of intersection
"""
function raytrace(ray::Ray, plane::Plane)
    dist::Float64 = dot( plane.point - ray.origin, plane.normal) / dot( ray.direction, plane.normal)
    # Check if parallel
    if dist < 0 || dist == Inf
        return false, NaN
    end
    return true, dist
end

# Ray - 3D Cone Intersection
# Returns hit, distance
# hit: a boolean indicating if an intersection occurred (false if parallel or negative)
# dist: distance to closest intersection point
function raytrace(ray::Ray, cone::Cone)

	cos_theta_squared::Float64 = (cos(cone.theta))^2
	CO::Coord = ray.origin - cone.tip

	a::Float64 = dot(ray.direction, cone.axis)^2 - cos_theta_squared
	b::Float64 = 2.0 * (dot(ray.direction, cone.axis) * dot(CO, cone.axis) - dot(ray.direction, CO) * cos_theta_squared)
	c::Float64 = dot(CO, cone.axis)^2 - dot(CO, CO) * cos_theta_squared

	determinant::Float64 = b^2 - 4.0*a*c

	if determinant < 0
		return false, NaN
	elseif determinant == 0
		return true, -b/(2*a)
	end

	# Now we need to verify we are not intersecting with the shadow cone
	one_over_two_a::Float64 = 1.0 / (2.0 * a)
	t1::Float64 = (-b - sqrt(determinant)) * one_over_two_a
	t2::Float64 = (-b + sqrt(determinant)) * one_over_two_a

	p1::Coord = ray.origin + ray.direction * t1
	p2::Coord = ray.origin + ray.direction * t2

	valid1::Bool = false
	valid2::Bool = false

	# Check against shadow cone & ensure intersection is in front of ray
	if t1 >= 0 && dot((p1-cone.tip),cone.axis) > 0
		valid1 = true
	end
	if t2 >= 0 && dot((p2-cone.tip),cone.axis) > 0
		valid2 = true
	end

	# If both points hit real cone, return closer one
	if valid1 && valid2
		if t1 < t2
			return true, t1
		else
			return true, t2
		end
	end

	if valid1 && !valid2
		return true, t1
	end

	if !valid1 && valid2
		return true, t2
	end

	return false, NaN
end

# Ray - 3D Sphere Intersection
# Returns hit, distance
# hit: a boolean indicating if an intersection occurred (false if parallel or negative)
# dist: distance to closest intersection point
function raytrace(ray::Ray, sphere::Sphere)
    d::Coord = ray.origin - sphere.center
    t::Float64 = -dot(ray.direction, d)
    discriminant::Float64 = t^2
    discriminant -= magnitude(d)^2
    discriminant += sphere.radius^2

    # If the discriminant is less than zero, they don't hit
    if discriminant < 0
        return false, Inf
    end
    sqrt_val::Float64 = sqrt(discriminant)
    pos::Float64 = t - sqrt_val
    neg::Float64 = t + sqrt_val

    if pos < 0  && neg < 0
        return false, NaN
    end
    if pos < 0 && neg > 0
        return true, neg
    end
    if pos < neg && pos > 0
        return true, pos
    end

    return true, neg
end

# Ray - 3D Infinite Cylinder Intersection (works for cylinder direction)
# Returns hit, distance
# hit: a boolean indicating if an intersection occurred (false if parallel or negative)
# dist: distance to closest intersection point
function raytrace(ray::Ray, infcylinder::InfCylinder)
    A = infcylinder.center
    # Generate point new point in cylinder for math
    B = infcylinder.center + infcylinder.normal
    O = ray.origin
    r = infcylinder.radius
    AB = B - A
    AO = O - A
    AOxAB = cross(AO, AB)
    VxAB = cross(ray.direction, AB)
    ab2::Float64 = dot(AB, AB)
    a::Float64 = dot(VxAB, VxAB)
    b::Float64 = 2.0 * dot(VxAB, AOxAB)
    c::Float64 = dot(AOxAB, AOxAB) - (r*r * ab2)

    # Check Determininant
    det::Float64 = b^2 - 4.0 * a * c
    if det < 0
        return false, Inf
    end

    pos::Float64 = (-b + sqrt(det)) / (2.0 * a)
    neg::Float64 = (-b - sqrt(det)) / (2.0 * a)

    if pos < 0
        if neg < 0
            return false, NaN
        end
        return true, neg
    end
    if neg < 0
        return true, pos
    end
    if pos < neg
        return true, pos
    else
        return true, neg
    end
end

"""
    reflect(ray::Ray, plane::Plane)

Reflects a ray off a plane.

# Return
* `Ray`: A new ray with the same origin as input, but with the new reflected direction
"""
function reflect(ray::Ray, plane::Plane)
    a = dot(ray.direction, plane.normal)
    b = plane.normal * (2.0 * a)
    c = ray.direction - b
    reflected_ray::Ray = Ray(ray.origin, c)
    return reflected_ray
end

"""
    generate_random_ray(box::Box)

Returns a randomly sampled ray from within an axis aligned bounding box.
"""
function generate_random_ray(box::Box)
    ray = Ray(Coord(0.0, 0.0, 0.0), Coord(0.0, 0.0, 0.0))

    # Sample Origin
    width = Coord(box.upper_right.x - box.lower_left.x, box.upper_right.y - box.lower_left.y, box.upper_right.z - box.lower_left.z)
    ray.origin.x = box.lower_left.x + rand(Float64)*width.x
    ray.origin.y = box.lower_left.y + rand(Float64)*width.y
    ray.origin.z = box.lower_left.z + rand(Float64)*width.z

    # Sample Direction From Sphere
    theta::Float64 = rand(Float64) * 2.0 * pi
    z::Float64 = -1.0 + 2.0 * rand(Float64)
    zo::Float64 = sqrt(1.0 - z*z)
    ray.direction.x = zo * cos(theta);
    ray.direction.y = zo * sin(theta);
    ray.direction.z = z;

    # Normalize Direction
    ray.direction = unitize(ray.direction)

    return ray
end

"""
    find_intersection(ray::Ray, regions::Array{Region})

Performs ray tracing on an array of regions.

# Return
* `Ray`: A new Ray that has been moved just accross the point of intersection.
* `Int64`: The surface id that was hit.
* `String`: The boundary condition of the surface that was hit.
"""
function find_intersection(ray::Ray, regions::Array{Region})
    BUMP::Float64 = 1.0e-9
    min::Float64 = 1e30
    id::Int64 = -1
    for i = 1:length(regions)
        hit, dist = raytrace(ray, regions[i].surface)
        if hit == true
            if dist < min
                min = dist
                id = i
            end
        end
    end

    new_ray::Ray = Ray(ray.origin + ray.direction * (min + BUMP), ray.direction)

    if regions[id].surface.reflective == true
        new_ray = reflect(new_ray, regions[id].surface)
        new_ray.origin = new_ray.origin + new_ray.direction * (2.0 * BUMP)
		return new_ray, id, "reflective"
    end

	if regions[id].surface.vacuum == true
		return new_ray, id, "vacuum"
    end

    return new_ray, id, "transmission"

end

"""
    find_intersection(ray::Ray, geometry::Geometry)

Performs ray tracing on a Geometry

# Return
* `Ray`: A new Ray that has been moved just accross the point of intersection.
* `Int64`: The surface id that was hit.
* `String`: The boundary condition of the surface that was hit.
"""
function find_intersection(ray::Ray, geometry::Geometry)
	cell_id = find_cell_id(ray.origin, geometry)
	regions::Array{Region} = geometry.cells[cell_id].regions
	return find_intersection(ray, regions)
end


text:
const chainrules_fallback = which(rrule, Tuple{Any})

"""
  has_chain_rrule(T)

For a type-tuple `T` e.g. `Tuple{typeof(f), Int, Float64}`, checks if there is a `rrule` defined for it.
Excluding the generic fallback.
The first return value is `true` if the `rrule` exists, `false` otherwise.
If it does not, then the second argument is a list of edges to attach to the CodeInfo for a generated function,
such that if a suitable rule is defined later, the generated function will recompile.
"""
function has_chain_rrule(T)
  m = meta(Tuple{typeof(rrule),T.parameters...})
  if m.method !== chainrules_fallback
    # found a rrule, no need to add any edges
    return true, nothing
  end

  return false, m.instance
end

"""
    is_kwfunc(sigt...)

Determines if `sigt` is the type signature of a kwfunction.
Each element of `sigt` should be a type.
Either the first 3 types are a kwfunc type, a NamedTuple and the matching base function type,
or the first argument is the base function type and it is not a kwfunction.
the remaining types in `sigt` are the types of the argument.

"""
is_kwfunc(::Vararg{Any}) = false
is_kwfunc(k, ::Type{<:NamedTuple}, f, args...) = k===Core.kwftype(f)

"""
    wrap_chainrules_output(x)

Convert `x` from the differentials types ChainRules uses  to the format Zygote uses internally
(including conjugating complex gradients).
"""
@inline wrap_chainrules_output(x) = conj(unthunk(x))  # For now we are just not going to deal with thunks
@inline wrap_chainrules_output(x::Tuple) = map(wrap_chainrules_output, x)
@inline wrap_chainrules_output(x::ChainRules.AbstractZero) = nothing
for T_outer in (:Tuple, :NamedTuple)
  # we create separate methods rather than using a `Union` + an `if` so that we avoid a
  # branch that changes output type, because nested AD on that kinda thing makes Zygote less
  # than happy.
  @eval @inline function wrap_chainrules_output(x::ChainRules.Composite{P, T}) where {P, T<:$T_outer}
    xp = map(wrap_chainrules_output, x)
    convert($T_outer, xp)
  end
end

"""
    wrap_chainrules_input(x)

Convert `x` from the format  Zygote uses internally (including conjugated complex gradients)
to differentials types ChainRules uses.
"""
@inline wrap_chainrules_input(x) = conj(x)
@inline wrap_chainrules_input(::Nothing) = ChainRules.Zero()
@inline function wrap_chainrules_input(xs::Union{Tuple, NamedTuple})
  xp = map(wrap_chainrules_input, xs)
  ChainRules.Composite{Any, typeof(xp)}(xp)
end

"""
  ZBack{F}(back) <: Function

Wrapper for a ChainRules pullback `back`, that causes it to follow Zygote conventions.
(A functor here is used rather than a closure to avoid boxing issues);
"""
struct ZBack{F} <: Function
  back::F
end
@inline (s::ZBack)(dy) = wrap_chainrules_output(s.back(wrap_chainrules_input(dy)))
# `nothing->nothing` can be deleted after https://github.com/FluxML/Zygote.jl/issues/603
# though it might be worth keeping as a performance optimization (benchmarking pending)
@inline (s::ZBack)(::Nothing) = nothing

"""
    chain_rrule(f, args...)

Returns a the (primal) value of `f(args...)` and a pullback, by invoking `ChainRulesCore.rrule(f, args...)`.
The pullback is appropriately wrapped up to follow Zygote conventions.
"""
@inline function chain_rrule(f, args...)
  y, back = rrule(f, args...)
  return y, ZBack(back)
end

"""
  chain_rrule_kw(kwf, kwargs, f, args...)

As per [`chain_rrule`](@ref) but with support for kwargs.
`kwf` should be the kwfunc matching to `f`, and `kwargs` are a `NamedTuple` of keyword arguments.
"""
@inline function chain_rrule_kw(kwf, kwargs, f, args...)
  y, back = rrule(f, args...; kwargs...)
  kw_zpullback(dy) = (nothing, nothing, ZBack(back)(dy)...)  # first two nothings are for kwfunc and kwargs
  return y, kw_zpullback
end


text:
#############################################################################
# polyhedra.jl
# Handles constraints specifying the feasible region to be polyhedra
# associated with a submodular function f: 2ⱽ → ℝ, where V is the base set.
# The included polyhedra are:
#
## The submodular polyhedron: SubmodPoly = {s ∈ ℝᵖ, ∀ A ⊆ V, s(A) ≤ F(A)}.
## The base polyhedron: BasePoly = {s ∈ ℝᵖ, s(V) = F(V), ∀ A ⊆ V, s(A) ≤ F(A)}
##                               = SubmodPoly ∩ {s(V) = F(v)}.
## The positive polyhedron: PosPoly = {s ∈ ℝᵖ₊, ∀ A ⊆ V, s(A) ≤ F(A)}
##                                  = SubmodPoly ∩ ℝᵖ₊.
## The symmetric submodular polyhedron: SymPoly = {s ∈ ℝᵖ, ∀ A ⊆ V, |s|(A) ≤ F(A)}
##                                              = {s ∈ ℝᵖ, |s| ∈ SubmodPoly}.
# Here p = |V|, i.e. the cardinality of the base set.
#############################################################################

import Base.in

export AssocPoly
export SubmodPoly, BasePoly, PosPoly, SymPoly
export in, fenchel

abstract type AssocPoly <: ContiSet end

mutable struct SubmodPoly{T} <: AssocPoly
  head::Symbol
  F::T              # the submodular function
  V::AbstractArray # indexes of the base set
end

function SubmodPoly(F::SubmodFunc)
  @assert evaluate(F, [])[1] == 0
  x = get_sv(F)[1]
  V = x.baseset
  return(SubmodPoly(:subpoly, F, V))
end

mutable struct BasePoly{T} <: AssocPoly
  head::Symbol
  F::T
  V::AbstractArray
end

function BasePoly(F::SubmodFunc)
  @assert evaluate(F, [])[1] == 0
  x = get_sv(F)[1]
  V = x.baseset
  return(BasePoly(:basepoly, F, V))
end

mutable struct PosPoly{T} <: AssocPoly
  head::Symbol
  F::T
  V::AbstractArray
end

function PosPoly(F::SubmodFunc)
  @assert evaluate(F, [])[1] == 0
  x = get_sv(F)[1]
  V = x.baseset
  return(PosPoly(:pospoly, F, V))
end

mutable struct SymPoly{T} <: AssocPoly
  head::Symbol
  F::T
  V::AbstractArray
end

function SymPoly(F::SubmodFunc)
  @assert evaluate(F, [])[1] == 0
  x = get_sv(F)[1]
  V = x.baseset
  return(SymPoly(:sympoly, F, V))
end

function get_sv(p::AssocPoly)
  return get_sv(p.F)
end

function get_cv(p::AssocPoly)
  return Variable[]
end

# check if w is in p

function in(w::AbstractArray, p::AssocPoly, Tol::Float64 = 1e-3)
  n = length(w)
  @assert length(p.V) == n
  x = prox(p, w, Tol)
  if sum(abs.(x - w)) < Tol
    return true
  else
    return false
  end
end

# compute min w^x st x in p

# fenchel(p, w) = greedy(p.F, -w, p.V), w is nonpositve,
#               = -Inf,                 w is not nonpositive
function fenchel(p::SubmodPoly, w::AbstractArray)
  n = length(w)
  @assert length(p.V) == n
  if any(x -> x>0, w)
    return -Inf
  else
    return greedy(p.F, -w)
    # return greedy(p.F, -w, p.V)
  end
end

# fenchel(p, w) = greedy(p.F, -w, p.V)
function fenchel(p::BasePoly, w::AbstractArray)
  n = length(w)
  @assert length(p.V) == n
  return greedy(p.F, -w)
  # return greedy(p.F, -w, p.V)
end

# fenchel(p, w) = -sign(w_{-}) * greedy(p.F, -w_{-}, p.V)
function fenchel(p::PosPoly, w::AbstractArray)
  n = length(w)
  @assert length(p.V) == n
  u = 0.5 * (w - abs.(w))
  v = greedy(p.F, -u)
  # v = greedy(p.F, -u, p.V)
  v = -sign.(u) .* v
  return v
end

# fenchel(p, w) = -sign(w) .* greedy(p.F, |w|, p.V)
function fenchel(p::SymPoly, w::AbstractArray)
  n = length(w)
  @assert length(p.V) == n
  v = greedy(p.F, abs.(w))
  # v = greedy(p.F, abs.(w), p.V)
  v = -sign.(w) .* v
  return v
end


text:
module AdversarialGridworld
    using POMDPs
    using POMDPModelTools
    using StaticArrays
    using Parameters
    using Base.Cartesian
    using Random
    using Compose
    using ColorSchemes
    using LocalApproximationValueIteration
    using GridInterpolations
    using LocalFunctionApproximation
    using Serialization
    using POMDPSimulators
    using Statistics
    using Cairo, Fontconfig

    export GWPos, TwoAgentPos, aind, dir, syma, AdversarialGridworldMDP, DoubleAdversarialGridworldMDP, valid_pos,
            ThreeAgentPos, aind_ego, aind_adv, syma_ego, syma_adv,
            random_valid_pos, ego_pos, adversary_pos, agents_overlap
    include("mdp.jl")
    include("mdp_2adv.jl")

    export load_items, convert_policies_to_tables, solve_for_policies,
            generate_task_set, generate_task, solve_for_policy, value_to_table,
            action_values, policy_to_table, average_performance
    include("task_generator.jl")

end


text:
using CELFI
using Test

@testset "CELFI.jl" begin
    # Write your tests here.
end


text:
"""
  Constructs and returns the energy system model.
"""

function buildmodel(input)

    println("\nBuilding model...")

    @unpack REGION, PLANT, PLANTFACT, HOUR, numregions, CO2Cap, load, maxcap, assum, discountrate, cf, hydro_inflow = input

    m = Model(Gurobi.Optimizer)

    @variables m begin
        Electricity[r in REGION, p in PLANT, h in HOUR]       >= 0        # MW
        InstalledCapacity[r in REGION, p in PLANT]            >= 0        # MW
        ReservoirContent[h in HOUR]                           >= 0        # TWh
        Emission                                              >= 0        # Mton

        # # Exercise 2
        BatteryContent[r in REGION, h in HOUR]                >= 0        # MWh / MW
        BatteryInput[r in REGION, h in HOUR]                  >= 0        # MW
        #
        # # Exercise 3
        TransmissionCapacity[r in REGION, R in REGION]        >= 0        # MW
        Transmission[r in REGION, R in REGION, h in HOUR]     >= 0        # MW, r: transmission in, R: transmission out
    end #variables
    #print(Capacity)

    #Variable bounds
    for r in REGION, p in PLANT
        set_upper_bound(InstalledCapacity[r, p], maxcap[r, p]) # MW
    end

    function a(d,p)
        1-(1/(1+d)^assum[:Lifetime,p])
    end

    # test expression TODO: check all expressions
    @expression(m, Runningcost[r in REGION], sum(Electricity[r,p,h].*assum[:RunCost,p] for p in PLANT, h in HOUR)) # €
    @expression(m, Investmentcost[r in REGION],
        sum(InstalledCapacity[r,p].*(assum[:InvestmentCost,p].*(10^3)).*discountrate./a(discountrate,p) for p in PLANT)) # €
    @expression(m, Fuelcost[r in REGION], sum(Electricity[r,p,h].*assum[:FuelCost,p]./assum[:Efficiency,p] for h in HOUR, p in PLANT)) #€
    @expression(m, Systemcost[r in REGION], Runningcost[r] + Investmentcost[r] + Fuelcost[r]) # €

    @expression(m, GeneratedElectricity[r in REGION, h in HOUR], sum(Electricity[r,p,h] for p in [:Wind, :PV, :Gas, :Hydro]))#, :Nuclear])) # MW
    # @expression(m, GeneratedElectricityNet[r in REGION, h in HOUR], sum(Electricity[r,p,h] for p in PLANT)-BatteryInput[r,h])
    @expression(m, GeneratedElectricityNet[r in REGION, h in HOUR],
        sum(Electricity[r,p,h] for p in PLANT)-BatteryInput[r,h]-sum(Transmission[R, r, h] for R in REGION)./assum[:Efficiency,:Transmission]) # transmission input in Electricity[r, :Transmission,h]

    @expression(m, CapacityPerHour[r in REGION, p in PLANT, h in HOUR], sum(InstalledCapacity[r, p].*cf[r, p, h])) # MW

    @constraints m begin
        Generation[r in REGION, p in PLANT, h in HOUR], # name of constraint
            Electricity[r, p, h] <= CapacityPerHour[r, p, h]

        Syscost[r in REGION], # name of constraint
            Systemcost[r] >= 0 # sum of all annualized costs

        CapacityConstraints[r in REGION, p in PLANT],
            InstalledCapacity[r, p] <= maxcap[r,p]

        #GenElec[r in REGION, h in HOUR],
        #   GeneratedElectricity[r, h] >= load[r, h]
        GenElec[r in REGION, h in HOUR],
            GeneratedElectricityNet[r, h] >= load[r, h]

        ReservoirBalance[h in HOUR[1:end-1]],
            ReservoirContent[h+1] <= ReservoirContent[h] + hydro_inflow[h] - Electricity[:SE, :Hydro, h]./(10^6)

        ReservoirEndOfYear,
            ReservoirContent[1] == ReservoirContent[end]

        ReservoirLimit[h in HOUR],
            ReservoirContent[h] <= 33 # TWh

        # #Exercise 2:
        MinEmission,
            Emission >= sum((Electricity[r, :Gas, h]./assum[:Efficiency, :Gas]).*assum[:EmissionFactor, :Gas] for r in REGION, h in HOUR)./(10^6)

        MaxEmission,
            Emission <= CO2Cap #Mton CO2

        BatteryBalance[r in REGION, h in HOUR[1:end-1]],
            BatteryContent[r, h+1] <= BatteryContent[r, h] - Electricity[r, :Batteries, h]./assum[:Efficiency, :Batteries] + BatteryInput[r, h] #TODO: h-1 correct?

        BatteryInputMax[r in REGION, h in HOUR],
            BatteryInput[r, h] <= sum(Electricity[r,p,h] for p in PLANT) # can get input from batteries, but would be inefficient

        BatteryMax[r in REGION, h in HOUR],
            BatteryContent[r, h] <= InstalledCapacity[r, :Batteries]

        # #Exercise 3
        TransmissionCapacityEquality[r in REGION, R in REGION],
            TransmissionCapacity[r, R] == TransmissionCapacity[R, r]

        TransmissionCapacityMax[r in REGION, R in REGION],
            TransmissionCapacity[r, R] <= maxcap[r, :Transmission]

            # Kolla så att blev rätt
        TransmissionInstalledCapacity[r in REGION],
            sum(TransmissionCapacity[r, R] for R in REGION)./2 == InstalledCapacity[r, :Transmission]

        TransmissionMax[r in REGION, R in REGION, h in HOUR],
            Transmission[r, R, h] <= TransmissionCapacity[r, R]

        TransmissionProduction[r in REGION, h in HOUR], # transmission in
            sum(Transmission[r, R, h] for R in REGION) == Electricity[r, :Transmission, h]

    end #constraints

    @objective m Min begin
        sum(Systemcost[r] for r in REGION)
    end # objective

    return (; m, InstalledCapacity, Electricity, Emission, GeneratedElectricity, TransmissionCapacity)

end # buildmodel

function calculate_generation(input)

end # calculate_generation


text:
module BilevelBenchmark

if VERSION < v"0.7.0"
    Pkg.installed("BinDeps") == nothing && Pkg.add("BinDeps")
    Cvoid = Void
end

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("This test function suit is not properly installed. Please run Pkg.build(\"BilevelBenchmark\")")
end

export bilevel_leader, bilevel_follower, bilevel_solutions
export bilevel_settings, bilevel_ranges
export TP_settings,TP_leader,TP_follower,TP_test, TP_optimum
export SMD_settings, SMD_solutions, SMD_ranges, SMD_leader, SMD_follower, SMD_optimum
export PMM_settings, PMM_Ψ, PMM_leader, PMM_follower, PMM_test, PMM_ranges, TP_ranges, PMM_optimum
export SMD_Ψ
export get_problem

include("TP.jl")
include("SMD.jl")
include("PMM.jl")
include("deprecated.jl")

function get_problem(benchmark, fnum; D_ul = 2, D_ll = 3)
    if benchmark == :SMD
        F = (x,y) -> SMD_leader(x,y,fnum)
        f = (x,y) -> SMD_follower(x,y,fnum)
        Ψ = x -> SMD_Ψ(x, D_ll, fnum)
        bounds_ul, bounds_ll = SMD_ranges(D_ul, D_ll, fnum)
        z = SMD_optimum(fnum)
    elseif benchmark == :PMM
        F = (x,y) -> PMM_leader(x,y,fnum)
        f = (x,y) -> PMM_follower(x,y,fnum)
        Ψ = x -> PMM_Ψ(x, D_ll, fnum)
        bounds_ul, bounds_ll = PMM_ranges(D_ul, D_ll, fnum)
        z = PMM_optimum(fnum)
    elseif benchmark == :TP
        @warn "TP not implement the Ψ mapping"
        F = (x,y) -> TP_leader(x,y,fnum)
        f = (x,y) -> TP_follower(x,y,fnum)
        Ψ = x -> zeros(D_ll)
        bounds_ul, bounds_ll = TP_ranges(fnum)
        z = TP_optimum(fnum)
    else
        error("Only :SMD, :PMM and :TP are defined.")
        return
    end

    return F, f, bounds_ul, bounds_ll, Ψ, z
end

end # module


text:
module M
    using CxxWrap
    using CxxWrap.CxxWrapCore: SafeCFunction

    struct JLClient
        fuck::Int32
        # fuckCb::Ptr{Nothing}
        fuckCb::SafeCFunction
    end

    function fuckCb(x)
        println(x)
        println(x)
        println(x)
        return nothing
    end

    function Client()
        jl_client = JLClient(
            32,
            # @cfunction(fuckCb, Nothing, (Int32, )),
        )
        return TestCppClient(jl_client)
    end

    @wrapmodule(joinpath(@__DIR__, "../../result/lib", "libIBApiWrapper"))

    function __init__()
        @initcxx
    end
end


text:
using DSP
function ExtractFlattedGroundRoll(AShot,datainfo)
    ### Input:
	### AShot: one shot gather (nt,nrx,nry)
	### datainfo: a Dictionary contains info for one shot gather
	### Output:
	### A shot gather with only the flatted ground roll estimatation
	nrx = datainfo["nrx"]; #number of receiver x
	nry = datainfo["nry"]; #number of receiver y
	Rgrid = datainfo["Rgrid"]; #receiver grid space
	dt = datainfo["dt"]; # number of time
	# Pad zeros in time domain for the shot data
	AShot = [zeros(300,nrx,nry);AShot;zeros(100,nrx,nry)]

	vel = 1415; #velocity sand(dry) 0.2-1 Km/s, sand(water saturated 1.5-2 Km/s)

	# find the center position of receiver = source position
	CentX = (nrx+1)/2;
	CentY = (nry+1)/2;

	# calculate the shift grid for each receiver and save them in the matrix shiftT
	a = 1:nrx;
	b = 1:nry;
	row = reshape(repeat(a,outer=nry),nrx,nry);
	col = reshape(repeat(b,inner=nry),nrx,nry);
	shiftT = round.(sqrt.((row .-CentX).^2 .+ (col .-CentY).^2) .*(Rgrid/vel/(dt*1e-3)));

	Shift_Data = zeros(size(AShot));
	nt1 =  datainfo["nt"] + 300 + 100;
	for i = 1:nrx
	    for j = 1:nry
		Shift_Data[1:Int(nt1-shiftT[i,j]),i,j] = AShot[Int(shiftT[i,j]+1):nt1,i,j];
	    end
	end

	# define a filter in time domain to mute the signal besides groundroll
	A = zeros(nt1,1);
	band = 10;
	low = 230;
	high = 540;

	A[low+band+1:high-band-1] .= DSP.Windows.tukey(high-band-1-(low+band+1)+1,0.6;padding=0,zerophase=false)

	Shift_taper_Data = zeros(size(Shift_Data));
	for i = 1:nrx
	    for j = 1:nry
		Shift_taper_Data[:,i,j] = Shift_Data[:,i,j] .* A;
	    end
	end
	return Shift_taper_Data;
end


text:
# ---
# title: 146. LRU Cache
# id: problem146
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design
# link: <https://leetcode.com/problems/lru-cache/description/>
# hidden: true
# ---
#
# Design a data structure that follows the constraints of a **[Least Recently
# Used (LRU)
# cache](https://en.wikipedia.org/wiki/Cache_replacement_policies#LRU)**.
#
# Implement the `LRUCache` class:
#
#   * `LRUCache(int capacity)` Initialize the LRU cache with **positive** size `capacity`.
#   * `int get(int key)` Return the value of the `key` if the key exists, otherwise return `-1`.
#   * `void put(int key, int value)` Update the value of the `key` if the `key` exists. Otherwise, add the `key-value` pair to the cache. If the number of keys exceeds the `capacity` from this operation, **evict** the least recently used key.
#
# **Follow up:**
# Could you do `get` and `put` in `O(1)` time complexity?
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["LRUCache", "put", "put", "get", "put", "get", "put", "get", "get", "get"]
#     [[2], [1, 1], [2, 2], [1], [3, 3], [2], [4, 4], [1], [3], [4]]
#     **Output**
#     [null, null, null, 1, null, -1, null, -1, 3, 4]
#
#     **Explanation**
#     LRUCache lRUCache = new LRUCache(2);
#     lRUCache.put(1, 1); // cache is {1=1}
#     lRUCache.put(2, 2); // cache is {1=1, 2=2}
#     lRUCache.get(1);    // return 1
#     lRUCache.put(3, 3); // LRU key was 2, evicts key 2, cache is {1=1, 3=3}
#     lRUCache.get(2);    // returns -1 (not found)
#     lRUCache.put(4, 4); // LRU key was 1, evicts key 1, cache is {4=4, 3=3}
#     lRUCache.get(1);    // return -1 (not found)
#     lRUCache.get(3);    // return 3
#     lRUCache.get(4);    // return 4
#
#
#
#
# **Constraints:**
#
#   * `1 <= capacity <= 3000`
#   * `0 <= key <= 3000`
#   * `0 <= value <= 104`
#   * At most `3 * 104` calls will be made to `get` and `put`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
greeting(name) = "Hello, $name"

bye(name) = "Bye, $(name)!"

max_stars_repo_name: JuliaTagBot/Cincinnati

max_issues_repo_name: JuliaTagBot/Cincinnati

max_forks_repo_name: JuliaTagBot/Cincinnati


text:
# ------------ Value alias and APIs ------------ #

const Value = MLIR.API.MlirValue

get_type(v::Value) = MLIR.API.mlirValueGetType(v)
get_owner_block(v::Value) = MLIR.API.mlirBlockArgumentGetOwner(v)
get_owner_block_arg_index(v::Value) = MLIR.API.mlirBlockArgumentGetArgNumber(v)
get_owner_op(v::Value) = MLIR.API.mlirOpResultGetOwner(v)
get_owner_op_result_index(v::Value) = MLIR.API.mlirOpResultGetResultNumber(v)
is_null(v::Value) = MLIR.API.mlirValueIsNull(v)
Base.:(==)(v1::Value, v2::Value) = MLIR.API.mlirValueEqual(v1, v2)
is_block_args(v::Value) = MLIR.API.mlirValueIsABlockArgument(v)
is_op_result(v::Value) = MLIR.API.mlirValueIsAOpResult(v)
dump(v::Value) = MLIR.API.mlirValueDump(v)
set_type!(v::Value, t::Type) = MLIR.API.mlirBlockArgumentSetType(v, t)

@doc(
"""
const Value = MLIR.API.MlirValue
""", Value)


text:
using TopologicalDataAnalysis
using Test

@testset "TopologicalDataAnalysis.jl" begin
    # Write your tests here.
end


text:
function cholesky_mod(beta, H)
    @assert beta > 0
    if minimum(diag(H)) > 0
        tau = 0
    else
        tau = -minimum(diag(H)) + beta
    end
    while true
        candidate = H + tau * eye(H)
        check_posdef(candidate) && return candidate
        tau = max(2*tau, beta)
    end
end

function check_posdef(A)
    try
        L = chol(A)
        return true
    catch
        return false
    end
end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# A* shortest-path algorithm

function a_star_impl!(
        g::AGraphOrDiGraph,# the graph
        s::Integer,
        t::Integer, # the end vertex
        colormap::Vector{Int},  # an (initialized) color-map to indicate status of vertices
        distmx::AEdgeMap{T},
        heuristic::Function    # heuristic fn (under)estimating distance to target
    ) where T<:Number
    V = vertextype(g)
    colormap[s] = 1
    frontier = PriorityQueue{Tuple{T,Vector{Edge{V}},Int},T}()
    frontier[(zero(T), Vector{Edge{V}}(), s)] = zero(T)

    path = Vector{Edge{V}}()
    foundpath = false
    while !isempty(frontier)
        (cost_so_far, path, u) = dequeue!(frontier)
        if u == t
            foundpath = true
            break
        end
        for e in out_edges(g, u)
            v = dst(e)
            if colormap[v] < 2
                dist = distmx[e]

                colormap[v] = 1
                new_path = cat(path, Edge{V}(u,v),dims=1)
                path_cost = cost_so_far + dist
                enqueue!(frontier,
                        (path_cost, new_path, v),
                        path_cost + heuristic(v))
            end
        end
        colormap[u] = 2
    end
    return foundpath ? path : Vector{Edge{V}}()
end

"""
    a_star(g, s, t, distmx=weights(g), heuristic = n->0)

Computes the shortest path between vertices `s` and `t` using the
[A* search algorithm](http://en.wikipedia.org/wiki/A%2A_search_algorithm). An
optional heuristic function and edge distance matrix may be supplied.
Returns an empty path if there are no such paths.
"""
function a_star(
        g::AGraphOrDiGraph,  # the graph
        s::Integer,                       # the start vertex
        t::Integer,                       # the end vertex
        distmx::AEdgeMap=weights(g),
        heuristic::Function = n -> 0
    )
            # heuristic (under)estimating distance to target
    colormap = zeros(Int, nv(g))
    a_star_impl!(g, s, t, colormap, distmx, heuristic)
end


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

mutable struct IoK8sApiRbacV1AggregationRule <: SwaggerModel
    clusterRoleSelectors::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1LabelSelector} } # spec name: clusterRoleSelectors

    function IoK8sApiRbacV1AggregationRule(;clusterRoleSelectors=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1AggregationRule, Symbol("clusterRoleSelectors"), clusterRoleSelectors)
        setfield!(o, Symbol("clusterRoleSelectors"), clusterRoleSelectors)
        o
    end
end # type IoK8sApiRbacV1AggregationRule

const _property_map_IoK8sApiRbacV1AggregationRule = Dict{Symbol,Symbol}(Symbol("clusterRoleSelectors")=>Symbol("clusterRoleSelectors"))
const _property_types_IoK8sApiRbacV1AggregationRule = Dict{Symbol,String}(Symbol("clusterRoleSelectors")=>"Vector{IoK8sApimachineryPkgApisMetaV1LabelSelector}")
Base.propertynames(::Type{ IoK8sApiRbacV1AggregationRule }) = collect(keys(_property_map_IoK8sApiRbacV1AggregationRule))
Swagger.property_type(::Type{ IoK8sApiRbacV1AggregationRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1AggregationRule[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1AggregationRule }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1AggregationRule[property_name]

function check_required(o::IoK8sApiRbacV1AggregationRule)
    true
end

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


text:
module Chemistry

using JuLIP.FIO: load_dict

export atomic_number,
       chemical_symbol,
       atomic_mass,
       symmetry,
       lattice_constant,
       rnn

data = load_dict(@__DIR__()[1:end-3] * "data/asedata.json")

const _rnn = [Float64(d) for d in data["rnn"]]
const _symbols = [Symbol(s) for s in data["symbols"]]
const _masses = [Float64(m) for m in data["masses"]]
const _refstates = Dict{String, Any}[ d == nothing ? Dict{String, Any}() : d
                                      for d in data["refstates"]]

const _numbers = Dict{Symbol, Int16}()
for (n, sym) in enumerate(_symbols)
   _numbers[sym] = Int16(n - 1)
end

atomic_number(sym::Symbol) = _numbers[sym]

chemical_symbol(z::Integer) = _symbols[z+1]

atomic_mass(z::Integer) = _masses[z+1]
atomic_mass(sym::Symbol) = atomic_mass(atomic_number(sym))

element_name(z::Integer) = _names[z+1]
element_name(sym::Symbol) = element_name(atomic_number(sym))

rnn(z::Integer) = _rnn[z+1]       # >>> TODO: move to utils ?? bulk as well??
rnn(sym::Symbol) = rnn(atomic_number(sym))

symmetry(z::Integer) = Symbol(_refstates[z+1]["symmetry"])
symmetry(sym::Symbol) = symmetry(atomic_number(sym))

lattice_constant(z::Integer) = Float64( _refstates[z+1]["a"] )
lattice_constant(sym::Symbol) = lattice_constant(atomic_number(sym))

refstate(z::Integer) = _refstates[z+1]
refstate(sym::Symbol) = refstate(atomic_number(sym))

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -19.736879506320445 <= q <= 22.217442927624916)

                     add_NL_constraint(m, :(gelu(-0.3136061325511519 + 0.2821341474705279*gelu(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])) + 0.4396178028948987*gelu(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])) + -0.24849525966922226*gelu(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])) + 0.5972646324288218*gelu(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])) + -0.8620507270647968*gelu(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))) + gelu(-0.9705953045797062 + -0.6259685253929317*gelu(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])) + 0.41918630672941015*gelu(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])) + -0.9388772793667424*gelu(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])) + -0.18849282727638839*gelu(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])) + -0.7819718693068158*gelu(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))) + gelu(0.8285261526161638 + 0.6863822914423507*gelu(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])) + -0.7060930690780838*gelu(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])) + -0.2956343589188948*gelu(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])) + 0.32353098534638214*gelu(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])) + 0.6963022628505415*gelu(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))) + gelu(-0.2246585184074279 + 0.8591320330102441*gelu(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])) + -0.5037522661683527*gelu(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])) + 0.9621050832972355*gelu(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])) + -0.1995700954953512*gelu(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])) + 0.442251462608958*gelu(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))) + gelu(0.5761249271362421 + -0.042724357144423575*gelu(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])) + -0.4537006876708345*gelu(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])) + -0.12519521024691826*gelu(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])) + -0.7426009177864223*gelu(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])) + -0.8758186284143128*gelu(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function compare_to(obj::Enum, arg0::Enum)
    return jcall(obj, "compareTo", jint, (Enum,), arg0)
end

function describe_constable(obj::Enum)
    return jcall(obj, "describeConstable", Optional, ())
end

function equals(obj::Enum, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_declaring_class(obj::Enum)
    return jcall(obj, "getDeclaringClass", Class, ())
end

function get_m_hz_frequency(obj::Frequency)
    return jcall(obj, "getMHzFrequency", jdouble, ())
end

function get_name(obj::Frequency)
    return jcall(obj, "getName", JString, ())
end

function get_ratio(obj::Frequency)
    return jcall(obj, "getRatio", jdouble, ())
end

function get_satellite_system(obj::Frequency)
    return jcall(obj, "getSatelliteSystem", SatelliteSystem, ())
end

function get_wavelength(obj::Frequency)
    return jcall(obj, "getWavelength", jdouble, ())
end

function hash_code(obj::Enum)
    return jcall(obj, "hashCode", jint, ())
end

function name(obj::Enum)
    return jcall(obj, "name", JString, ())
end

function ordinal(obj::Enum)
    return jcall(obj, "ordinal", jint, ())
end

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

function value_of(::Type{Enum}, arg0::Class, arg1::JString)
    return jcall(Enum, "valueOf", Enum, (Class, JString), arg0, arg1)
end

function value_of(::Type{Frequency}, arg0::JString)
    return jcall(Frequency, "valueOf", Frequency, (JString,), arg0)
end

function values(::Type{Frequency})
    return jcall(Frequency, "values", Vector{Frequency}, ())
end


text:
# Load dependencies
using NeuralQuantum, QuantumLattices
using Logging, Printf, ValueHistories

# Select the numerical precision
T      = Float64
# Select how many sites you want
Nsites = 7

# Create the lattice as [Nx, Ny, Nz]
lattice = SquareLattice([Nsites],PBC=true)
# Create the lindbladian for the QI model
lind = quantum_ising_lind(lattice, g=1.0, V=2.0, γ=1.0)
# Create the Problem (cost function) for the given lindbladian
# targeting the Steady State, using a memory efficient encoding and
# minimizing |Lρ|^2 as a variance, which is more efficient.
prob = SteadyStateProblem(T, lind, operators=true, variance=true)

#-- Observables
# Define the local observables to look at.
Sx  = QuantumLattices.LocalObservable(lind, sigmax, Nsites)
Sy  = QuantumLattices.LocalObservable(lind, sigmay, Nsites)
Sz  = QuantumLattices.LocalObservable(lind, sigmaz, Nsites)
# Create the problem object with all the observables to be computed.
oprob = ObservablesProblem(Sx, Sy, Sz)

# Define the Neural Network. A NDM with N visible spins and αa=2 and αh=1
#alternative vectorized rbm: net  = RBMSplit(Complex{T}, Nsites, 6)
net  = NDM(T, Nsites, 1, 2)
# Create a cached version of the neural network for improved performance.
cnet = cached(net)
# Chose a sampler. Options are FullSumSampler() which sums over the whole space
# ExactSampler() which does exact sampling or MCMCSamler which does a markov
# chain.
sampl = MCMCSampler(Metropolis(Nsites), 3000, burn=50)
# Chose a sampler for the observables.
osampl = FullSumSampler()

# Chose the gradient descent algorithm (alternative: Gradient())
# for more information on options type ?SR
algo  = SR(ϵ=T(0.01), use_iterative=true)
# Optimizer: how big the steps of the descent should be
optimizer = Optimisers.Descent(0.02)

# Create a multithreaded Iterative Sampler.
is = MTIterativeSampler(cnet, sampl, prob, algo)
ois = MTIterativeSampler(cnet, osampl, oprob)

# Create the structure to store all output data
minimization_data = MVHistory()
Δw = grad_cache(cnet)

# Solve iteratively the problem
for i=1:110
    # Sample the gradient
    grad_data  = sample!(is)
    obs_data = sample!(ois)

    # Logging
    @printf "%4i -> %+2.8f %+2.2fi --\t \t-- %+2.5f\n" i real(grad_data.L) imag(grad_data.L) real(obs_data.ObsAve[1])
    push!(minimization_data, :loss, grad_data.L)
    for (name,val)=zip(obs_data.ObsNames, obs_data.ObsAve)
        push!(minimization_data, Symbol(name), val)
    end

    succ = precondition!(Δw.tuple_all_weights, algo , grad_data, i)
    !succ && break
    Optimisers.update!(optimizer, cnet, Δw)
end

# Optional: compute the exact solution
ρ   = last(steadystate.master(lind)[2])
ESx = real(expect(SparseOperator(Sx), ρ))
ESy = real(expect(SparseOperator(Sy), ρ))
ESz = real(expect(SparseOperator(Sz), ρ))
exacts = Dict("Sx"=>ESx, "Sy"=>ESy, "Sz"=>ESz)
## - end Optional

using Plots
data = minimization_data

iter_cost, cost = get(minimization_data[:loss])
pl1 = plot(iter_cost, real(cost), yscale=:log10)

iter_mx, mx = get(minimization_data[:obs_1])
pl2 = plot(iter_mx, real(mx))
hline!(pl2, [ESx,ESx]);

plot(pl1, pl2, layout=(2,1))


text:
"""
    function generate(
        saved_files_dir::String,
        schema_paths::Vector{String};
        generate_types::Bool = true,
        generate_functions::Bool = true,
    )

Generate Julia code files for GraphQL types and functions.

"graphqlgen_types.jl": contains all the GraphQL types

"graphqlgen_functions.jl": contains all the GraphQL functions (mutations, queries, subscriptions)

* `saved_files_dir`: directory where the generated files will be saved
* `schema_paths`: list of paths to GraphQL schema files. This can be a file or a directory. If it's a directory, it will be recursively searched for GraphQL schema files.
* `generate_types`: whether to generate "types.jl"
* `generate_functions`: whether to generate "functons.jl"
"""
function generate(
    saved_files_dir::String,
    schema_paths::Vector{String};
    generate_types::Bool = true,
    generate_functions::Bool = true,
)
    schema = ""
    for p in schema_paths
        if isfile(p)
            _, ext = splitext(p)
            if ext in (".graphql", ".schema")
                @info "reading in schema file" file = p
                str = String(read(fp))
                schema *= str
                schema *= "\n"
            end
        else
            for (root, dirs, files) in walkdir(p)
                for f in files
                    fp = joinpath(root, f)
                    ".git" in split(fp, Base.Filesystem.path_separator) && continue
                    _, ext = splitext(fp)
                    if ext in (".graphql", ".schema")
                        @info "reading in schema file" file = f
                        str = String(read(fp))
                        schema *= str
                        schema *= "\n"
                    end
                end
            end
        end
    end

    # generate types and functions
    types, functions = GraphQLGen.tojl(GraphQLGen.parse(schema))

    dir = abspath(saved_files_dir)
    try
        mkdir(dir)
    catch e
        if e isa Base.IOError && e.code == -17
            # ignore
        else
            rethrow()
        end
    end

    if generate_types
        filename = joinpath(dir, "graphqlgen_types.jl")
        open("$filename", "w") do f
            GraphQLGen.print(f, types)
        end
        @info "Generated Julia GraphQL types" path = filename
    end

    if generate_functions
        filename = joinpath(dir, "graphqlgen_functions.jl")
        open("$filename", "w") do f
            GraphQLGen.print(f, functions)
        end
        @info "Generated Julia GraphQL functions" path = filename
    end

    return nothing
end

function generate(
    saved_files_dir::String,
    schema_path::String;
    generate_types::Bool = true,
    generate_functions::Bool = true,
)
    generate(saved_files_dir, [schema_path]; generate_types, generate_functions)
end


text:
"""
    T_transform(Yc, Xc, Xp, x) (stationnary)

"""
function T_transform(Yc::ContinuousUnivariateDistribution,  # Distribution des données observées pour la période de calibration
                     Xc::ContinuousUnivariateDistribution,  # Distribution des données du modèle pour la période de calibration
                     Xp::ContinuousUnivariateDistribution,  # Distribution des données du modèle pour la période projetée
                     x::Real)

    F_Xp = cdf(Xp, x)  # Évaluation en x de la CDF des données du modèle pour la période projeté

    Fm1_Xc = quantile(Xc, F_Xp)  # Quantile correspondant de la distribution des données du modèle pour la période de calibration

    F_Yp = cdf(Yc,  Fm1_Xc)  # Évaluation de la CDF des données observées pour la période de calibration au point correspondant

    return F_Yp  # Valeur de la CDF évaluée au point x
end

"""
    T_transform(Yx, Xc, Xc_cov, Xp, Xp_cov, x)  (non-stationnary)

"""
function T_transform(Yc::EGPpower,  # Distribution des données observées pour la période de calibration
                     Xc::nonstatEGPpower,  # Distribution des données du modèle pour la période de calibration
                     Xc_cov::Real,  # covariable pour la période de calibration (ex : [CO₂] de l'année 1980)
                     Xp::nonstatEGPpower,  # Distribution des données du modèle pour la période projetée
                     Xp_cov::Real,  # covariable pour la période projetée (ex : [CO₂] de l'année 2050)
                     x::Real)

    F_Xp = cdf(Xp, x, Xp_cov)  # Évaluation en x de la CDF des données du modèle pour la période projeté

    Fm1_Xc = quantile(Xc, F_Xp, Xc_cov)  # Quantile correspondant de la distribution des données du modèle pour la période de calibration

    F_Yp = cdf(Yc,  Fm1_Xc)  # Évaluation de la CDF des données observées pour la période de calibration au point correspondant

    return F_Yp  # Valeur de la CDF évaluée au point x
end

## Bias correction

"""
    bias_correct(obs, ref, covariate_ref; fut, covariate_fut, censoring_obs, censoring_sim, wet_thresh)

"""
function bias_correct(obs::Array{Float64, 1},
                ref::Array{Float64, 1},
                covariate_ref::Array{Float64, 1};
                fut::Array{Float64, 1}=Float64[],
                covariate_fut::Array{Float64, 1}=Float64[],
                censoring_obs::Real=5.0,
                censoring_sim::Real=5.0,
                wet_thresh::Float64=1.0)

    # Transformation des valeurs en fonction du seuil
    obs_t = obs[obs .> wet_thresh] .- wet_thresh

    ref_t = ref[ref .> wet_thresh] .- wet_thresh
    covariate_ref_t = covariate_ref[ref .> wet_thresh]

    fut_t = fut[fut .> wet_thresh] .- wet_thresh
    covariate_fut_t = covariate_fut[fut .> wet_thresh]

    # Ajustement du modèle sur les observations
    fm_obs = EGPpowerfit(obs_t, censoring=censoring_obs)

    # Fusion des prédictions et des projections pour l'ajustement
    sim = vcat(ref_t, fut_t)
    covariate_sim = vcat(covariate_ref_t, covariate_fut_t)

    # Ajustement du modèle sur les simulations
    fm_sim, fms = EGPpowerfit(sim, censoring=censoring_sim, covariate=covariate_sim, scalecov=false, uppertailcov=true, lowertailcov=false)

    # Correction des prédictions
    output_ref = similar(ref)
    output_ref[ref .<= wet_thresh] .= 0.0
    output_ref[ref .> wet_thresh] = quantile.(fm_obs, cdf.(EGPpower.(fm_sim, covariate_ref_t), ref_t)) .+ wet_thresh

    if !isempty(fut)
        # Récupération de la CDF des obs pour le futur (Kallache)
        Yp = CDFT.(Ref(fm_obs), Ref(EGPpower(fm_sim, last(covariate_ref))), EGPpower.(fm_sim, covariate_fut_t))
        p = cdf.(getfield.(Yp, :Xp), fut_t)

        # Correction des projections
        output_fut = similar(fut)
        output_fut[fut .<= wet_thresh] .= 0.0
        output_fut[fut .> wet_thresh] = quantile.(Yp, p) .+ wet_thresh
        append!(output_ref, output_fut)

    end
    return output_ref
end

function bias_correct(obs::Array{Float64, 1},
                ref::Array{Float64, 1},
                covariate_ref::Array{Float64, 1},
                fut::Array{Float64, 1},
                covariate_fut::Array{Float64, 1};
                censoring_obs::Real=5.0,
                censoring_sim::Real=5.0,
                wet_thresh::Float64=1.0)

    return bias_correct(obs, ref, covariate_ref, fut=fut, covariate_fut=covariate_fut, censoring_obs=censoring_obs, censoring_sim=censoring_sim, wet_thresh=wet_thresh)
end

function bias_correct(obs::TimeArray{Float64,1,T,Array{Float64,1}},
                ref::TimeArray{Float64,1,T,Array{Float64,1}},
                covariate_ref::TimeArray{Float64,1,T,Array{Float64,1}},
                fut::TimeArray{Float64,1,T,Array{Float64,1}},
                covariate_fut::TimeArray{Float64,1,T,Array{Float64,1}};
                censoring_obs::Real=5.0,
                censoring_sim::Real=5.0,
                wet_thresh::Float64=1.0) where {T <: TimeType}

    x = bias_correct(values(obs), values(ref), values(covariate_ref), fut=values(fut), covariate_fut=values(covariate_fut), censoring_obs=censoring_obs,
      censoring_sim=censoring_sim, wet_thresh=wet_thresh)
    t = vcat(timestamp(ref), timestamp(fut))

    return TimeArray(t, x)
end

function bias_correct(obs::TimeArray{Float64,1,T,Array{Float64,1}},
                ref::TimeArray{Float64,1,T,Array{Float64,1}},
                covariate_ref::TimeArray{Float64,1,T,Array{Float64,1}};
                censoring_obs::Real=5.0,
                censoring_sim::Real=5.0,
                wet_thresh::Float64=1.0) where {T <: TimeType}

    x = bias_correct(values(obs), values(ref), values(covariate_ref), censoring_obs=censoring_obs, censoring_sim=censoring_sim, wet_thresh=wet_thresh)
    t = timestamp(ref)

    return TimeArray(t, x)
end


text:
@testset "Fminbox" begin
    # Quadratic objective function
    # For (A*x-b)^2/2
    function quadratic!(g, x, AtA, Atb, tmp)
        mul!(tmp, AtA, x)
        v = dot(x,tmp)/2 + dot(Atb,x)
        if g !== nothing
            g .= tmp .+ Atb
        end
        return v
    end

    Random.seed!(1)
    N = 8
    boxl = 2.0
    outbox = false
    # Generate a problem where the bounds-free solution lies outside of the chosen box
    local _objective
    while !outbox
        A = randn(N,N)
        AtA = A'*A
        b = randn(N)
        initial_x = randn(N)
        tmp = similar(initial_x)
        funcg! = (g, x) -> quadratic!(g, x, AtA, A'*b, tmp)
        _objective = OnceDifferentiable(x->funcg!(nothing, x), (g, x)->funcg!(g, x), funcg!, initial_x)
        results = optimize(_objective, initial_x, ConjugateGradient())
        results = optimize(_objective, Optim.minimizer(results), ConjugateGradient())  # restart to ensure high-precision convergence
        @test Optim.converged(results)
        opt_x = Optim.minimizer(results)
        g = similar(opt_x)
        @test funcg!(g, opt_x) + dot(b,b)/2 < 1e-8
        @test norm(g) < 1e-4
        outbox = any(t -> abs(t) > boxl, opt_x)
    end

    # fminbox
    l = fill(-boxl, N)
    u = fill(boxl, N)
    initial_x = (rand(N) .- 0.5) .* boxl
    for _optimizer in (ConjugateGradient(), GradientDescent(), LBFGS(), BFGS(), NGMRES(), OACCEL())
        debug_printing && printstyled("Solver: ", summary(_optimizer), "\n", color=:green)
        results = optimize(_objective, l, u, initial_x, Fminbox(_optimizer))
        if _optimizer isa NGMRES && Sys.iswindows() && Sys.WORD_SIZE == 64
            @test_broken Optim.converged(results)
        else
            @test Optim.converged(results)
        end
        @test summary(results) == "Fminbox with $(summary(_optimizer))"
        opt_x = Optim.minimizer(results)
        NLSolversBase.gradient!(_objective, opt_x)
        g = NLSolversBase.gradient(_objective)
        # check first-order constrained optimality conditions
        for i = 1:N
            @test abs(g[i]) < 3e-3 || (opt_x[i] < -boxl+1e-3 && g[i] > 0) || (opt_x[i] > boxl-1e-3 && g[i] < 0)
        end
    end

    # Throw ArgumentError when initial guess is outside the box
    @test_throws ArgumentError optimize(_objective, l, u, 2u, Fminbox(GradientDescent()))

    # tests for #180
    results = optimize(_objective, l, u, initial_x, Fminbox(), Optim.Options(outer_iterations = 2))
    @test Optim.iterations(results) == 2
    @test Optim.minimum(results) == _objective.f(Optim.minimizer(results))

    # Warn when initial condition is not in the interior of the box
    initial_x = rand([-1,1],N)*boxl
    @test_logs (:warn, "Initial position cannot be on the boundary of the box. Moving elements to the interior.\nElement indices affected: [1, 2, 3, 4, 5, 6, 7, 8]") optimize(_objective, l, u, initial_x, Fminbox(), Optim.Options(outer_iterations = 1))

    # might fail if changes are made to Optim.jl
    # TODO: come up with a better test
    #results = Optim.optimize(_objective, initial_x, l, u, Fminbox(); optimizer_o = Optim.Options(iterations = 2))
    #@test Optim.iterations(results) == 470
    @testset "simple input" begin
        function exponential(x)
            return exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
        end

        function exponential_gradient!(storage, x)
            storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
            storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
            storage
        end

        function exponential_gradient(x)
            storage = similar(x)
            storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
            storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
            storage
        end

        initial_x = [0.0, 0.0]
        optimize(exponential, exponential_gradient!, initial_x, BFGS())
        lb = fill(-0.1, 2)
        ub = fill(1.1, 2)
        od = OnceDifferentiable(exponential, initial_x)
        optimize(od, lb, ub, initial_x, Fminbox())
        od_forward = OnceDifferentiable(exponential, initial_x; autodiff = :forward)
        optimize(od_forward, lb, ub, initial_x, Fminbox())
        optimize(exponential, lb, ub, initial_x, Fminbox())
        optimize(exponential, exponential_gradient!, lb, ub, initial_x, Fminbox())
        optimize(od, lb, ub, initial_x)
        optimize(od_forward, lb, ub, initial_x)
        optimize(exponential, lb, ub, initial_x)
        optimize(exponential, exponential_gradient!, lb, ub, initial_x)
        @testset "inplace and autodiff keywords #616" begin
            optimize(exponential, lb, ub, initial_x, Fminbox())
            optimize(exponential, lb, ub, initial_x, Fminbox(); autodiff = :finite)
            optimize(exponential, lb, ub, initial_x, Fminbox(); autodiff = :forward)
            optimize(exponential, exponential_gradient, lb, ub, initial_x, Fminbox(), inplace = false)
        end
        @testset "error for second order methods #616" begin
            @test_throws ArgumentError optimize(x->x, (G,x)->x, rand(1),rand(1),rand(1), Fminbox(Newton()))
            @test_throws ArgumentError optimize(x->x, (G,x)->x, rand(1),rand(1),rand(1), Fminbox(NewtonTrustRegion()))
        end
        @testset "allow for an Optim.Options to be passed #623" begin
            optimize(exponential, lb, ub, initial_x, Fminbox(), Optim.Options())
            optimize(exponential, exponential_gradient!, lb, ub, initial_x, Fminbox(), Optim.Options())
            @test_broken optimize(exponential, exponential_gradient, lb, ub, initial_x, Optim.Options())
        end
    end
end


text:
"""
    (uvels, vvels, wvels) = streamfunction_xz()

Define a two-dimensional flow field in the x-z directions from a simple streamfunction.
"""
function streamfunction_xz()
    scal = 2e-1
    f(x, y, z) = scal*sin(x*2π/128)*sin(z*π/128) #stream function

    ϕcorners=[f(x,0.,z) for x in 0:128, z in -128:0]
    ϕcenters=[f(x,0.,z) for x in 0.5:128, z in -128:-0.5]

    uu=-diff(ϕcorners,dims=2)[1:end-1,:]
    ww=diff(ϕcorners,dims=1)
    uu=reshape(uu,(128,1,128))
    ww=reshape(ww,(128,1,129))
    uu = reverse(uu, dims=3)
    ww = reverse(ww, dims=3)

    uvels = fill(uu, 2)
    vvels = fill(0*uu, 2)
    wvels = fill(ww, 2)
    uvels = cat(uvels..., dims=4)
    vvels = cat(vvels..., dims=4)
    wvels = cat(wvels..., dims=4)

    return uvels, vvels, wvels, ϕcenters
end

"""
    (uvels, vvels, wvels) = streamfunction_xy()

Define a two-dimensional flow field in the x-y directions from a simple streamfunction
as explained at https://shaddenlab.berkeley.edu/uploads/LCS-tutorial/examples.html#Sec7.1
"""
function streamfunction_xy()
    scal = 3e-1
    f(x, y, z) = scal*(0.3*sin(x*1π/128)*sin(y*2π/128)+0.7*sin(x*2π/128)*sin(y*1π/128)) #stream function

    ϕcorners=[f(x,y,0.) for x in 0:128, y in 0:128]
    ϕcenters=[f(x,y,0.) for x in 0.5:128, y in 0.5:128]

    uu=-diff(ϕcorners,dims=2)[1:end-1,:]
    vv=diff(ϕcorners,dims=1)[:,1:end-1]
    uu=reshape(uu,(128,128,1))
    vv=reshape(vv,(128,128,1))
    ww=zeros(128,128,2)

    uvels = fill(uu, 2)
    vvels = fill(vv, 2)
    wvels = fill(ww, 2)
    uvels = cat(uvels..., dims=4)
    vvels = cat(vvels..., dims=4)
    wvels = cat(wvels..., dims=4)

    return uvels, vvels, wvels, ϕcenters
end


text:
module CompGrids

using ParallelStencil
using ImplicitGlobalGrid
using MPI
using OffsetArrays

@init_parallel_stencil(Threads, Float64, 3);

export backend,
       BackendPETSc,  BackendParallelStencil, BackendNone,

       #
       initialize_backend, initialize_grid,

       # Boundary topology

       # Grids
       RegularRectilinearCollocatedGrid, Bounded, Ghosted, Periodic

# Declare the backend we are employing:
include("Backends.jl")
include("PETSc_backend.jl")
include("ParallelStencil_backend.jl")
include("Grids.jl")

# Declare the backend we are employing:

# Define different grid types

end # module


text:
function EvalAnalDani( x, y, rc, mm, mc )
# ---------------------------------------------------------------------------
# ANALYTICAL SOLUTION - PRESSURE AND VELOCITY AROUND A CIRCULAR INCLUSION:
#
# BASED ON DANI SCHMID'S 2002 CYL_P_MATRIX.M
# FAR FIELD FLOW - VISCOSITIES - GEOMETRY
#
# ---------------------------------------------------------------------------

# INPUT:
gr  = 0;                        # Simple shear: gr=1, er=0
er  = -1;                       # Strain rate
A   =   mm*(mc-mm)/(mc+mm);
i   =   1im;

# --------------------------------------------------------------
# PRESSURE CALCULATION OUTSIDE OF AN INCLUSION IN THE y-PLANE
# --------------------------------------------------------------

# println(x)
# println(y)
# println(rc)

# INSIDE CLAST
if sqrt(x^2.0 + y^2.0)<=rc

    Z       =   x + i*y;
    P       =   0;

    # VELOCITY
    V_tot      =  (mm/(mc+mm))*(i*gr+2*er)*conj(Z)-(i/2)*gr*Z;
    vx         =  real(V_tot);
    vy         =  imag(V_tot);
    eta        =  mc;

    A   = mm*(mc-mm)/(mc+mm);
    v1x = (- 3.0*A*er*rc^4*x^4 - 6.0*A*gr*rc^4*x^3.0*y + 18.0*A*er*rc^4*x^2.0*y^2.0 + 6.0*A*gr*rc^4*x*y^3.0 - 3.0*A*er*rc^4*y^4 + 2.0*A*er*rc^2.0*x^6.0 + 4*A*gr*rc^2.0*x^5*y - 10*A*er*rc^2.0*x^4*y^2.0 - 10*A*er*rc^2.0*x^2.0*y^4 - 4*A*gr*rc^2.0*x*y^5 + 2.0*A*er*rc^2.0*y^6.0 + er*mm*x^8 + 4*er*mm*x^6.0*y^2.0 + 6.0*er*mm*x^4*y^4 + 4*er*mm*x^2.0*y^6.0 + er*mm*y^8)/(mm*(x^2.0 + y^2.0)^4);
    v1y = (3.0*A*gr*rc^4*x^4 - 24.0*A*er*rc^4*x^3.0*y - 18.0*A*gr*rc^4*x^2.0*y^2.0 + 24.0*A*er*rc^4*x*y^3.0 + 3.0*A*gr*rc^4*y^4 - 4*A*gr*rc^2.0*x^6.0 + 24.0*A*er*rc^2.0*x^5*y + 8*A*gr*rc^2.0*x^4*y^2.0 + 16*A*er*rc^2.0*x^3.0*y^3.0 + 12*A*gr*rc^2.0*x^2.0*y^4 - 8*A*er*rc^2.0*x*y^5 + 2.0*gr*mm*x^8 + 8*gr*mm*x^6.0*y^2.0 + 12*gr*mm*x^4*y^4 + 8*gr*mm*x^2.0*y^6.0 + 2.0*gr*mm*y^8)/(2.0*mm*(x^2.0 + y^2.0)^4);
    v2x = (A*rc^2.0*(3.0*gr*rc^2.0*x^4 - 24.0*er*rc^2.0*x^3.0*y - 18.0*gr*rc^2.0*x^2.0*y^2.0 + 24.0*er*rc^2.0*x*y^3.0 + 3.0*gr*rc^2.0*y^4 + 8*er*x^5*y + 12*gr*x^4*y^2.0 - 16*er*x^3.0*y^3.0 + 8*gr*x^2.0*y^4 - 24.0*er*x*y^5 - 4*gr*y^6.0))/(2.0*mm*(x^2.0 + y^2.0)^4);
    v2y = -(- 3.0*A*er*rc^4*x^4 - 6.0*A*gr*rc^4*x^3.0*y + 18.0*A*er*rc^4*x^2.0*y^2.0 + 6.0*A*gr*rc^4*x*y^3.0 - 3.0*A*er*rc^4*y^4 + 2.0*A*er*rc^2.0*x^6.0 + 4*A*gr*rc^2.0*x^5*y - 10*A*er*rc^2.0*x^4*y^2.0 - 10*A*er*rc^2.0*x^2.0*y^4 - 4*A*gr*rc^2.0*x*y^5 + 2.0*A*er*rc^2.0*y^6.0 + er*mm*x^8 + 4*er*mm*x^6.0*y^2.0 + 6.0*er*mm*x^4*y^4 + 4*er*mm*x^2.0*y^6.0 + er*mm*y^8)/(mm*(x^2.0 + y^2.0)^4);
    Txy = mm*(v1y + v2x);
    Txx = mm*(v1x + v1x);
    Tyy = mm*(v2y + v2y);
    sxx = real(-P + Txx)
    syy = real(-P + Tyy)
    sxy = real(Txy)
else
     # OUTSIDE CLAST, RESP. MATRIX
     Z              =   x + i*y;
     # PRESSURE
     P          =   -2.0*mm.*(mc-mm)./(mc+mm).*real(rc^2.0/Z.^2.0*(i*gr+2*er));

     # VELOCITY
     phi_z          = -(i/2)*mm*gr*Z-(i*gr+2*er)*A*rc^2*Z^(-1);
     d_phi_z        = -(i/2)*mm*gr + (i*gr+2*er)*A*rc^2/Z^2;
     conj_d_phi_z   = conj(d_phi_z);
     psi_z          = (i*gr-2*er)*mm*Z-(i*gr+2*er)*A*rc^4*Z^(-3);
     conj_psi_z     = conj(psi_z);

     V_tot          = (phi_z- Z*conj_d_phi_z - conj_psi_z) / (2*mm);
     vx         =  real(V_tot);
     vy         =  imag(V_tot);
     eta        = mm;

    # Evaluate stresses (Valid one)
    v1x = (2.0*er*mm)/(mm + mc);
    v1y = (gr*(3.0*mm + mc))/(2.0*(mm + mc));
    v2x = (gr*(mm - mc))/(2.0*(mm + mc));
    v2y = -(2.0*er*mm)/(mm + mc);
    Txy = mc*(v1y + v2x);
    Txx = mc*(v1x + v1x);
    Tyy = mc*(v2y + v2y);
    sxx = -P + Txx
    syy = -P + Tyy
    sxy = Txy
end

return vx, vy, P, sxx, syy, sxy, v1x, v1y, v2x, v2y, eta
end

function Tractions( x, y, rc, mm, mc, phase )

# ---------------------------------------------------------------------------
# ANALYTICAL SOLUTION - PRESSURE AND VELOCITY AROUND A CIRCULAR INCLUSION:
#
# BASED ON DANI SCHMID'S 2002 CYL_P_MATRIX.M
# FAR FIELD FLOW - VISCOSITIES - GEOMETRY
#
# ---------------------------------------------------------------------------

# INPUT:
gr  = 0;                        # Simple shear: gr=1, er=0
er  = -1;                       # Strain rate
A   =   mm*(mc-mm)/(mc+mm);
i   =   1im;

# --------------------------------------------------------------
# PRESSURE CALCULATION OUTSIDE OF AN INCLUSION IN THE y-PLANE
# --------------------------------------------------------------

# INSIDE CLAST
if phase == 2

    Z       =   x + i*y;
    P       =   0;   # if you want you can add NaN, but according to Schmid's thesis it's zero inside

    # VELOCITY
    V_tot      =  (mm/(mc+mm))*(i*gr+2*er)*conj(Z)-(i/2)*gr*Z;
    vx         =  real(V_tot);
    vy         =  imag(V_tot);
    eta        =  mc;

    # Evaluate stresses
    v1x = (2*er*mm)/(mm + mc);
    v1y = (gr*(3*mm + mc))/(2*(mm + mc));
    v2x = (gr*(mm - mc))/(2*(mm + mc));
    v2y = -(2*er*mm)/(mm + mc);

    Txy = mm*(v1y + v2x);
    Txx = mm*(v1x + v1x);
    Tyy = mm*(v2y + v2y);
    sxx = real(-P + Txx)
    syy = real(-P + Tyy)
    sxy = real(Txy)
end
if phase == 1
    # OUTSIDE CLAST, RESP. MATRIX
    Z              =   x + i*y;
    # PRESSURE
    P          =   -2.0*mm.*(mc-mm)./(mc+mm).*real(rc^2.0/Z.^2.0*(i*gr+2*er));

    # VELOCITY
    phi_z          = -(i/2)*mm*gr*Z-(i*gr+2*er)*A*rc^2*Z^(-1);
    d_phi_z        = -(i/2)*mm*gr + (i*gr+2*er)*A*rc^2/Z^2;
    conj_d_phi_z   = conj(d_phi_z);
    psi_z          = (i*gr-2*er)*mm*Z-(i*gr+2*er)*A*rc^4*Z^(-3);
    conj_psi_z     = conj(psi_z);

    V_tot          = (phi_z- Z*conj_d_phi_z - conj_psi_z) / (2*mm);
    vx         =  real(V_tot);
    vy         =  imag(V_tot);
    eta        = mm;

    # Evaluate stresses (Valid one)
    A = mm*(mc-mm)/(mc+mm);
    v1x = (- 3*A*er*rc^4*x^4 - 6*A*gr*rc^4*x^3*y + 18*A*er*rc^4*x^2*y^2 + 6*A*gr*rc^4*x*y^3 - 3*A*er*rc^4*y^4 + 2*A*er*rc^2*x^6 + 4*A*gr*rc^2*x^5*y - 10*A*er*rc^2*x^4*y^2 - 10*A*er*rc^2*x^2*y^4 - 4*A*gr*rc^2*x*y^5 + 2*A*er*rc^2*y^6 + er*mm*x^8 + 4*er*mm*x^6*y^2 + 6*er*mm*x^4*y^4 + 4*er*mm*x^2*y^6 + er*mm*y^8)/(mm*(x^2 + y^2)^4);
    v1y = (3*A*gr*rc^4*x^4 - 24*A*er*rc^4*x^3*y - 18*A*gr*rc^4*x^2*y^2 + 24*A*er*rc^4*x*y^3 + 3*A*gr*rc^4*y^4 - 4*A*gr*rc^2*x^6 + 24*A*er*rc^2*x^5*y + 8*A*gr*rc^2*x^4*y^2 + 16*A*er*rc^2*x^3*y^3 + 12*A*gr*rc^2*x^2*y^4 - 8*A*er*rc^2*x*y^5 + 2*gr*mm*x^8 + 8*gr*mm*x^6*y^2 + 12*gr*mm*x^4*y^4 + 8*gr*mm*x^2*y^6 + 2*gr*mm*y^8)/(2*mm*(x^2 + y^2)^4);
    v2x = (A*rc^2*(3*gr*rc^2*x^4 - 24*er*rc^2*x^3*y - 18*gr*rc^2*x^2*y^2 + 24*er*rc^2*x*y^3 + 3*gr*rc^2*y^4 + 8*er*x^5*y + 12*gr*x^4*y^2 - 16*er*x^3*y^3 + 8*gr*x^2*y^4 - 24*er*x*y^5 - 4*gr*y^6))/(2*mm*(x^2 + y^2)^4);
    v2y = -(- 3*A*er*rc^4*x^4 - 6*A*gr*rc^4*x^3*y + 18*A*er*rc^4*x^2*y^2 + 6*A*gr*rc^4*x*y^3 - 3*A*er*rc^4*y^4 + 2*A*er*rc^2*x^6 + 4*A*gr*rc^2*x^5*y - 10*A*er*rc^2*x^4*y^2 - 10*A*er*rc^2*x^2*y^4 - 4*A*gr*rc^2*x*y^5 + 2*A*er*rc^2*y^6 + er*mm*x^8 + 4*er*mm*x^6*y^2 + 6*er*mm*x^4*y^4 + 4*er*mm*x^2*y^6 + er*mm*y^8)/(mm*(x^2 + y^2)^4);
    Txy = mc*(v1y + v2x);
    Txx = mc*(v1x + v1x);
    Tyy = mc*(v2y + v2y);
    sxx = -P + Txx
    syy = -P + Tyy
    sxy = Txy
end

return P, real(v1x), real(v1y), real(v2x), real(v2y)
end


text:
"""
    measure(mac_or_fun::Union{Expr, Function}, tr::Trace; normalize=false,
            threshold=0, explore_worst=false)`

Apply `mac_or_fun` to every function call in `tr`, with a `green->red` color scheme.
Remove all function calls whose result is below `threshold`. `measure` runs `tr()` first
to get rid of the compile time, but that is by not always sufficient. Call `measure`
twice for best accuracy (esp. with `@allocated`).

When `explore_worst=true`, `measure` will only expand the worst child of each call.
This is a useful way of quickly exploring a trace that would be too large to profile
otherwise.
"""
function measure(mac_or_fun::Union{Expr, Function}, trace::Trace; normalize=false,
                 threshold=0, explore_worst=false)
    f = handle_mac(mac_or_fun)
    trace() # run it once, to get (at least some of) the JIT behind us
    top_value = f(trace)
    if normalize && top_value == 0
        warn("The normalization is 0 - disabling normalization")
        normalize = false
    end
    norm(x) = normalize ? (x / top_value) : x
    shorten(x) = x
    shorten(x::AbstractFloat) = signif(x, 4)
    function trav(tr, res)
        new_called = Trace[]
        called_res = map(norm∘f, tr.called)
        worst = isempty(tr.called) ? nothing : maximum(called_res)
        for (t, t_res) in zip(tr.called, called_res)
            if number(t_res) >= threshold
                push!(new_called,
                      ((!explore_worst || t_res == worst) ?
                       trav(t, t_res) :
                       Trace(t.func, t.args, t.kwargs, [], shorten(t_res))))
            end
        end
        return Trace(tr.func, tr.args, tr.kwargs, new_called, shorten(res))
    end
    tr2 = trav(trace, norm(top_value))
    # number(x) is to handle TrialEstimate
    return greenred(tr2; map=(value(tr2)==0 ? identity : x->number(x)/number(value(tr2))))
end

################################################################################
# BenchmarkTools

const judgement_colors = Dict(:invariant=>:black, :regression=>:red, :improvement=>:green)

@require BenchmarkTools begin
    using BenchmarkTools: @benchmark, Trial, TrialEstimate, prettytime, prettymemory
    using BenchmarkTools: memory, ratio, judge, Benchmark, TrialEstimate, Trial, params
    using BenchmarkTools: loadparams!, run, TrialJudgement

    export benchmark
    TraceCalls.show_val(io::IO, mime, t::Trial) =
        # The default `show` method only shows the time.
        print(io, "Trial(", prettytime(time(t)), ", ",
              prettymemory(memory(t)), ")")
    TraceCalls.show_val(io::IO, mime, t::TrialEstimate) =
        print(io, "TrialEstimate(", prettytime(time(t)), ", ",
              prettymemory(memory(t)), ")")

    """ `benchmark(tr::Trace; warmup=true)` calls `BenchmarkTools.@benchmark` on `tr`. If
    `tr` is already a benchmark, then we reuse that benchmark's parameters, so the new
    results are more comparable and faster-computed. """
    function benchmark(tr::Trace; warmup=true)
        if value(tr) isa Union{Benchmark, Trial, TrialEstimate}
            benchmarkable = apply_macro(:@benchmarkable, tr)
            loadparams!(benchmarkable, params(value(tr)))
            if warmup;
                # See https://discourse.julialang.org/t/benchmarktools-theory-and-practice/5728/2
                BenchmarkTools.warmup(benchmarkable) end
            run(benchmarkable)
        else
            apply_macro(:@benchmark, tr)
        end
    end
    number(x::TrialEstimate) = time(x)
    # That's type-piracy. Remove it once https://github.com/JuliaCI/BenchmarkTools.jl/pull/73
    # gets decided. Could also be fixed by creating our own "divide" function.
    Base.:/(a::TrialEstimate, b::TrialEstimate) = ratio(a, b)

    BenchmarkTools.ratio(group1::Group, group2::Group) =
        if_not_error(()->apply(tr->ratio(value(tr), value(group2)), group1),
                     group1, group2)
    add_color(judge::TrialJudgement) = FontColor(judgement_colors[time(judge)], judge)
    BenchmarkTools.judge(new_group::Group, old_group::Group) =
        if_not_error(()->apply(tr->add_color(judge(value(tr), value(old_group))), new_group),
                     new_group, old_group)
end

################################################################################
# TraceGroupBenchmark

struct GroupBenchmark
    estimator::Function
    groups::Vector{Group}
end
Base.getindex(gb::GroupBenchmark, ind) = gb.groups[ind]
Base.length(gb::GroupBenchmark) = length(gb.groups)
Base.endof(gb::GroupBenchmark) = length(gb)
Base.sort(gb::GroupBenchmark; kwargs...) =
    GroupBenchmark(gb.estimator, sort(gb.groups; kwargs...))

function show_gb(io::IO, mime, gb::GroupBenchmark)
    write(io, "Benchmark\n")
    for group in gb.groups
        write(io, " ")
        show(io, mime, group)
        write(io, "\n")
    end
end
Base.show(io::IO, mime::MIME"text/html", gb::GroupBenchmark) = show_gb(io, mime, gb)
Base.show(io::IO, mime::MIME"text/plain", gb::GroupBenchmark) = show_gb(io, mime, gb)
Base.time(grp::Group) = value(grp) isa Exception ? 10000 : length(grp) * time(value(grp))
signature(grp::Group) = signature(grp[1])

@require BenchmarkTools begin
    using BenchmarkTools: run

    trace_benchmark(trace::Trace; estimator=median) =
        sort(run(GroupBenchmark(estimator, group_by(signature, map(tr->:not_run, trace)))),
             by=time, rev=true)
    trace_benchmark(code::Function, to_trace; kwargs...) =
        trace_benchmark(trace(code, to_trace); kwargs...)
    trace_benchmark(file_to_include::String, to_trace; kwargs...) =
        define_benchmark(()->include(file_to_include), to_trace; kwargs...)
    BenchmarkTools.run(gb::GroupBenchmark) =
        GroupBenchmark(gb.estimator, map_groups(gb.estimator∘benchmark, gb.groups))
    BenchmarkTools.ratio(gb1::GroupBenchmark, gb2::GroupBenchmark) =
        GroupBenchmark(gb1.estimator, map(ratio, gb1.groups, gb2.groups))
    function BenchmarkTools.judge(new_gb::GroupBenchmark, old_gb::GroupBenchmark;
                                  sort=true)
        @assert(map(signature, new_gb.groups) == map(signature, old_gb.groups),
                "The two GroupBenchmarks are not in the same order. The new benchmark should be computed by calling `run(old_benchmark)`")
        gb = GroupBenchmark(new_gb.estimator, map(judge, new_gb.groups, old_gb.groups))
        by(grp) = (value(grp) isa Exception ? -1000000 : -(time(ratio(value(grp)))-1))
        return sort ? Base.sort(gb; by=by) : gb
    end
end

################################################################################

"""
    compilation_times(to_trace, trace::Trace; warm_up=true,
                      comp_fun=Base.precompile)

Compute the time to compile each distinct call signature in `trace`, then aggregate
it by method to return a report of how much total time is spent compiling each function
called by `trace`.
"""
function compilation_times(to_trace, trace::Trace; warm_up=true,
                           comp_fun=Base.precompile)
    # My understanding is that `precompile(f)` triggers compilation (or at least,
    # inference) on some of the functions that `f` calls. Thus, we iterate
    # `collect(trace)[end:-1:1]` so that the leaves are precompiled first.
    # It would potentially be more exact if I used:
    # https://github.com/JuliaLang/julia/blob/903644385b91ed8d95e5e3a5716c089dd1f1b08a/base/reflection.jl#L710-L714
    signatures = OrderedSet(signature(tr) for tr in collect(trace)[end:-1:1])
    do_comp(signature) = comp_fun(signature[1], signature[2:end])
    # Recommended, since otherwise compilation times are reported much higher.
    # My guess is that `do_comp(f, type_tuple)` spends a lot of time generating
    # code for `type_tuple`. Not 100% sure, though.
    if warm_up; foreach(do_comp, signatures) end
    @trace to_trace nothing   # reset Julia's method table for these functions
    dict = Dict{Any, Float64}()
    for signature in signatures
        method = which(signature[1], signature[2:end])
        dict[method] = get(dict, method, 0) + @elapsed(do_comp(signature))
    end
    return OrderedDict(sort(collect(dict), by=last, rev=true)...)
end

""" `@compilation_times to_trace expr` (similar to `@trace to_trace expr`) executes
`expr` while tracing the modules/functions/files in `to_trace`, then returns an estimate
of how much time it takes to compile each of the method specializations called within
`expr`. """
macro compilation_times(to_trace, expr)
    @gensym ev_to_trace trace
    esc(quote
        $ev_to_trace = $to_trace
        $trace = $TraceCalls.@trace($ev_to_trace, $expr)
        $TraceCalls.compilation_times($ev_to_trace, $trace)
        end)
end


text:
"""Julia program to implement Minimum XOR value.
Given an array of integers, we have to find the minimum XOR value, a pair in the array makes."""

function minimumXOR(arr, n)
    # Sort the array, so that elements with minimum bit difference are present consecutive.
    sort(arr)
    minXOR = typemax(Int64)
    val = 0

    # Compute the XOR values of the consecutive elements of the array.
    for i in 1:(n-1)
        val = arr[i] ⊻ arr[i+1]
        minXOR = min(minXOR, val)
    end

    return minXOR
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("No numbers are present in the array!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
res = minimumXOR(arr, n)
print("The minimum XOR value present in the given array is $res")

"""
SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers are present in the array? 8
Enter the numbers: 25 45 36 458 651 89 52 95
The minimum XOR value present in the given array is 9

SAMPLE II

How many numbers are present in the array? 5
Enter the numbers: 5 4 3 2 1
The minimum XOR value present in the given array is 1

"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
"""
Defines universal physical constants and unit conversions
Significant digits/uncertainties NOT defined
Certainly not more than 3-4 sig figs since G is not known better than 4 sig figs

Sources: Google, Wikipedia, NIST
Jeff Kaplan <jeffkaplan@caltech.edu>
February 2, 2013
"""

#################################
# Physical constants
#
# G - Gravitational constant
# C - Speed of light in vacuum
# AMU - Atomic mass unit
# H - Plank's constant
# HBAR - Plank's constant/2pi
# K  -  Boltzmann constant

#################################
# Astronomical constants
#
# MSUN - Mass of the sun
# RSUN - Radius of the sun
# LSUN - Luminosity of the sun
# TSUN - Surface temperature of the sun
# M_EARTH - Mass of the earth
# R_EARTH - Radius of the earth
# AU  - Astronomical unit
# PC  - Parsec
# LYR - Light year

##===================================================================
## Meters kilograms seconds (SI)
##===================================================================
MKS_G = 6.67384e-11 # m^3 kg^-1 s^-2
MKS_C = 299792458.0 # m s^-1
MKS_AMU = 1.660538921e-27 # kg
MKS_H = 6.62606957e-34     # J s = m^2 kg / s
MKS_HBAR = 1.054571726e-34 # J s = m^2 kg / s
MKS_K = 1.3806488e-23  # J K^-1

MKS_MSUN = 1.98855e30  # kg
MKS_RSUN = 6.955e8     # m
MKS_LSUN = 3.933e26    # W -Watt
MKS_TSUN = 5.780e3     # K
MKS_M_EARTH = 5.972e24 # kg
MKS_R_EARTH = 6478100.0 # m
MKS_AU   = 1.49597871e11 # m
MKS_PC   = 3.08567758e16  # m
MKS_LYR  = 9.4605284e15  # M

#################################
# Conversion constants
#
# EV - Electron volt
# CM - Centimeter: cgs length unit
# GRAM - Gram:   cgs mass unit
# ERG - Erg:   cgs energy unit
# DYNE - Dyne: cgs force unit
# GRAM_PER_CM_CUBED: density conversion to cgs
# DYNE_PER_CM_SQUARED: pressure conversion to cgs
MKS_EV = 1.60217646-19 # J
MKS_CM = 0.01          # m
MKS_GRAM = 0.001       # kg
MKS_ERG  = 1.0e-7      # J
MKS_DYNE = 1.0e-5      # N = m  kg s^-2
MKS_GRAM_PER_CM_CUBED = 1000.0 # kg m ^-3
MKS_DYNE_PER_CM_SQUARED = 0.1 # Pa = N m^-2  = kg m^-1 s^-2

##===================================================================
## Centimeters grams seconds
##===================================================================
CGS_G = 6.67384e-8 # cm^3 g^-1 s^-2
CGS_C = 29979245800. # cm s^-1
CGS_AMU = 1.660538921e-24 # g
CGS_H = 6.62606957e-27     # erg s = cm^2 g / s
CGS_HBAR = 1.054571726e-27 # erg s = cm^2 g / s
CGS_K = 1.3806488e-16  # erg K^-1

CGS_MSUN = 1.98855e33  # g
CGS_RSUN = 6.955e10     # cm
CGS_LSUN = 3.933e33    # erg s^-1
CGS_TSUN = 5.780e3     # K
CGS_M_EARTH = 5.972e27 # g
CGS_R_EARTH = 647810000. # cm
CGS_AU   = 1.49597871e13 # cm
CGS_PC   = 3.08567758e18  # cm
CGS_LYR  = 9.4605284e17  # cm
CGS_sigmaT = 6.6524e-25 # cm^2
mp = 1.6726231e-24 # g

#################################
# Conversion constants
#
# EV - Electron volt
# M - Meter: mks length unit
# KG - Kilogram:   mks mass unit
# J - Joule:   mks energy unit
# N - Newton: mks force unit
# KG_PER_M_CUBED: density conversion to mks
# N_PER_M_SQUARED: pressure conversion to mks
CGS_EV = 1.60217646-12 # erg
CGS_M = 100.0          # cn
CGS_KG = 1000.0       # g
CGS_J  = 1.0e+7      # erg
CGS_N = 1.0e+5      # dyne = cm  g s^-2
CGS_KG_PER_M_CUBED = 0.001 # g cm ^-3
CGS_N_PER_M_SQUARED = 10.0 # ba = dyne cm^-2  = g cm^-1 s^-2

##===================================================================
## Augmented geometric units G=c=Msun=1
#  Units are dimensionless, except Boltzmann K which has Kelvin^-1
#  S - second     W - Watt
##===================================================================
AGEO_G = 1.0
AGEO_C = 1.0
AGEO_AMU = 8.3505012e-58
AGEO_K  = 7.7229851e-71   # K^-1

AGEO_MSUN = 1.0
AGEO_RSUN = 470934.0
AGEO_LSUN = 4.4659113e-16

#################################
# Conversion constants
AGEO_EV = 8.9621524e-67
AGEO_M = 6.77116916e-4
AGEO_CM = 6.77116916e-6
AGEO_S  = 202994.545
AGEO_KG = 5.02739933e-31
AGEO_GRAM =5.02739933e-34
AGEO_KG_PER_M_CUBED = 1.6193935e-21
AGEO_GRAM_PER_CM_CUBED = 1.6193935e-18
AGEO_N_PER_M_SQUARED = 1.80181827e-38
AGEO_DYNE_PER_CM_SQUARED = 1.80181827e-37
AGEO_J = 5.59373614e-48
AGEO_ERG = 5.59373614e-55
AGEO_W  = 1.13549741e-42
AGEO_ERGS_PER_S = 1.13549741e-49

#################################
# Inverse conversion constants
# Gives you a feel for what '1' is for quantities in AGEO units
AGEO_LENGTH_IN_M = 1.47684983e3  # m
AGEO_LENGTH_IN_CM = 1.47684983e5  # cm Approx 1.5 km
AGEO_TIME_IN_S = 4.92624076e-3   # s   Approx 5 ms
AGEO_DENSITY_IN_MKS = 6.17515138e20  # kg m^-3
AGEO_DENSITY_IN_CGS = 6.17515138e17  # g cm^-3  Approx 6000x nuclear density
AGEO_PRESSURE_IN_MKS = 5.54994928e37 # Pa
AGEO_PRESSURE_IN_CGS = 5.54994928e36 # ba = dyne cm^-2
AGEO_ENERGY_IN_MKS = 1.78771393e47 # Joules
AGEO_ENERGY_IN_CGS = 1.78771393e54 # ergs
AGEO_POWER_IN_MKS  = 8.80671314e41 # W
AGEO_POWER_IN_CGS  = 8.80671314e48 # erg s^-1


text:
#
# Interaction of GrooupManifold with a ValidationaManifold
#
array_value(e::Identity) = e

array_point(p) = ValidationMPoint(p)
array_point(p::ValidationMPoint) = p

const ValidationGroup{𝔽} = ValidationManifold{𝔽,G} where {G<:AbstractGroupManifold}

function adjoint_action(M::ValidationGroup, p, X; kwargs...)
    is_point(M, p, true; kwargs...)
    eM = Identity(M.manifold)
    is_vector(M, eM, X, true; kwargs...)
    Y = ValidationTVector(adjoint_action(M.manifold, array_value(p), array_value(X)))
    is_vector(M, eM, Y, true; kwargs...)
    return Y
end

function adjoint_action!(M::ValidationGroup, Y, p, X; kwargs...)
    is_point(M, p, true; kwargs...)
    eM = Identity(M.manifold)
    is_vector(M, eM, X, true; kwargs...)
    adjoint_action!(M.manifold, array_value(Y), array_value(p), array_value(X))
    is_vector(M, eM, Y, true; kwargs...)
    return Y
end

Identity(M::ValidationGroup) = array_point(Identity(M.manifold))
identity_element!(M::ValidationGroup, p) = identity_element!(M.manifold, array_value(p))

function Base.inv(M::ValidationGroup, p; kwargs...)
    is_point(M, p, true; kwargs...)
    q = array_point(inv(M.manifold, array_value(p)))
    is_point(M, q, true; kwargs...)
    return q
end

function inv!(M::ValidationGroup, q, p; kwargs...)
    is_point(M, p, true; kwargs...)
    inv!(M.manifold, array_value(q), array_value(p))
    is_point(M, q, true; kwargs...)
    return q
end

function lie_bracket(M::ValidationGroup, X, Y)
    eM = Identity(M.manifold)
    is_vector(M, eM, X, true)
    is_vector(M, eM, Y, true)
    Z = ValidationTVector(lie_bracket(M.manifold, array_value(X), array_value(Y)))
    is_vector(M, eM, Z, true)
    return Z
end

function lie_bracket!(M::ValidationGroup, Z, X, Y)
    eM = Identity(M.manifold)
    is_vector(M, eM, X, true)
    is_vector(M, eM, Y, true)
    lie_bracket!(M.manifold, array_value(Z), array_value(X), array_value(Y))
    is_vector(M, eM, Z, true)
    return Z
end

function compose(M::ValidationGroup, p, q; kwargs...)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    x = array_point(compose(M.manifold, array_value(p), array_value(q)))
    is_point(M, x, true; kwargs...)
    return x
end

function compose!(M::ValidationGroup, x, p, q; kwargs...)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    compose!(M.manifold, array_value(x), array_value(p), array_value(q))
    is_point(M, x, true; kwargs...)
    return x
end

function translate(M::ValidationManifold, p, q, conv::ActionDirection; kwargs...)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    x = array_point(translate(M.manifold, array_value(p), array_value(q), conv))
    is_point(M, x, true; kwargs...)
    return x
end

function translate!(M::ValidationManifold, x, p, q, conv::ActionDirection; kwargs...)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    translate!(M.manifold, array_value(x), array_value(p), array_value(q), conv)
    is_point(M, x, true; kwargs...)
    return x
end

function inverse_translate(M::ValidationManifold, p, q, conv::ActionDirection; kwargs...)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    x = array_point(inverse_translate(M.manifold, array_value(p), array_value(q), conv))
    is_point(M, x, true; kwargs...)
    return x
end

function inverse_translate!(
    M::ValidationManifold,
    x,
    p,
    q,
    conv::ActionDirection;
    kwargs...,
)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    inverse_translate!(M.manifold, array_value(x), array_value(p), array_value(q), conv)
    is_point(M, x, true; kwargs...)
    return x
end

function translate_diff(M::ValidationManifold, p, q, X, conv::ActionDirection; kwargs...)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    is_vector(M, q, X, true; kwargs...)
    Y = ValidationTVector(
        translate_diff(M.manifold, array_value(p), array_value(q), array_value(X), conv),
    )
    pq = translate(M, p, q, conv)
    is_vector(M, pq, Y, true; kwargs...)
    return Y
end

function translate_diff!(
    M::ValidationManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection;
    kwargs...,
)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    is_vector(M, q, X, true; kwargs...)
    translate_diff!(
        M.manifold,
        array_value(Y),
        array_value(p),
        array_value(q),
        array_value(X),
        conv,
    )
    pq = translate(M, p, q, conv)
    is_vector(M, pq, Y, true; kwargs...)
    return Y
end

function inverse_translate_diff(
    M::ValidationManifold,
    p,
    q,
    X,
    conv::ActionDirection;
    kwargs...,
)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    is_vector(M, q, X, true; kwargs...)
    Y = ValidationTVector(
        inverse_translate_diff(
            M.manifold,
            array_value(p),
            array_value(q),
            array_value(X),
            conv,
        ),
    )
    pinvq = inverse_translate(M, p, q, conv)
    is_vector(M, pinvq, Y, true; kwargs...)
    return Y
end

function inverse_translate_diff!(
    M::ValidationManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection;
    kwargs...,
)
    is_point(M, p, true; kwargs...)
    is_point(M, q, true; kwargs...)
    is_vector(M, q, X, true; kwargs...)
    inverse_translate_diff!(
        M.manifold,
        array_value(Y),
        array_value(p),
        array_value(q),
        array_value(X),
        conv,
    )
    pinvq = inverse_translate(M, p, q, conv)
    is_vector(M, pinvq, Y, true; kwargs...)
    return Y
end

function exp_lie(M::ValidationManifold, X; kwargs...)
    is_vector(
        M,
        Identity(M.manifold),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    q = array_point(exp_lie(M.manifold, array_value(X)))
    is_point(M, q, true; kwargs...)
    return q
end

function exp_lie!(M::ValidationManifold, q, X; kwargs...)
    is_vector(
        M,
        Identity(M.manifold),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    exp_lie!(M.manifold, array_value(q), array_value(X))
    is_point(M, q, true; kwargs...)
    return q
end

function log_lie(M::ValidationManifold, q; kwargs...)
    is_point(M, q, true; kwargs...)
    X = ValidationTVector(log_lie(M.manifold, array_value(q)))
    is_vector(
        M,
        Identity(M.manifold),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    return X
end

function log_lie!(M::ValidationManifold, X, q; kwargs...)
    is_point(M, q, true; kwargs...)
    log_lie!(M.manifold, array_value(X), array_value(q))
    is_vector(
        M,
        Identity(M.manifold),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    return X
end


text:
# ogm_ls.jl

using MIRT: ogm_ls

using LinearAlgebra: norm, opnorm, I
using Random: seed!
using Plots: plot, plot!, scatter!, gui
using LaTeXStrings
using Test: @test, @test_throws, @inferred

seed!(0); M = 30; N = 6; A = randn(M,N); y = randn(M)
a2 = opnorm(A)^2
reg = 0.1 * a2
xh = (A'A + reg*I) \ A'y
cost = (x) -> 1/2 * norm(A * x - y)^2 + reg/2 * norm(x)^2
fun = (x,iter) -> (cost(x) - cost(xh), norm(x - xh) / norm(xh), time())
grad1 = (x) -> A' * (A * x - y) + reg * x
curv1 = (x) -> a2 + reg

B = [A, I] # matrix blocks
gradf = [u -> u - y, v -> reg * v] # f functions gradients
curvf = [v -> 1, v -> reg]

niter = 90
#x1, out1 = @inferred ogm_ls(   grad1, curv1, zeros(N), niter=niter, fun=fun) todo: fails
x1, out1 = ogm_ls(   grad1, curv1, zeros(N), niter=niter, fun=fun)
x2, out2 = ogm_ls(B, gradf, curvf, zeros(N), niter=niter, fun=fun)

#@test x1 ≈ xh # no, converges too slowly
@test x2 ≈ xh

if isinteractive()
	lf = x -> log10(max(x,1e-17))
	costk = out -> lf.([out[k][1] for k=1:niter+1])
	errk = out -> lf.([out[k][2] for k=1:niter+1])
	allk = out -> (costk(out), errk(out))
	cost1, err1 = allk(out1)
	cost2, err2 = allk(out2)

	k = 0:niter
	plot(xlabel="k", ylabel=L"\log(\Psi(x_k) - \Psi(x_*))")
	scatter!(k, cost1, color=:blue, label="cost1")
	scatter!(k, cost2, color=:red, marker=:x, label="cost2")
	p1 = plot!()

	plot(xlabel="k", ylabel=L"\log(\|x_k - x_*\|/\|x_*\|)")
	scatter!(k, err1, color=:blue, label="NRMSD1")
	scatter!(k, err2, color=:red, marker=:x, label="NRMSD2")
	p2 = plot!()
	plot(p1, p2)
	gui()
end


text:
# An additive zero mean Gaussian noise Unscented Kalman Filter example,
#shamelessly taken from a paper by Rambabu Kandepu, Bjarne Foss and Lars Imsland
#"Applying the unscented Kalman filter for nonlinear state estimation"
#Link to the pdf version of the paper is here:
#http://folk.ntnu.no/bjarnean/pubs/journals/journal-47.pdf
#The example is the Van der Pol oscillator (Section 4.1 on page 6 of the paper)
#so see this paper for the full details of the problem. Here we only state
#enough to set up the Kalman Filter problem and solve it.

#Here we import the required modules
using StateSpace
using Distributions
using Gadfly
using DataFrames
using Colors

################################################################################
#Section: Set the Unscented Kalman Filter Parameters
#-------------------------------------------------------------------------------
#Here we're going to set the process model and the observation model functions
#for the oscillator along with the corresponding covarience matrices.
Δt = 0.1 # Set the time step

#Set the process function
function processFunction(x::Vector, Δt::Float64=0.1, μ::Float64=0.3)
    x1 = zeros(x)
    x1[1] = x[1] + Δt * -x[2]
    x1[2] = x[2] + Δt * (-μ * (1 - x[1]^2) * x[2] + x[1])
    x1
end

#Set the observation function.
observationFunction(x::Vector) = x

#Set process noise covariance matrix
processCovariance = 1e-2*[0.1 0;
    0 1e-3]

#Set the observation noise covariance
observationCovariance = 1e-1*eye(2)

#Create additive noise UKF model
ukfStateModel = NonlinearGaussianSSM(processFunction, processCovariance,
observationFunction, observationCovariance)
#End Section: Set the Unscented Kalman Filter Parameters
################################################################################

################################################################################
#Section: Generate noisy observations
#-------------------------------------------------------------------------------
#Set the true initial state of the oscillator
trueInitialState = [1.4, 0.0]

#Set the number of observations and generate the true values and the noisy
#observations.
numObs = 200
trueState = zeros(2,numObs)
noisyObs = zeros(2,numObs)
trueState[:,1] = trueInitialState
noisyObs[:,1] = trueInitialState + sqrt(observationCovariance)*randn(2)
for i in 2:numObs
    trueState[:,i] = processFunction(trueState[:,i-1])
    noisyObs[:,i] = trueState[:,i] + sqrt(observationCovariance)*randn(2)
end
#End Section: Generate noisy observations
################################################################################

################################################################################
#Section: Set guess of initial state
#-------------------------------------------------------------------------------
initial_guess = MvNormal([0.0,5.0], 5.0*eye(2))
#End Section: Set guess of initial state
################################################################################

################################################################################
#Section: Execute the Additive noise Unscented Kalman Filter
#-------------------------------------------------------------------------------
filtered_state = filter(ukfStateModel, noisyObs, initial_guess, UKF())
#End Section: Execute the Additive noise Unscented Kalman Filter
################################################################################

################################################################################
#Section: Plot results
#-------------------------------------------------------------------------------
x_data = 1:numObs
x1_array = Vector{Float64}(numObs)
x1Var_array = Vector{Float64}(numObs)
x1_Guess = initial_guess.μ[1]
x1Var_Guess = 2*sqrt(initial_guess.Σ.mat[1,1])
for i in x_data
    current_state = filtered_state.state[i]
    x1_array[i] = current_state.μ[1]
    x1Var_array[i] = 2*sqrt(current_state.Σ.mat[1,1])
end
df_fs = DataFrame(
    x = x_data*Δt,
    y = x1_array,
    ymin = x1_array - x1Var_array,
    ymax = x1_array + x1Var_array,
    f = "Filtered values"
    )

n = 3
getColors = distinguishable_colors(n, Color[LCHab(70, 60, 240)],
                                   transform=c -> deuteranopic(c, 0.5),
                                   lchoices=Float64[65, 70, 75, 80],
                                   cchoices=Float64[0, 50, 60, 70],
                                   hchoices=linspace(0, 330, 24))
oscillatorx1_state_plot = plot(
    layer(x=x_data*Δt, y=noisyObs[1,:], Geom.point, Theme(default_color=getColors[2])),
    layer(x=x_data*Δt, y=trueState[1,:], Geom.line, Theme(default_color=getColors[3])),
    layer(df_fs, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Time (seconds)"), Guide.ylabel("x1"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Unscented Kalman Filter (Additive) Example")
    )
display(oscillatorx1_state_plot)
#End Section: Plot results
################################################################################