text:
"""
Fusion of Geometry and Color Information for Scene Segmentation
"""
module GCF
export Config

using Statistics
using RGBDSegmentation
using Images
using ImageSegmentationEvaluation
using ImageSegmentation
import RGBDSegmentation:
    clusterize,
    CDNImage,
    rgb2lab!,
    colors,
    distances,
    number_of_pixels
struct Config
    nev::Integer
    lambdas
    proportion_of_pixels::Float64
    function Config(; nev::Integer=5,
                    lambdas::Array=[0.3,0.4],
                    proportion_of_pixels::Float64=0.001)
        if proportion_of_pixels <0 || proportion_of_pixels >1
            throw("Invalid Proportion")
        end
        return new(nev, lambdas, proportion_of_pixels)
    end
end
using Clustering
using SpectralClustering
using Distances

function weight(i::Integer,
                j::Vector{<:Integer},
                pixel_i,
                neighbors_data)
    dist = Distances.colwise(SqEuclidean(),
                             pixel_i[3:8],
                             neighbors_data[3:8, :])

    return exp.(-(dist))
end
function process_labels(labels::SegmentedImage, min_size=350)
    deletion_rule = i -> segment_pixel_count(labels, i) < min_size
    replacement_rule = (i, j) -> -segment_pixel_count(labels, j)
    return prune_segments(labels, deletion_rule, replacement_rule)
end
function clusterize(cfg::Config, img::CDNImage, k::Integer)
    n_of_pixels = number_of_pixels(img)
    best_clus = []
    cluster_resu = []
    best_perf = -Inf
    curr_perf = -5000
    for lambda in cfg.lambdas
        @info "$lambda $curr_perf $best_perf"
        n_of_samples = round(Int,
                             n_of_pixels * cfg.proportion_of_pixels)
        img_a = copy(img)
        rgb2lab!(img_a)
        std_per_dimension = dropdims(std(img_a, dims=[2,3]), dims=(2,3))
        img_a[1:3, :, :] .*= 3/(sum(std_per_dimension[1:3]))
        img_a[4:6, :, :] .*= lambda*(3/(sum(std_per_dimension[4:6])))

        nystrom = NystromMethod(
            EvenlySpacedLandmarkSelection(),
            n_of_samples,
            weight,
            cfg.nev,
            true)
        cluster_resu = SpectralClustering.clusterize(
            nystrom,
            KMeansClusterizer(k),
            img_a)

        labels = Images.label_components(
            reshape(
                assignments(cluster_resu),
                (size(img_a, 2), size(img_a, 3))))

        segmented_image = SegmentedImage(img, labels)
        labels = process_labels(segmented_image)
        curr_perf = ImageSegmentationEvaluation.evaluate(
            FRCRGBD(),
            colors(img_a),
            z(img_a),
            labels)
        if (curr_perf > best_perf)
            best_perf = curr_perf
            best_clus = segmented_image
        end
    end
    return best_clus
end
end


text:
# Plot text labels at arbitrary coordinate points
function PlotText(text::String, x::Real, y::Real; vertical::Bool = false,
                  pixel_offset::ImVec2 = ImVec2(0,0))
    LibCImPlot.PlotText(text, x, y, vertical, pixel_offset)
end

function PlotImage(user_texture_id, bounds_min::ImVec2 = ImVec2(0,0),
                   bounds_max::ImVec2 = ImVec2(1,1),
                   uv0::ImVec2 = (0,0), uv1::ImVec2=(1,1), tint_col::ImVec4=(1,1,1,1))
    LibCImPlot.PlotImage(user_texture_id, bounds_min, bounds_max, uv0, uv1, tint_col)
end

ext: jl


text:
using .Scattering
import Base: length

abstract type Method end
phaseshift(k₀::Real, m, P::Potential, method::Method; kwargs...) = method(k₀, m, P; kwargs...)
phaseshift(k₀::AbstractArray, m, P::Potential, method::Method; kwargs...) = [method(k, m, P; kwargs...) for k in k₀]
length(::Method) = 1
crossection(δₗ::Real, k::Real) = 4π*sin(δₗ)^2/k^2
crossection(δₗ::Vector{<:Real}, k::Vector{<:Real}) = @. 4π*sin(δₗ)^2/k^2

max_stars_repo_name: ErlendLima/Scattering

max_issues_repo_name: ErlendLima/Scattering

max_forks_repo_name: ErlendLima/Scattering


text:
module TMP

using Gridap
using ChainRulesCore
import ChainRulesCore: rrule

domain = (0,1,0,1)
cells = (10,10)
model = CartesianDiscreteModel(domain,cells)

order = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

# Manufactured solution
u(x) = x[1]+x[2]

# Weak form functions
a(p,u,v) = ∫( (p*p)*∇(v)⋅∇(u) )dΩ
l(v) = ∫( -v*Δ(u) )dΩ

# Space for the design dofs
reffe_S = ReferenceFE(lagrangian,Float64,0)
S = TestFESpace(model,reffe_S)
R = TrialFESpace(S)

# Space for the material params
reffe_Q = ReferenceFE(lagrangian,Float64,1)
Q = TestFESpace(model,reffe_Q)
P = TrialFESpace(Q)

# Space for the primal solution
reffe_V = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe_V,dirichlet_tags="boundary")
U = TrialFESpace(V,u)

"""
    f(uvec) -> Number

Objective function
"""
function f(uvec)
  function α(uh)
    eh = u-uh
    eh*eh
  end
  uh = FEFunction(U,uvec)
  sum(∫(α(uh))dΩ)
end

"""
    uvec = g(pvec)

Direct problem from material parameters
"""
function g(pvec)
  ph = FEFunction(P,pvec)
  op = AffineFEOperator((u,v)->a(ph,u,v),l,U,V)
  uh = solve(op)
  uvec = get_free_values(uh)
  uvec
end

"""
    pvec = h(rvec)

Material parameters from design dofs `rvec`.
"""
function h(rvec)
  rh = FEFunction(R,rvec)
  op = AffineFEOperator((p,q)->h1(p,q),q->h1(rh,q),P,Q)
  ph = solve(op)
  pvec = get_free_values(ph)
  pvec
end
h1(p,q) = ∫( ∇(p)⋅∇(q) + p*q )dΩ

"""
   F(rvec) -> Number

Objective function from design dofs
"""
function F(rvec)
  pvec = h(rvec) # Smoothing
  uvec = g(pvec) # Forward PDE
  f(uvec) # eval objective
end

# reverse rules

function rrule(::typeof(f),uvec)
  function f_pullback(dFdF)
    NO_FIELDS, dFdF*df_du(uvec)
  end
  f(uvec), f_pullback
end

function df_du(uvec)
  function dα(uh,du)
    eh = uh-u
    de = du
    ∫( 2*eh*de )dΩ
  end
  uh = FEFunction(U,uvec)
  assemble_vector(du->dα(uh,du),V)
end

function rrule(::typeof(g),pvec)
  uvec = g(pvec)
  function g_pullback(dFdu)
    NO_FIELDS, dg_dp(dFdu,uvec,pvec)
  end
  uvec, g_pullback
end

function dg_dp(dFdu,uvec,pvec)

  # Adjoint problem
  ph = FEFunction(P,pvec)
  A = assemble_matrix((λ,v)->a(ph,v,λ),U,V)
  λvec = A\dFdu

  # Total derivative wrt pvec
  da(p,u,v,dp) = ∫( (2*p*dp)*∇(v)⋅∇(u) )dΩ
  uh = FEFunction(U,uvec)
  λh = FEFunction(V,λvec)
  -1*assemble_vector((dp)->da(ph,uh,λh,dp),Q)
end

function rrule(::typeof(h),rvec)
  pvec = h(rvec)
  function h_pullback(dFdp)
    NO_FIELDS, dh_dr(dFdp)
  end
  pvec, h_pullback
end

function dh_dr(dFdp)
  A = assemble_matrix((p,q)->h1(q,p),P,Q)
  ϕvec = A\dFdp

  ϕh = FEFunction(S,ϕvec)
  assemble_vector((dr)->h1(dr,ϕh),S)
end

rvec = rand(num_free_dofs(R))
@show F(rvec)

pvec, h_pullback = rrule(h,rvec)
uvec, g_pullback = rrule(g,pvec)
obj, f_pullback = rrule(f,uvec)

dFdF = 1
_, dFdu = f_pullback(dFdF)
_, dFdp = g_pullback(dFdu)
_, dFdr = h_pullback(dFdp)

display(dFdr)

@show num_free_dofs(U)
@show num_free_dofs(P)
@show num_free_dofs(R)

using Zygote
dFdr2, = Zygote.gradient(F,rvec)

display(dFdr2)

end # module

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
NumericalTricks.volume(eM::eMesh{T1,Nothing}) where {T1} = error("meshes without tets have no volume")
function NumericalTricks.volume(eM::eMesh{T1,Tet}) where {T1}
    vol = 0.0
    for k = 1:n_tet(eM)
        vol += volume(eM.point[eM.tet[k]])
    end
    return vol
end

NumericalTricks.area(eM::eMesh{Nothing,T2}) where {T2} = error("meshes without triangles have no area")
function NumericalTricks.area(eM::eMesh{Tri,T2}) where {T2}
    area_cum = 0.0
    for k = 1:n_tri(eM)
        area_cum += area(vertex_pos_for_tri_ind(eM, k))
    end
    return area_cum
end

NumericalTricks.area(a::BB)   where {BB <: BoundingBox} = 8 * dot(a.e, SVector{3,Float64}(a.e[2], a.e[3], a.e[1]))
NumericalTricks.volume(a::BB) where {BB <: BoundingBox} = 8 * prod(a.e)


text:
function calc_h_para(reactor, R_0, B_0)
  Q_term = ( reactor.Q + 5 ) / ( reactor.Q * 5 )

  Q_term * reactor.P_F * B_0 / R_0
end


text:
module TestExamples2DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "MHD" begin
  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.000112085405989522, 5.908049046758828e-6, 5.908049046772597e-6, 8.493815951190774e-6, 1.297306491315482e-6, 1.2316941466595997e-6, 1.2316941466702522e-6, 1.8368487700950752e-6, 4.4903164608485606e-7],
      linf = [0.0002693844074639351, 1.6289386154083596e-5, 1.6289386155082797e-5, 2.747802438275715e-5, 5.344948833307939e-6, 8.137608476288527e-6, 8.137608476843639e-6, 1.2121292898431557e-5, 2.306551665486275e-6])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.726186766315203e-6, 8.427427215633836e-7, 8.427427215591145e-7, 1.2159504573391751e-6, 1.0090006943334665e-6, 1.0388026398212963e-6, 1.0388026398262758e-6, 1.2752164698175327e-6, 4.4590538253070487e-7],
      linf = [9.843039574297663e-6, 7.349930834232854e-6, 7.349930834218976e-6, 1.0468740074986993e-5, 4.9775926587170005e-6, 7.752322894960528e-6, 7.752322895182573e-6, 1.0724791524654997e-5, 2.275855084896829e-6],
      volume_flux = (flux_derigs_etal, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"),
      l2   = [3.7762324533854616e-6, 1.5534623833573546e-6, 1.4577234868196855e-6, 1.7647724628707057e-6, 1.4831911814574333e-6, 1.456369119716533e-6, 1.4115666913995062e-6, 1.804758237422838e-6, 8.320469738087189e-7],
      linf = [3.670661330201774e-5, 1.530289442645827e-5, 1.3592183785327006e-5, 1.5173897443654383e-5, 9.43771379136038e-6, 1.0906323046233624e-5, 1.0603954940346938e-5, 1.5900499596113726e-5, 5.978772247650426e-6],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.03637302248881514, 0.043002991956758996, 0.042987505670836056, 0.02574718055258975, 0.1621856170457943, 0.01745369341302589, 0.017454552320664566, 0.026873190440613117, 5.336243933079389e-16],
      linf = [0.23623816236321427, 0.3137152204179957, 0.30378397831730597, 0.21500228807094865, 0.9042495730546518, 0.09398098096581875, 0.09470282020962917, 0.15277253978297378, 4.307694418935709e-15])
  end

  @trixi_testset "elixir_mhd_orszag_tang.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
      l2   = [0.21671818699492046, 0.26364533451882755, 0.31388339884931576, 0.0, 0.5122729051056656, 0.2291828683723837, 0.34308184804821307, 0.0, 0.003210081660131638],
      linf = [1.2648775120640154, 0.673836824735164, 0.8585376516993447, 0.0, 2.8148181376188246, 0.657323372510843, 0.9595541188847632, 0.0, 0.0525390951527749],
      tspan = (0.0, 0.09))
  end

  @trixi_testset "elixir_mhd_orszag_tang.jl with flux_hll" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
      l2   = [0.10793703703273708, 0.20177434827232335, 0.2296817363329907, 0.0, 0.2994119440979309, 0.15675567943351448, 0.24281245338083307, 0.0, 0.003548624701792127],
      linf = [0.560160833798496, 0.5098933680011135, 0.6566913038761877, 0.0, 0.9905839416293134, 0.39936698379939284, 0.6734754381366532, 0.0, 0.12739518123975463],
      tspan = (0.0, 0.06), surface_flux = (flux_hll, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [7.144325530681224e-17, 2.123397983547417e-16, 5.061138912500049e-16, 3.6588423152083e-17, 8.449816179702522e-15, 3.9171737639099993e-16, 2.445565690318772e-16, 3.6588423152083e-17, 9.971153407737885e-17],
      linf = [2.220446049250313e-16, 8.465450562766819e-16, 1.8318679906315083e-15, 1.1102230246251565e-16, 1.4210854715202004e-14, 8.881784197001252e-16, 4.440892098500626e-16, 1.1102230246251565e-16, 4.779017148551244e-16],
      maxiters = 1,
      initial_condition = initial_condition_constant,
      atol = 2.0e-13)
  end

  @trixi_testset "elixir_mhd_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
      l2   = [1.2428577757311665, 1.7993079220010408, 1.6900088072938513, 0.0, 2.262655937452063, 0.21271823802469103, 0.23326763151925287, 0.0, 0.0034077510132434945],
      linf = [10.399670938190301, 14.061104237417954, 15.555424931057713, 0.0, 16.71324802175841, 1.3218953536100904, 1.4137334529340777, 0.0, 0.0882662999481618],
      tspan = (0.0, 0.05))
  end

  @trixi_testset "elixir_mhd_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_blast_wave.jl"),
      l2   = [0.17569823349539196, 3.853292514951796, 2.474036084054808, 0.0, 355.36545703316915, 2.3525087027248084, 1.394705056983077, 0.0, 0.029910308624236454],
      linf = [1.5831869462980066, 44.20237543674303, 12.866364462538552, 0.0, 2237.8614138686, 13.066894956593798, 8.984965484247244, 0.0, 0.5226498664960756],
      tspan = (0.0, 0.003))
  end
end

end # module


text:
module AlgebraicNumbers

using Nemo
import PolynomialRoots

export AlgebraicNumber
export *,+,-,/,^,root,==,inv
export sqrt,cbrt
export exp_alg,cos_alg,sin_alg
export log_alg,acos_alg,asin_alg
export pow2

include("algebraic.jl")
include("trig.jl")
include("promote.jl")
include("newton.jl")
include("inv_totient.jl")

end


text:
const allowed_res_name_non_alphanumeric_chars = ['_'];

#
# An String representing the Readeble name of the result repo
#
struct ResName
    res_name::AbstractString
    function ResName(resname::AbstractString)
        validate_res_name(resname);
        new(resname);
    end
end

Base.string(rn::ResName) = string(rn.res_name);


text:
using Plots
f1(x,L)=exp(-x/L)
L=2
x=collect(range(1, 10, length=1000))
y=collect(range(1, 101, length=1000))
phi=( (f1.(x,L*0.5)*20).+8 ) .- ((f1.(x,L*0.5)*5).^2.5 );
plot( phi ,-y*40,label=false )

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
abstract type AbstractRangeNumber end

Number(rn::AbstractRangeNumber) = rn.v

#struct NoRangeNumber{T<:Number} <: AbstractRangeNumber
#    v::T
#    function NoRangeNumber(v::Number)
#        return new{typeof(v)}(v)
#    end
#end

"""
    RangeNumber(v, r)

# Example
a = RangeNumber(3, 1:7)
"""
struct RangeNumber{T<:Number, R} <: AbstractRangeNumber
    v::T

    function RangeNumber(v::Number, r::UnitRange)
        v in r || error("$v not in $r")
        return new{typeof(v), r}(v)
    end
end

RangeNumber{T, R}(n) where {T, R} = RangeNumber(n, R)

Base.:+(a::RangeNumber{T, R}, b::RangeNumber{T, R}) where {T, R} = RangeNumber(a.v+b.v, R)
Base.:+(a::Number, b::RangeNumber{T, R}) where {T, R} = RangeNumber(a+b.v, R)
Base.:+(a::RangeNumber{T, R}, b::Number) where {T, R} = RangeNumber(a.v+b, R)

Base.:(==)(a::RangeNumber{T, R}, b::RangeNumber{T, R}) where {T, R} = a.v == b.v
Base.:(==)(a::Number, b::RangeNumber{T, R}) where {T, R} = a == b.v
Base.:(==)(a::RangeNumber{T, R}, b::Number) where {T, R} = a.v == b


text:
(var"##MTIIPVar#390", var"##MTKArg#388")->begin
        @inbounds begin
                begin
                    (ModelingToolkit.fill_array_with_zero!)(var"##MTIIPVar#390")
                    let (x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₁₇, x₁₈, x₁₉, x₂₀, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆, x₂₇, x₂₈, x₂₉, x₃₀, x₃₁, x₃₂) = (var"##MTKArg#388"[1], var"##MTKArg#388"[2], var"##MTKArg#388"[3], var"##MTKArg#388"[4], var"##MTKArg#388"[5], var"##MTKArg#388"[6], var"##MTKArg#388"[7], var"##MTKArg#388"[8], var"##MTKArg#388"[9], var"##MTKArg#388"[10], var"##MTKArg#388"[11], var"##MTKArg#388"[12], var"##MTKArg#388"[13], var"##MTKArg#388"[14], var"##MTKArg#388"[15], var"##MTKArg#388"[16], var"##MTKArg#388"[17], var"##MTKArg#388"[18], var"##MTKArg#388"[19], var"##MTKArg#388"[20], var"##MTKArg#388"[21], var"##MTKArg#388"[22], var"##MTKArg#388"[23], var"##MTKArg#388"[24], var"##MTKArg#388"[25], var"##MTKArg#388"[26], var"##MTKArg#388"[27], var"##MTKArg#388"[28], var"##MTKArg#388"[29], var"##MTKArg#388"[30], var"##MTKArg#388"[31], var"##MTKArg#388"[32])
                        var"##MTIIPVar#390"[1] = (*)((*)(1.0, (+)((*)((*)(1, x₁), 1), (*)((*)(1, x₁₂), 1))), 1.0)
                        var"##MTIIPVar#390"[2] = (*)((*)(1.0, (*)((*)((+)((+)((+)(x₂₁, false), 0), 0), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[3] = (*)((*)(1.0, (*)((*)((+)((+)((+)(x₂₂, false), 0), 0), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[4] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(1, x₂₉)), 0), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[5] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₃, x₂₉)), 0), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[6] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₄, x₂₉)), 0), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[7] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[8] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[9] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[10] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[11] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), 1)), 1.0)
                        var"##MTIIPVar#390"[13] = (*)((*)(1.0, (+)((*)((*)(1, x₂), 1), (*)((*)((+)((+)((+)(x₂₁, false), 0), 0), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[14] = (*)((*)(1.0, (*)((*)((+)((+)((+)(x₂₂, false), 0), 0), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[15] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(1, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[16] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₃, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[17] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₄, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[18] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[19] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[20] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[21] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[22] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₁, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[25] = (*)((*)(1.0, (+)((*)((*)(1, x₃), 1), (*)((*)((+)((+)((+)(x₂₂, false), 0), 0), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[26] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(1, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[27] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₃, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[28] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₄, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[29] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[30] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[31] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[32] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[33] = (*)((*)(1.0, (*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)((+)(x₂₂, false), 0), 0))), 1.0)
                        var"##MTIIPVar#390"[37] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₄), 1), (*)((*)((+)((+)(0, (*)(1, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0))), (*)((*)(1, x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[38] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₃, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0)), (*)((*)((+)((+)((+)(x₂₃, false), 0), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[39] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₄, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0)), (*)((*)((+)((+)((+)(x₂₄, false), 0), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[40] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0)), (*)((*)((+)((+)((+)(x₂₅, false), 0), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[41] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₁₅), 1), (*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0))), (*)((*)((+)((+)(0, (*)(1, x₃₂)), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[42] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₂)), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[43] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[44] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), 1))), 1.0)
                        var"##MTIIPVar#390"[49] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₅), 1), (*)((*)((+)((+)(0, (*)(x₂₃, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0))), (*)((*)((+)((+)((+)(x₂₃, false), 0), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[50] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₄, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0)), (*)((*)((+)((+)((+)(x₂₄, false), 0), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[51] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₁₆), 1), (*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0))), (*)((*)((+)((+)((+)(x₂₅, false), 0), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[52] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(1, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[53] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₁₇), 1), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0))), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[54] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[55] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₃, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₃, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[61] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₆), 1), (*)((*)((+)((+)(0, (*)(x₂₄, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(x₂₄, x₂₉)), 0))), (*)((*)((+)((+)((+)(x₂₄, false), 0), 0), x₁₃), (+)((+)((+)(x₂₄, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[62] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(x₂₄, x₂₉)), 0)), (*)((*)((+)((+)((+)(x₂₅, false), 0), 0), x₁₃), (+)((+)((+)(x₂₄, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[63] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₄, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(1, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₄, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[64] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₄, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₄, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[65] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₁₈), 1), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₄, x₂₉)), 0))), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₄, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[66] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₄, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₄, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[73] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₇), 1), (*)((*)((+)((+)(0, (*)(x₂₅, x₂₉)), 0), x₁₂), (+)((+)(0, (*)(x₂₅, x₂₉)), 0))), (*)((*)((+)((+)((+)(x₂₅, false), 0), 0), x₁₃), (+)((+)((+)(x₂₅, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[74] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₅, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(1, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₅, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[75] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₅, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₅, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[76] = (*)((*)(1.0, (+)((*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₅, x₂₉)), 0)), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₅, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[77] = (*)((*)(1.0, (+)((+)((*)((*)(1, x₁₉), 1), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₅, x₂₉)), 0))), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)((+)(x₂₅, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[85] = (*)((*)(1.0, (+)((+)((+)((*)((*)(1, x₈), 1), (*)((*)((+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉)))), (*)((*)((+)((+)(0, (*)(1, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(1, x₃₂)), 0))), (*)((*)(1, x₁₄), 1))), 1.0)
                        var"##MTIIPVar#390"[86] = (*)((*)(1.0, (+)((+)((*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉))), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(1, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₆, false), 0), 0), x₁₄), 1))), 1.0)
                        var"##MTIIPVar#390"[87] = (*)((*)(1.0, (+)((+)((*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉))), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(1, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₇, false), 0), 0), x₁₄), 1))), 1.0)
                        var"##MTIIPVar#390"[88] = (*)((*)(1.0, (+)((+)((*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(1, x₃₀)), (*)((*)(1, x₃₂), x₂₉))), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(1, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₈, false), 0), 0), x₁₄), 1))), 1.0)
                        var"##MTIIPVar#390"[97] = (*)((*)(1.0, (+)((+)((+)((*)((*)(1, x₉), 1), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)))), (*)((*)((+)((+)(0, (*)(x₂₆, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(x₂₆, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₆, false), 0), 0), x₁₄), (+)((+)((+)(x₂₆, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[98] = (*)((*)(1.0, (+)((+)((*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉))), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(x₂₆, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₇, false), 0), 0), x₁₄), (+)((+)((+)(x₂₆, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[99] = (*)((*)(1.0, (+)((+)((+)((*)((*)(1, x₂₀), 1), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₆, x₃₀)), (*)((*)(x₂₆, x₃₂), x₂₉)))), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(x₂₆, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₈, false), 0), 0), x₁₄), (+)((+)((+)(x₂₆, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[109] = (*)((*)(1.0, (+)((+)((+)((*)((*)(1, x₁₀), 1), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉)))), (*)((*)((+)((+)(0, (*)(x₂₇, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(x₂₇, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₇, false), 0), 0), x₁₄), (+)((+)((+)(x₂₇, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[110] = (*)((*)(1.0, (+)((+)((*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₇, x₃₀)), (*)((*)(x₂₇, x₃₂), x₂₉))), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(x₂₇, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₈, false), 0), 0), x₁₄), (+)((+)((+)(x₂₇, false), 0), 0)))), 1.0)
                        var"##MTIIPVar#390"[121] = (*)((*)(1.0, (+)((+)((+)((*)((*)(1, x₁₁), 1), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)), x₁₂), (+)((+)(0, (*)(x₂₈, x₃₀)), (*)((*)(x₂₈, x₃₂), x₂₉)))), (*)((*)((+)((+)(0, (*)(x₂₈, x₃₂)), 0), x₁₃), (+)((+)(0, (*)(x₂₈, x₃₂)), 0))), (*)((*)((+)((+)((+)(x₂₈, false), 0), 0), x₁₄), (+)((+)((+)(x₂₈, false), 0), 0)))), 1.0)
                    end
                end
            end
        nothing
    end


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

const _DEFAULT_VALUE_AT_MARGIN = 0.1

function tolerance(
    x;
    bounds = (0, 0),
    margin = 0,
    sigmoid = quadratic,
    value_at_margin = _DEFAULT_VALUE_AT_MARGIN,
)
    lo, up = bounds
    lo <= up || throw(DomainError("Lower bound must be <= upper bound"))
    margin >= 0 || throw(DomainError("`margin` must be > 0"))
    _tolerance(x, lo, up, margin, sigmoid, value_at_margin)
end

function _tolerance(x::Real, lo, up, margin, sigmoid, vmargin)
    in_bounds = lo <= x <= up
    if iszero(margin)
        return ifelse(in_bounds, one(x), zero(x))
    else
        d = ifelse(x < lo, lo - x, x - up) / margin
        value = in_bounds ? one(x) :  sigmoid(d, vmargin)
        return convert(typeof(x), value)
    end
end

function _tolerance(xs::AbstractArray, args...)
    map(x -> _tolerance(x, args...), xs)
end

function quadratic(x::Number, value_at_1::Number)
    scale = sqrt(1 - value_at_1)
    scaled_x = scale * x
    ifelse(abs(scaled_x) < 1, 1 - scaled_x ^ 2, zero(typeof(x)))
end
quadratic(xs::AbstractArray, value_at_1::Number) = map(x -> quadratic(x, value_at_1), xs)
quadratic!(xs::AbstractArray, value_at_1::Number) = xs .= quadratic.(xs, value_at_1)

function _check_value_at_1(x::Number)
    0 < x < 1 || throw(DomainError("`value_at_1` must be in range [0, 1)"))
end
function _check_value_at_1_strict(x::Number)
    0 < x < 1 || throw(DomainError("`value_at_1` must be in range (0, 1)"))
end


text:
using AxisArrays
using JuMP
using MathOptInterface
using SCIP


text:
#Function to help with plotting wave functions
#24/07/2017

function getalpha(atom)

  if atom=="87Rb"
    alpha_c = 9.0760
  end

  return alpha_c
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")
include("EnergiesT.jl")
include("PotentialsT.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_MGC_cg.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("gen_dr.jl")
include("calc_strfact_v1.jl")
include("calc_ewald_v1.jl")

include("diag_lobpcg.jl")
include("KS_solve_SCF.jl")

include("andersonmix.jl")
include("KS_solve_SCF_andersonmix.jl")

include("pulaymix.jl")
include("KS_solve_SCF_pulaymix.jl")

include("get_ub_lb_lanczos.jl")
include("KS_solve_ChebySCF.jl")
include("chebyfilt.jl")
include("norm_matrix_induced.jl")

function test_main( Ns; method="SCF" )

    LatVecs = gen_lattice_sc(16.0)

    pw = PWGrid( Ns, LatVecs )

    Ω  = pw.Ω
    r  = pw.r
    G  = pw.gvec.G
    G2 = pw.gvec.G2
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    actual = Npoints/Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    Sf = calc_strfact( Xpos, G )

    E_nn = calc_ewald( pw, Xpos, Sf )

    Vg = zeros(ComplexF64,Npoints)
    prefactor = -4*pi/Ω
    for ig=2:Npoints
        Vg[ig] = prefactor/G2[ig]
    end
    V_ionic = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    #
    Nstates = 1
    Focc = [1.0]

    if method == "CG"
        psi, Energies, Potentials = KS_solve_Emin_cg(
            pw, V_ionic, Focc, Nstates, NiterMax=1000, E_NN=E_nn )
        #
        Y = ortho_gram_schmidt(psi)
        mu = Y' * op_H( pw, Potentials, Y )
        evals, evecs = eigen(mu)
        psi = Y*evecs
    elseif method == "SCF"
        Energies, Potentials, psi, evals = KS_solve_SCF( pw, V_ionic, Focc, Nstates, E_NN=E_nn )
    elseif method == "ChebySCF"
        Energies, Potentials, psi, evals = KS_solve_ChebySCF( pw, V_ionic, Focc, Nstates, β=0.8, E_NN=E_nn )
    elseif method == "SCF_andersonmix"
        Energies, Potentials, psi, evals = KS_solve_SCF_andersonmix( pw, V_ionic, Focc, Nstates, β=0.5, E_NN=E_nn )
    elseif method == "SCF_pulaymix"
        Energies, Potentials, psi, evals = KS_solve_SCF_pulaymix( pw, V_ionic, Focc, Nstates, β=0.5, E_NN=E_nn)
    else
        println("ERROR: unknown method: ", method)
        exit()
    end

    for st = 1:Nstates
        @printf("State # %d, Energy = %18.10f\n", st, real(evals[st]))
    end

    print_Energies(Energies)

end

@time test_main( [66,66,66], method="CG" )
#@time test_main( [64,64,64], method="SCF" )
#@time test_main( [64,64,64], method="ChebySCF" )
#@time test_main( [64, 64, 64], method="SCF_andersonmix" )
#@time test_main( [64, 64, 64], method="SCF_pulaymix")


text:
## date handling functions
##------------------------

function getLastDayOfMonth(dt::Date)
    isendofmonth = x -> Dates.lastdayofmonth(x) == x
    Dates.tonext(isendofmonth, dt)
end

function getLastDayOfQuarter(dt::Date)
    quart = Dates.quarterofyear(dt)
    mon, day = 3, 31
    if quart == 2
        mon, day = 6, 30
    elseif quart == 3
        mon, day = 9, 30
    elseif quart == 4
        mon, day = 12, 31
    end
    newDate = Date(year(dt), mon, day)
end

function formatDate(dat::ASCIIString)
    if length(dat) == 4 # yearly data
        newDat = Date(string(dat, "-12-31"))
    else
        # monthly data
        if dat[5] == 'M'
            someDate = Date(string(dat[1:4], "-", dat[6:7], "-01"))
            newDat = getLastDayOfMonth(someDate)

            ## quaterly data
        elseif dat[5:6] == "Q1"
            someDate = Date(string(dat[1:4], "-01-01"))
            newDat = getLastDayOfQuarter(someDate)
        elseif dat[5:6] == "Q2"
            someDate = Date(string(dat[1:4], "-04-01"))
            newDat = getLastDayOfQuarter(someDate)
        elseif dat[5:6] == "Q3"
            someDate = Date(string(dat[1:4], "-07-01"))
            newDat = getLastDayOfQuarter(someDate)
        elseif dat[5:6] == "Q4"
            someDate = Date(string(dat[1:4], "-10-01"))
            newDat = getLastDayOfQuarter(someDate)
        end
    end
    return newDat
end

## database style select functions
##--------------------------------

import Base.get
function get(df::DataFrame, sb::Symbol, f::Function)
    ## pick observation based on single entry in row
    nObs, nVars = size(df)
    logInd = Array(Bool, nObs)
    for ii=1:nObs
        logInd[ii] = f(df[sb][ii])
    end
    return df[logInd, :]
end

function get(df::DataFrame, f::Function)
    ## pick observation based on complete row
    nObs, nVars = size(df)
    logInd = Array(Bool, nObs)
    for ii=1:nObs
        logInd[ii] = f(df[ii, :])
    end
    return df[logInd, :]
end

## conversion functions
##---------------------

function tofloat(f::String)
     try
         return float(f)
     catch
         return NA
     end
end

function tofloat(Void)
    return NA
end

function clean_and_push!(da::DataArray, Void)
    push!(da, NA)
end

function clean_and_push!(da::DataArray, x::String)
    if (x == "NA") | (x == "")
        push!(da, NA)
    else
        push!(da, x)
    end
end

## function clean_entry(x::Union(String,Nothing))
##     if typeof(x) == Nothing # also equals Void
##         return "NA"
##     else
##         return x
##     end
## end

## function clean_append!(vals::Union(Array{UTF8String,1},Array{ASCIIString,1}),val::Union(UTF8String,ASCIIString,Nothing))
##     append!(vals,[clean_entry(val)])
## end


text:
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of Intel Corporation nor the names of its contributors
#       may be used to endorse or promote products derived from this software
#       without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

using FactCheck

include("fusion.jl")
include("test_snapshot.jl")
include("transforms/test_fixed_dims.jl")
include("transforms/test_neuron.jl")
include("transforms/test_tree_cleaners.jl")
include("stdlib/test_accuracy.jl")
include("stdlib/test_concat.jl")
include("stdlib/test_convolution.jl")
include("stdlib/test_dropout.jl")
include("stdlib/test_embed_id.jl")
include("stdlib/test_gru.jl")
include("stdlib/test_hdf5-data.jl")
include("stdlib/test_inner-product.jl")
include("stdlib/test_lstm.jl")
include("stdlib/test_math.jl")
include("stdlib/test_memory-data-layer.jl")
include("stdlib/test_pooling.jl")
include("stdlib/test_relu.jl")
include("stdlib/test_reshape.jl")
include("stdlib/test_rnn.jl")
include("stdlib/test_softmax.jl")
include("stdlib/test_tanh.jl")
include("stdlib/test_transform.jl")

FactCheck.exitstatus()


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.generator" begin

@testset "pathgraph" begin
    @test_throws DomainError pathgraph(1)
    P2 = pathgraph(2)
    @test nodecount(P2) == edgecount(P2) + 1
    P100 = pathgraph(100)
    @test nodecount(P100) == edgecount(P100) + 1
    @test distance(P100, 1, 100) == 99
end

@testset "cyclegraph" begin
    @test_throws DomainError cyclegraph(2)
    C3 = cyclegraph(3)
    @test nodecount(C3) == edgecount(C3)
    C100 = cyclegraph(100)
    @test nodecount(C100) == edgecount(C100)
    @test distance(C100, 1, 100) == 1
    @test distance(C100, 1, 50) == 49
end

@testset "completebipartite" begin
    @test_throws DomainError completebipartite(0, 1)
    @test_throws DomainError completebipartite(1, 0)
    K1_1 = completebipartite(1, 1)
    @test nodecount(K1_1) == 2
    @test edgecount(K1_1) == 1
    K3_3 = completebipartite(3, 3)
    @test nodecount(K3_3) == 6
    @test edgecount(K3_3) == 9
end

@testset "completegraph" begin
    @test_throws DomainError completegraph(-1)
    K0 = completegraph(0)
    @test nodecount(K0) == 0
    @test edgecount(K0) == 0
    K1 = completegraph(1)
    @test nodecount(K1) == 1
    @test edgecount(K1) == 0
    K2 = completegraph(2)
    @test nodecount(K2) == 2
    @test edgecount(K2) == 1
    K3 = completegraph(3)
    @test nodecount(K3) == 3
    @test edgecount(K3) == 3
    K5 = completegraph(5)
    @test nodecount(K5) == 5
    @test edgecount(K5) == 10
    K20 = completegraph(20)
    @test nodecount(K20) == 20
    @test edgecount(K20) == 190
end

@testset "laddergraph" begin
    @test_throws DomainError laddergraph(0)
    L3 = laddergraph(3)
    @test edgecount(L3) == 7 # 3n - 2
    L100 = laddergraph(100)
    @test edgecount(L100) == 298 # 3n - 2
    @test distance(L100, 1, 100) == 50
end

@testset "circularladder" begin
    @test_throws DomainError circularladder(2)
    CL3 = circularladder(3)
    @test edgecount(CL3) == 9 # 3n
    CL100 = circularladder(100)
    @test edgecount(CL100) == 300 # 3n
    @test distance(CL100, 1, 50) == 25
    @test distance(CL100, 1, 200) == 2
end

@testset "moebiusladder" begin
    @test_throws DomainError moebiusladder(2)
    ML3 = moebiusladder(3)
    @test edgecount(ML3) == 9 # 3n
    ML100 = moebiusladder(100)
    @test edgecount(ML100) == 300 # 3n
    @test distance(ML100, 1, 50) == 25
    @test distance(ML100, 1, 200) == 1
end

@testset "squaregrid" begin
    @test_throws DomainError squaregrid(1, 0)
    @test_throws DomainError squaregrid(0, 1)
    single = squaregrid(1, 1)
    @test nodecount(single) == 1
    @test edgecount(single) == 0
    hundred = squaregrid(10, 10)
    @test nodecount(hundred) == 100
    @test edgecount(hundred) == 180
end

@testset "petersen" begin
    @test_throws DomainError generalizedpetersen(2, 1)
    @test_throws DomainError generalizedpetersen(3, 2)
    p = petersengraph()
    @test edgecount(p) == 15 # 3n
    d = dodecahedralgraph()
    @test edgecount(d) == 30 # 3n
end

end # graph.generator


text:
using OREnvironment
using ORInterface
using Test

@testset "VariableDomain" begin
  domain = OREnvironment.VariableDomain(1.0, 5.0)
  @test OREnvironment.get_lb(domain) == 1.0
  @test OREnvironment.get_ub(domain) == 5.0
  @test OREnvironment.get_middle_point(domain) == 3.0

  OREnvironment.set_lb!(domain, 2.0)
  OREnvironment.set_ub!(domain, 3.0)
  @test OREnvironment.get_middle_point(domain) == 2.5
  @test OREnvironment.get_lb(domain) == 2.0
  @test OREnvironment.get_ub(domain) == 3.0
  @test OREnvironment.is_value_within_the_domain(domain, 2)
  @test OREnvironment.is_value_within_the_domain(domain, 2.5)
  @test OREnvironment.is_value_within_the_domain(domain, 3)
  @test OREnvironment.is_value_within_the_domain(domain, 4.0) == false
  @test OREnvironment.is_value_within_the_domain(domain, 1) == false
end

function constructorSolution(numConstraints=2)
    status = OREnvironment.constructStatus(numConstraints);
    Tvariables = Int; sizeArray = 3;
    sol = OREnvironment.constructSolution(:FixedLengthArray, (Tvariables, sizeArray, status));
    return sol;
end

@testset "is_solution_within_bounds" begin
  d1 = OREnvironment.VariableDomain(1.0, 5.0)
  d2 = OREnvironment.VariableDomain(2.0, 4.0)
  d3 = OREnvironment.VariableDomain(3.0, 7.0)
  domain = [d1, d2, d3]
  sol = constructorSolution()
  @test OREnvironment.is_solution_within_bounds(domain, sol) == false

  OREnvironment.add_solution!(sol, 1, 2);
  OREnvironment.add_solution!(sol, 2, 3);
  OREnvironment.add_solution!(sol, 3, 6);
  @test OREnvironment.is_solution_within_bounds(domain, sol)
end

@testset "building problem" begin
  cost = collect(1.0:6.0);
  variables1 = [1, 3, 4, 6];
  variables2 = [1, 3, 5, 6];
  coefs1 = [2.3, 3.2, 3.1, 12.34];
  coefs2 = coefs1 .+ 1.0;
  constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
  constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
  constraints = [constraint1, constraint2];
  variablesConstraints = [[1,2], Int[], [1,2], [1], [2], [1,2]];
  domain = [OREnvironment.VariableDomain(0.0, 5.0) for i in 1:6]
  p = OREnvironment.DefaultProblem(cost, constraints, variablesConstraints, :max, domain);

  @test OREnvironment.get_number_of_constraints(p) == 2;
  @test OREnvironment.get_number_of_variables(p) == 6;

  @test OREnvironment.get_constraints(p) === constraints;
  for i in 1:length(cost)
    @test OREnvironment.get_cost(p, i) == cost[i];
    @test OREnvironment.get_lb_variable(p, i) == 0.0
    @test OREnvironment.get_ub_variable(p, i) == 5.0
    @test OREnvironment.get_middle_point_variable_domain(p, i) == 2.5
    @test OREnvironment.is_value_within_the_domain(p, i, 4.0) == true
  end
  for i in 1:length(variablesConstraints)
    @test OREnvironment.get_constraints_of_variable(p, i) == variablesConstraints[i];
  end
  for i in 1:length(cost)
    OREnvironment.set_cost!(p, i, 12.0*i);
    OREnvironment.set_lb_variable!(p, i, 1.0)
    OREnvironment.set_ub_variable!(p, i, 3.0)
    @test OREnvironment.is_value_within_the_domain(p, i, 4.0) == false
  end
  for i in 1:length(cost)
    @test OREnvironment.get_cost(p, i) == 12.0*i;
    @test OREnvironment.get_lb_variable(p, i) == 1.0
    @test OREnvironment.get_ub_variable(p, i) == 3.0
  end
  @test OREnvironment.get_obj_sense(p) == :max;
  OREnvironment.set_obj_sense!(p, :min);
  @test OREnvironment.get_obj_sense(p) == :min;

  @test OREnvironment.get_constraint(p, 1) === constraint1;
  @test OREnvironment.get_constraint(p, 1) == constraint1;
  @test OREnvironment.get_constraint(p, 2) === constraint2;
  @test OREnvironment.get_constraint(p, 2) == constraint2;

  @test OREnvironment.is_variable(OREnvironment.get_constraint(p,1), 1) == true;
  @test OREnvironment.is_variable(OREnvironment.get_constraint(p,1), 2) == false;
  @test OREnvironment.is_variable(OREnvironment.get_constraint(p,2), 4) == false;
  @test OREnvironment.is_variable(OREnvironment.get_constraint(p,2), 5) == true;
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,1)) == 15.0;
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,2)) == 9.0;
  @test OREnvironment.get_type(OREnvironment.get_constraint(p,2)) == :lessOrEq;
  @test OREnvironment.get_coefficient(OREnvironment.get_constraint(p,2), 5) == 4.1;
  @test OREnvironment.get_coefficient(OREnvironment.get_constraint(p,2), 4) == 0.0;
end

@testset "building problem as in practice" begin
  # read constraints one by one in a for loop and save them
  # obtain the variables constraints with the method
  # build the problem
  cost = rand(6);
  data_coefs = [rand(3), rand(2), rand(4)];
  data_vars  = [[1, 2, 3], [3, 6], [3, 4, 5, 6]];
  data_rhs = [12.3, 9.0, 33.0];
  data_symbol = [:lessOrEq, :equal, :greaterOrEq];
  # simulating the reading of constraints
  # check function read_constraints in ./test/Constraints.jl for actually
  # reading from a file!!!
  constraints = Array{Constraint, 1}();
  for i in 1:3
    vars = data_vars[i]; # in practice we read from file
    coefs = data_coefs[i];
    rhs = data_rhs[i];
    symbol = data_symbol[i];
    push!(constraints, OREnvironment.constructConstraint(rhs, symbol, vars, coefs));
  end
  domain = [OREnvironment.VariableDomain(0.0, 5.0) for i in 1:6]
  p = OREnvironment.constructProblem(cost, constraints, :max, domain);

  # test that everything works
  variablesConstraints = OREnvironment.get_relationship_variables_constraints(constraints, 6);
  @test length(constraints) == 3;
  @test length(p.constraints) == 3;
  for i in 1:3
    @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,i)) == data_rhs[i];
    @test OREnvironment.get_type(OREnvironment.get_constraint(p,i)) == data_symbol[i];
    @test OREnvironment.get_constraint(p, i) == constraints[i];
    @test OREnvironment.get_constraint(p, i) === constraints[i];
  end
  for i in 1:6
    @test OREnvironment.get_cost(p,i) == cost[i];
    @test OREnvironment.get_constraints_of_variable(p, i) == variablesConstraints[i];
  end

  # simulating adding more constraints one by one
  constraint4 = OREnvironment.constructConstraint(3.9, :equal, [5, 6], [1.4, 11.0]);
  OREnvironment.add_constraint!(p, constraint4);
  @test length(p.constraints) == 4;
  @test length(constraints) == 4; # this is because p has a pointer to constraints!!!
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,4)) == 3.9;
  @test OREnvironment.get_type(OREnvironment.get_constraint(p,4)) == :equal;
  @test OREnvironment.get_constraint(p, 4) == constraint4;
  @test OREnvironment.get_constraint(p, 4) === constraint4;
  variablesConstraints = OREnvironment.get_relationship_variables_constraints(constraints, 6);
  for i in 1:6
    @test OREnvironment.get_constraints_of_variable(p, i) == variablesConstraints[i];
  end

  constraint5 = OREnvironment.constructConstraint(5.5, :greaterOrEq, [3, 6], [2.4, 11.0]);
  OREnvironment.add_constraint!(p, constraint5);
  @test length(p.constraints) == 5;
  @test length(constraints) == 5; # this is because p has a pointer to constraints!!!
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,5)) == 5.5;
  @test OREnvironment.get_type(OREnvironment.get_constraint(p,5)) == :greaterOrEq;
  @test OREnvironment.get_constraint(p, 5) == constraint5;
  @test OREnvironment.get_constraint(p, 5) === constraint5;
  variablesConstraints = OREnvironment.get_relationship_variables_constraints(constraints, 6);
  for i in 1:6
    @test OREnvironment.get_constraints_of_variable(p, i) == variablesConstraints[i];
  end

  # simulating adding more constraints as a vector
  constraint6 = OREnvironment.constructConstraint(6.9, :equal, [5, 6], [1.4, 11.0]);
  constraint7 = OREnvironment.constructConstraint(7.9, :lessOrEq, [5, 6], [1.4, 11.0]);
  newConstraints = [constraint6, constraint7];
  OREnvironment.add_constraint!(p, newConstraints);
  @test length(p.constraints) == 7;
  @test length(constraints) == 7; # this is because p has a pointer to constraints!!!
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,6)) == 6.9;
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,7)) == 7.9;
  @test OREnvironment.get_type(OREnvironment.get_constraint(p,6)) == :equal;
  @test OREnvironment.get_type(OREnvironment.get_constraint(p,7)) == :lessOrEq;
  @test OREnvironment.get_constraint(p, 6) == constraint6;
  @test OREnvironment.get_constraint(p, 6) === constraint6;
  @test OREnvironment.get_constraint(p, 7) == constraint7;
  @test OREnvironment.get_constraint(p, 7) === constraint7;
  variablesConstraints = OREnvironment.get_relationship_variables_constraints(constraints, 6);
  for i in 1:6
    @test OREnvironment.get_constraints_of_variable(p, i) == variablesConstraints[i];
  end
end

function build_problem()
  cost = collect(1.0:6.0);
  data_coefs = [[1.1, 1.2, 1.3], [2.1, 2.2], [3.1, 3.2, 3.3, 3.4]];
  data_vars  = [[1, 2, 3], [3, 6], [3, 4, 5, 6]];
  data_rhs = [12.3, 9.0, 33.0];
  data_symbol = [:lessOrEq, :equal, :greaterOrEq];
  # simulating the reading of constraints
  constraints = Array{Constraint, 1}();
  for i in 1:3
    vars = data_vars[i]; # in practice we read from file
    coefs = data_coefs[i];
    rhs = data_rhs[i];
    symbol = data_symbol[i];
    push!(constraints, OREnvironment.constructConstraint(rhs, symbol, vars, coefs));
  end
  # variablesConstraints = OREnvironment.get_relationship_variables_constraints(constraints, 6);
  domain = [OREnvironment.VariableDomain(0.0, 5.0) for i in 1:6]
  p = OREnvironment.constructProblem(cost, constraints, :max, domain);
  return p;
end

@testset "scope of problem" begin
  # test that there is no issue with the scope of problem!
  p = build_problem();
  data_rhs = [12.3, 9.0, 33.0];
  data_symbol = [:lessOrEq, :equal, :greaterOrEq];
  @test length(p.constraints) == 3;
  for i in 1:3
    @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,i)) == data_rhs[i];
    @test OREnvironment.get_type(OREnvironment.get_constraint(p,i)) == data_symbol[i];
  end
  @test OREnvironment.get_coefficient(p.constraints[1], 1) == 1.1;
  @test OREnvironment.get_coefficient(p.constraints[1], 2) == 1.2;
  @test OREnvironment.get_coefficient(p.constraints[1], 3) == 1.3;
  @test OREnvironment.get_coefficient(p.constraints[1], 4) == 0.0;
  @test OREnvironment.get_coefficient(p.constraints[2], 3) == 2.1;
  @test OREnvironment.get_coefficient(p.constraints[2], 6) == 2.2;
  @test OREnvironment.get_coefficient(p.constraints[2], 5) == 0.0;
  @test OREnvironment.get_coefficient(p.constraints[3], 3) == 3.1;
  @test OREnvironment.get_coefficient(p.constraints[3], 4) == 3.2;
  @test OREnvironment.get_coefficient(p.constraints[3], 5) == 3.3;
  @test OREnvironment.get_coefficient(p.constraints[3], 6) == 3.4;
  for i in 1:6
    @test OREnvironment.get_cost(p,i) == i*1.0;
  end
  @test OREnvironment.get_constraints_of_variable(p, 1) == [1];
  @test OREnvironment.get_constraints_of_variable(p, 2) == [1];
  @test OREnvironment.get_constraints_of_variable(p, 3) == [1,2,3];
  @test OREnvironment.get_constraints_of_variable(p, 4) == [3];
  @test OREnvironment.get_constraints_of_variable(p, 5) == [3];
  @test OREnvironment.get_constraints_of_variable(p, 6) == [2,3];

  constraint4 = OREnvironment.constructConstraint(3.9, :equal, [5, 6], [1.4, 11.0]);
  OREnvironment.add_constraint!(p, constraint4);
  @test length(p.constraints) == 4;
  @test OREnvironment.get_rhs(OREnvironment.get_constraint(p,4)) == 3.9;
  @test OREnvironment.get_type(OREnvironment.get_constraint(p,4)) == :equal;
  @test OREnvironment.get_constraint(p, 4) == constraint4;
  @test OREnvironment.get_constraint(p, 4) === constraint4;
  @test OREnvironment.get_constraints_of_variable(p, 1) == [1];
  @test OREnvironment.get_constraints_of_variable(p, 2) == [1];
  @test OREnvironment.get_constraints_of_variable(p, 3) == [1,2,3];
  @test OREnvironment.get_constraints_of_variable(p, 4) == [3];
  @test OREnvironment.get_constraints_of_variable(p, 5) == [3, 4];
  @test OREnvironment.get_constraints_of_variable(p, 6) == [2,3,4];
end

@testset "is_feasible" begin
  typeVariables = Float64; numVariables = 6;  numConstraints = 2;
  status = OREnvironment.constructStatus(numConstraints);
  args = (typeVariables, numVariables, status);
  s = OREnvironment.constructSolution(:FixedLengthArray, args)
  cost = collect(1.0:6.0);
  variables1 = [1, 3, 4, 6];
  variables2 = [1, 3, 5, 6];
  coefs1 = [2.3, 3.2, 3.1, 12.34];
  coefs2 = coefs1 .+ 1.0;
  constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
  constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
  constraints = [constraint1, constraint2];
  domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
  p = OREnvironment.constructProblem(cost, constraints, :max, domain);
  @test OREnvironment.is_feasible(p,s)

  OREnvironment.add_solution!(s, 1, 1.2) # value 1.2 violates domain
  OREnvironment.update_constraint_consumption_and_feasibility!(s, p.constraints)
  @test OREnvironment.is_feasible(p,s) == false

  OREnvironment.add_solution!(s, 1, 0.5)
  OREnvironment.add_solution!(s, 2, 0.5)
  OREnvironment.add_solution!(s, 3, 0.5)
  OREnvironment.update_constraint_consumption_and_feasibility!(s, p.constraints)
  @test OREnvironment.is_feasible(p,s)

  OREnvironment.add_solution!(s,6,1.0) # we violate constraints
  OREnvironment.update_constraint_consumption_and_feasibility!(s, p.constraints)
  @test OREnvironment.is_feasible(p,s) == false
end

@testset "gap function" begin
  LB = 7.0;
  UB = 9.0;
  @test OREnvironment.gap(LB, UB, :min) == (UB-LB)/LB;
  @test OREnvironment.gap(LB, UB, :max) == (UB-LB)/UB;
end


text:
# Similar to JuMP/test/solvers.jl

using JuMP

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

mos = try_import(:MosekTools)
if mos
    mos_factory = with_optimizer(MosekTools.Mosek.Optimizer, QUIET=true)
else
    mos_factory = nothing
end
ismosek(factory) = factory === mos_factory
csd = try_import(:CSDP)
if csd
    csd_factory = with_optimizer(CSDP.Optimizer, printlevel=0)
else
    csd_factory = nothing
end
iscsdp(solver) = solver === csd_factory
sda = try_import(:SDPA)
if sda
    sda_factory = with_optimizer(SDPA.Optimizer)
else
    sda_factory = nothing
end
issdpa(solver) = solver === sda_factory
scs = false && try_import(:SCS) # It does not work
isscs(solver) = false
ipt = false && try_import(:Ipopt)

# Semidefinite solvers
sdp_factories = Any[]
mos && push!(sdp_factories, mos_factory)
csd && push!(sdp_factories, csd_factory)
sda && push!(sdp_factories, sda_factory)
#scs && push!(sdp_factories, scs_factory)

# Bilinear LP solvers
blp_factories = Any[]
ipt && push!(blp_factories, with_optimizer(Ipopt.Optimizer, print_level=0))


text:
"""
    find_subtours(n, X)

Compute sub-tours from tentative tour `X`.
"""
function find_subtours(n, X)
    (n, n) == size(X) || throw(DimensionMismatch(
        "n=$n but X has size $(size(X))."
    ))

    # List of all sub-tours in `X`
    # Empty if `X` is a tour
    sub_tours = Vector{Vector{Int}}()

    # `f[i]` is the city that follows `i` in the tour
    f = zeros(Int, n)
    @inbounds for j in 1:n
        for i in 1:n
            X[i, j] == 1.0 && (f[i] = j; break;)
        end
    end

    flags = falses(n)
    num_seen = 0
    while num_seen < n
        first = argmin(flags)  # first city in current tour
        flags[first] && break  # all cities have been visited
        subtour = [first]
        flags[first] = true
        num_seen += 1
        next = f[first]
        next == first && error("X[$first, $first] = $(X[first, first]).")
        while next != first
            push!(subtour, next)
            flags[next] = 1
            num_seen += 1
            next = f[next]
        end

        length(subtour) == n && return sub_tours  # X is a tour
        # @info "Found sub-tour: $subtour"
        push!(sub_tours, subtour)
    end

    @info "Found $(length(sub_tours)) sub-tours."
    return sub_tours

end


text:
using Unmarshal
using JSON
using Base.Test
import Base.==

#Tests for various type of composite structures, including Nullables
input = "{ \"bar\": { \"baz\": 17 }, \"foo\": 3.14 }"

immutable Bar
    baz::Int
end

immutable Foo
    bar::Bar
end

immutable Baz
    foo::Nullable{Float64}
    bar::Bar
end

immutable Qux
    baz::Nullable{String}
    bar::Bar
end

@test Unmarshal.unmarshal(Foo, JSON.parse(input)) === Foo(Bar(17))
@test Unmarshal.unmarshal(Baz, JSON.parse(input)) === Baz(Nullable(3.14), Bar(17))
@test Unmarshal.unmarshal(Qux, JSON.parse(input)) === Qux(Nullable{String}(),Bar(17))
@test_throws ArgumentError Unmarshal.unmarshal(Bar, JSON.parse(input))

#Test for structures of handling 1-D arrays
type StructOfArrays
        a1 :: Array{Float32, 1}
        a2 :: Array{Int, 1}
    end

function ==(A :: StructOfArrays, B :: StructOfArrays)
    A.a1 == B.a1 && A.a2 == B.a2
end

tmp = StructOfArrays([0,1,2], [1,2,3])
jstring = JSON.json(tmp)
@test Unmarshal.unmarshal(StructOfArrays, JSON.parse(jstring)) == tmp

#Test for handling 2-D arrays
type StructOfArrays2D
        a3 :: Array{Float64, 2}
        a4 :: Array{Int, 2}
    end

function ==(A :: StructOfArrays2D, B :: StructOfArrays2D)
    A.a3 == B.a3 && A.a4 == B.a4
end

tmp2 = StructOfArrays2D(ones(Float64, 2, 3), eye(Int, 2, 3))
jstring = JSON.json(tmp2)
@test Unmarshal.unmarshal(StructOfArrays2D, JSON.parse(jstring))  == tmp2

#Test for handling N-D arrays
tmp3 = randn(Float64, 2, 3, 4)
jstring = JSON.json(tmp3)
@test Unmarshal.unmarshal(Array{Float64, 3}, JSON.parse(jstring))  == tmp3

#Test for handling arrays of composite entities
tmp4 = Array{Array{Int,2}}(2)

tmp4[1] = ones(Int, 3, 4)
tmp4[2] = zeros(Int, 1, 2)
tmp4
jstring = JSON.json(tmp4)
@test Unmarshal.unmarshal(Array{Array{Int,2}}, JSON.parse(jstring)) == tmp4

# Test to check handling of complex numbers
tmp5 = zeros(Float32, 2) + 1im * ones(Float32, 2)
jstring = JSON.json(tmp5)
@test Unmarshal.unmarshal(Array{Complex{Float32}}, JSON.parse(jstring)) == tmp5

tmp6 = zeros(Float32, 2, 2) + 1im * ones(Float32, 2, 2)
jstring = JSON.json(tmp6)
@test Unmarshal.unmarshal(Array{Complex{Float32},2}, JSON.parse(jstring)) == tmp6

# Test to see handling of abstract types
type reconfigurable{T}
    x :: T
    y :: T
    z :: Int
end

function =={T1, T2}(A :: reconfigurable{T1}, B :: reconfigurable{T2})
    T1 == T2 && A.x == B.x && A.y == B.y && A.z == B.z
end

type higherlayer
    val :: reconfigurable
end

val = reconfigurable(1.0, 2.0, 3)
jstring = JSON.json(val)
@test Unmarshal.unmarshal(reconfigurable{Float64}, JSON.parse(jstring)) == reconfigurable{Float64}(1.0, 2.0, 3)

higher = higherlayer(val)
jstring = JSON.json(higher)
@test_throws ArgumentError Unmarshal.unmarshal(higherlayer, JSON.parse(jstring))

# Test string pass through
@test Unmarshal.unmarshal(String, JSON.parse(json("Test"))) == "Test"

# Test the verbose option
@test Unmarshal.unmarshal(Foo, JSON.parse(input), true) === Foo(Bar(17))
jstring = JSON.json(tmp3)
@test Unmarshal.unmarshal(Array{Float64, 3}, JSON.parse(jstring), true)  == tmp3

@test Unmarshal.unmarshal(String, JSON.parse(json("Test")), true) == "Test"

# Added test cases to attempt getting 100% code coverage
@test isequal(unmarshal(Nullable{Int64}, Void()), Nullable{Int64}())

@test_throws ArgumentError unmarshal(Nullable{Int64}, ones(Float64, 1))

# Test handling of Tuples
testTuples = ((1.0, 2.0, 3.0, 4.0), (2.0, 3.0))
jstring = JSON.json(testTuples)
@test Unmarshal.unmarshal(Tuple{Tuple{Float64}}, JSON.parse(jstring)) == testTuples
@test Unmarshal.unmarshal(Tuple{Array{Float64}}, JSON.parse(jstring)) == (([testElement...] for testElement in testTuples)...)
@test Unmarshal.unmarshal(Array{Tuple{Float64}}, JSON.parse(jstring)) == [testTuples...]
@test Unmarshal.unmarshal(Array{Array{Float64}}, JSON.parse(jstring)) == [([testElement...] for testElement in testTuples)...]
@test Unmarshal.unmarshal(Tuple{Tuple{Float64}}, JSON.parse(jstring), true) == testTuples

mutable struct DictTest
testDict::Dict{Int, String}
end

function ==(D1 :: DictTest, D2 :: DictTest)
    for iter in keys(D1.testDict)
        if !(D1.testDict[iter] == D2.testDict[iter])
          return false
        end
    end
    for iter in keys(D2.testDict)
        if !(D1.testDict[iter] == D2.testDict[iter])
          return false
        end
    end

    true
end

dictTest = DictTest(Dict{Int, String}(1 => "Test1", 2 => "Test2"))

#@show JSON.json(dictTest)
#@show JSON.parse(JSON.json(dictTest))
@test Unmarshal.unmarshal(DictTest, JSON.parse(JSON.json(dictTest)),true) == dictTest


text:
"""
Run a coarse thermal bubble simulation and save the output to NetCDF at the
10th time step. Then read back the output and test that it matches the model's
state.
"""
function run_thermal_bubble_netcdf_tests(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    closure = ConstantIsotropicDiffusivity(ν=4e-2, κ=4e-2)
    model = Model(architecture=arch, grid=grid, closure=closure)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    outputs = Dict("v"=>model.velocities.v,
                   "u"=>model.velocities.u,
                   "w"=>model.velocities.w,
                   "T"=>model.tracers.T,
                   "S"=>model.tracers.S)
    nc_writer = NetCDFOutputWriter(model, outputs,
                                   filename="dumptest.nc",
                                   frequency=10)
    push!(model.output_writers, nc_writer)

    xC_slice = 1:10
    xF_slice = 2:11
    yC_slice = 10:15
    yF_slice = 1
    zC_slice = 10
    zF_slice = 9:11
    nc_sliced_writer = NetCDFOutputWriter(model, outputs,
                                          filename="dumptestsliced.nc",
                                          frequency=10,
                                          xC=xC_slice, xF=xF_slice, yC=yC_slice,
                                          yF=yF_slice, zC=zC_slice, zF=zF_slice)
    push!(model.output_writers, nc_sliced_writer)

    time_step!(model, 10, Δt)

    close(nc_writer)
    close(nc_sliced_writer)

    u = read_output(nc_writer, "u")
    v = read_output(nc_writer, "v")
    w = read_output(nc_writer, "w")
    T = read_output(nc_writer, "T")
    S = read_output(nc_writer, "S")

    @test all(u .≈ Array(interiorparent(model.velocities.u)))
    @test all(v .≈ Array(interiorparent(model.velocities.v)))
    @test all(w .≈ Array(interiorparent(model.velocities.w)))
    @test all(T .≈ Array(interiorparent(model.tracers.T)))
    @test all(S .≈ Array(interiorparent(model.tracers.S)))

    u_sliced = read_output(nc_sliced_writer, "u")
    v_sliced = read_output(nc_sliced_writer, "v")
    w_sliced = read_output(nc_sliced_writer, "w")
    T_sliced = read_output(nc_sliced_writer, "T")
    S_sliced = read_output(nc_sliced_writer, "S")

    @test all(u_sliced .≈ Array(interiorparent(model.velocities.u))[xF_slice, yC_slice, zC_slice])
    @test all(v_sliced .≈ Array(interiorparent(model.velocities.v))[xC_slice, yF_slice, zC_slice])
    @test all(w_sliced .≈ Array(interiorparent(model.velocities.w))[xC_slice, yC_slice, zF_slice])
    @test all(T_sliced .≈ Array(interiorparent(model.tracers.T))[xC_slice, yC_slice, zC_slice])
    @test all(S_sliced .≈ Array(interiorparent(model.tracers.S))[xC_slice, yC_slice, zC_slice])
end

function run_jld2_file_splitting_tests(arch)
    model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)))

    u(model) = Array(model.velocities.u.data.parent)
    fields = Dict(:u => u)

    function fake_bc_init(file, model)
        file["boundary_conditions/fake"] = π
    end

    ow = JLD2OutputWriter(model, fields; dir=".", prefix="test", frequency=1,
                          init=fake_bc_init, including=[:grid],
                          max_filesize=200KiB, force=true)
    push!(model.output_writers, ow)

    # 531 KiB of output will be written which should get split into 3 files.
    time_step!(model, 10, 1)

    # Test that files has been split according to size as expected.
    @test filesize("test_part1.jld2") > 200KiB
    @test filesize("test_part2.jld2") > 200KiB
    @test filesize("test_part3.jld2") < 200KiB

    for n in string.(1:3)
        filename = "test_part" * n * ".jld2"
        jldopen(filename, "r") do file
            # Test to make sure all files contain structs from `including`.
            @test file["grid/Nx"] == 16

            # Test to make sure all files contain info from `init` function.
            @test file["boundary_conditions/fake"] == π
        end

        # Leave test directory clean.
        rm(filename)
    end
end

"""
Run two coarse rising thermal bubble simulations and make sure that when
restarting from a checkpoint, the restarted simulation matches the non-restarted
simulation numerically.
"""
function run_thermal_bubble_checkpointer_tests(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    closure = ConstantIsotropicDiffusivity(ν=4e-2, κ=4e-2)
    true_model = Model(architecture=arch, grid=grid, closure=closure)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    true_model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    checkpointed_model = deepcopy(true_model)

    time_step!(true_model, 9, Δt)

    checkpointer = Checkpointer(checkpointed_model; frequency=5, force=true)
    push!(checkpointed_model.output_writers, checkpointer)

    # Checkpoint should be saved as "checkpoint5.jld" after the 5th iteration.
    time_step!(checkpointed_model, 5, Δt)

    # Remove all knowledge of the checkpointed model.
    checkpointed_model = nothing

    restored_model = restore_from_checkpoint("checkpoint5.jld2")

    time_step!(restored_model, 4, Δt; init_with_euler=false)

    rm("checkpoint0.jld2", force=true)
    rm("checkpoint5.jld2", force=true)

    # Now the true_model and restored_model should be identical.
    @test all(restored_model.velocities.u.data      .≈ true_model.velocities.u.data)
    @test all(restored_model.velocities.v.data      .≈ true_model.velocities.v.data)
    @test all(restored_model.velocities.w.data      .≈ true_model.velocities.w.data)
    @test all(restored_model.tracers.T.data         .≈ true_model.tracers.T.data)
    @test all(restored_model.tracers.S.data         .≈ true_model.tracers.S.data)
    @test all(restored_model.timestepper.Gⁿ.u.data .≈ true_model.timestepper.Gⁿ.u.data)
    @test all(restored_model.timestepper.Gⁿ.v.data .≈ true_model.timestepper.Gⁿ.v.data)
    @test all(restored_model.timestepper.Gⁿ.w.data .≈ true_model.timestepper.Gⁿ.w.data)
    @test all(restored_model.timestepper.Gⁿ.T.data .≈ true_model.timestepper.Gⁿ.T.data)
    @test all(restored_model.timestepper.Gⁿ.S.data .≈ true_model.timestepper.Gⁿ.S.data)
end

@testset "Output writers" begin
    println("Testing output writers...")

    for arch in archs
         @testset "NetCDF [$(typeof(arch))]" begin
             println("  Testing NetCDF output writer [$(typeof(arch))]...")
             run_thermal_bubble_netcdf_tests(arch)
         end

        @testset "JLD2 [$(typeof(arch))]" begin
            println("  Testing JLD2 output writer [$(typeof(arch))]...")
            run_jld2_file_splitting_tests(arch)
        end

        @testset "Checkpointer [$(typeof(arch))]" begin
            println("  Testing Checkpointer [$(typeof(arch))]...")
            run_thermal_bubble_checkpointer_tests(arch)
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CBFlib_small_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CBFlib_small")
JLLWrappers.@generate_main_file("CBFlib_small", UUID("9b99e207-4e38-5950-946f-a3f2976e367a"))
end  # module CBFlib_small_jll


text:
# Uncomment the following to install the necessary packages for plotting.
# Pkg.add(["Conda", "PyCall"])
# using Conda
# Conda.add("cartopy")

using PyCall, PyPlot, ClimateBase
ccrs = pyimport("cartopy.crs")
LONLAT = ccrs.PlateCarree()
DEFPROJ = ccrs.Mollweide()

"""
    earthscatter(A::ClimArray, projection = ccrs.Mollweide(); kwargs...) → fig, ax, cb
Plot a `ClimArray` with space type `CoordinateSpace` as a scatter plot with the color
of the points being the value of `A` at these points. This requires that `A` has only
one dimension, the coordinates.

Keyword values are propagated to `scatter` and keys of interest are `cmap, vmin, vmax`.
"""
function earthscatter(A::ClimArray, projection = DEFPROJ; kwargs...)
    coords = dims(A, Coord)
    @assert length(dims(A)) == 1 "Input ClimArray must only have one `Coord` dimension."
    eqlon = [l[1] for l in coords]
    eqlat = [l[2] for l in coords]
    fig = figure()
    fig.tight_layout()
    ax = subplot(111, projection=projection)
    ax, cb = earthscatter!(ax, eqlon, eqlat, A.data, projection; kwargs...)
    if A.name ≠ Symbol("")
        ax.set_title(string(A.name))
    end
    return fig, ax, cb
end

function earthscatter!(ax, A, projection = DEFPROJ; kwargs...)
    coords = dims(A, Coord)
    lon = [l[1] for l in coords]
    lat = [l[2] for l in coords]
    earthscatter!(ax, lon, lat, A, projection; kwargs...)
end

function earthscatter!(ax, lon, lat, A, projection = DEFPROJ;
    coast = true, vmin = minimum(A), vmax = maximum(A), levels = 11,
    ticks = nothing, kwargs...)
    lon = Array(lon)
    lat = Array(lat)
    lvls = range(vmin, vmax, length = levels)
    cmap = matplotlib.cm.get_cmap(get(kwargs, :cmap, :viridis), length(lvls)-1)
    norm = matplotlib.colors.Normalize(vmin=lvls[1], vmax=lvls[end])
    cax, kw = matplotlib.colorbar.make_axes(ax,location="right",pad=0.02,shrink=0.8)
    cb = matplotlib.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm)

    ax.scatter(lon, lat; transform = LONLAT, c=A, kwargs..., cmap = cmap, vmin = vmin, vmax = vmax)
    coast && ax.coastlines()
    if projection == LONLAT
        ax.set_xticks(-180:60:180, crs=LONLAT)
        ax.set_yticks(-90:30:90, crs=LONLAT)
    end
    ax.gridlines(alpha = 0.25)
    cb.set_ticks(isnothing(ticks) ? lvls[1:max(1, levels÷5):end] : ticks)
    return ax, cb
end

"""
    earthsurface(A::ClimArray, projection = ccrs.Mollweide(); kwargs...) → fig, ax, cb
Plot a `ClimArray` with space type `OrthogonalSpace` as a surface plot with the color
of the surface being the value of `A`. This requires that `A` has exactly two dimensions,
`Lon, Lat`.

Keyword values are `coast=true` to enable coastlines and also `cmap, vmin, vmax, levels`.
"""
function earthsurface(A, projection = DEFPROJ;
    coast = true, kwargs...)
    fig = figure()
    fig.tight_layout()
    ax = subplot(111, projection=projection)
    ax, cb = earthsurface!(ax, A; kwargs...)
    if projection == LONLAT
        ax.set_xticks(-180:60:180, crs=LONLAT)
        ax.set_yticks(-90:30:90, crs=LONLAT)
    end
    if A.name ≠ Symbol("")
        ax.set_title(string(A.name))
    end
    return fig, ax, cb
end

function earthsurface!(ax, A::ClimArray; kwargs...)
    ax.gridlines(alpha = 0.25)
    lon, lat, s = _getwrappeddata(A)
    earthsurface!(ax, lon, lat, s; kwargs...)
end

function earthsurface!(ax, lon, lat, A;
    vmin = minimum(A), vmax = maximum(A), levels = 11,
    coast = true, specific_contour = nothing, kwargs...)

    lon, lat, A = _getwrappeddata(copy(Array(lon)), copy(Array(lat)), A)
    lvls = range(vmin, vmax, length = levels)

    # TODO: Add colorbar as a separate function
    lvls = range(vmin, vmax, length = levels)
    cmap = matplotlib.cm.get_cmap(get(kwargs, :cmap, :viridis), length(lvls)-1)
    # cmap.set_under("k")
    norm = matplotlib.colors.Normalize(vmin=lvls[1], vmax=lvls[end])
    cax, kw = matplotlib.colorbar.make_axes(ax,location="right",pad=0.02,shrink=0.8)
    cb = matplotlib.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm)
    coast && ax.coastlines(; color = "gray")

    ax.contourf(Array(lon), Array(lat), clamp.(A', vmin, vmax), lvls, cmap = cmap,
                transform=LONLAT, vmin=vmin, vmax = vmax)

    if specific_contour isa Array
        ax.contour(Array(lon), Array(lat), clamp.(A', vmin, vmax);
        levels = specific_contour, colors = "r", transform=LONLAT)
    end
    return ax, cb
end

function _getwrappeddata(A::ClimArray)
    lon = copy(Array(dims(A, Lon).val))
    lat = copy(Array(dims(A, Lat).val))
    s = copy(Array(A.data)) # assumes data have Lon x Lat dimensions
    _getwrappeddata(lon, lat, s)
end

function _getwrappeddata(lon, lat, s)
    if lon[end] == 359.5 && lon[1] == 0.5
        # add one more entry to have perfect overlap in plotting
        x = s[1, :]
        s = vcat(Array(s), x')
        push!(lon, 360.5)
    end
    return lon, lat, s
end


text:
FEM.linquadinterpmod()
using FEM.LinQuadInterpMod

module QuadTrigMod

using FEM
FEM.quadtriginterpmod()
using FEM.QuadTrigInterpMod

import FEM: createinterp, creategps, createstorage, get_field,
            Bmatrix, doftypes, get_ndofs, get_geotype, get_ref_area, dNdxmatrix

import FEM: AbstractMaterialStatus, AbstractElemStorage, AbstractFElement, FENode2
import FEM: Vertex6, Point2, GaussPoint2, Du, Dv, GeoTrig

export QuadTrig

type QuadTrigStorage <: AbstractElemStorage
    B::Matrix{Float64}
    DeBe::Matrix{Float64}
    Ke::Matrix{Float64}
    ɛ::Vector{Float64}
    f_int::Vector{Float64}
    u_field::Vector{Float64}
end

QuadTrigStorage() = QuadTrigStorage(zeros(4, 12), zeros(4,12),
                                  zeros(12,12), zeros(4), zeros(12), zeros(12))

type QuadTrig{T <: AbstractMaterialStatus} <: AbstractFElement{T}
    vertices::Vertex6
    gps::Vector{GaussPoint2}
    n::Int
    interp::QuadTrigInterp
    storage::QuadTrigStorage
    matstats::Vector{T}
    temp_matstats::Vector{T}
end
gausspoints(elem::QuadTrig) = elem.gps

# Constructor
function QuadTrig{T <: AbstractMaterialStatus}(vertices::Vertex6, n, interp::QuadTrigInterp,
                 storage::QuadTrigStorage, gps::Vector{GaussPoint2}, matstat::T)
    matstats = T[]
    temp_matstats = T[]
    for i in 1:length(gps)
        push!(matstats, copy(matstat))
        push!(temp_matstats, copy(matstat))
    end
    QuadTrig(vertices, gps, n, interp, storage, matstats, temp_matstats)
end

get_ref_area(::QuadTrig) = 0.5
get_geoelem(ele::QuadTrig) = GeoQTrig(ele.n, ele.vertices)
get_geotype(::QuadTrig) = GeoQTrig

createstorage(::Type{QuadTrig}) = QuadTrigStorage()
createinterp(::Type{QuadTrig}) = QuadTrigInterp()

function creategps(::Type{QuadTrig})
    p1 = 1/3
    p2 = 0.2
    p3 = 0.6
    [GaussPoint2(Point2(p1, p1), -0.281250000000000);
     GaussPoint2(Point2(p2, p3), 0.260416666666667);
     GaussPoint2(Point2(p2, p2), 0.260416666666667)
     GaussPoint2(Point2(p3, p2), 0.260416666666667)]
end

@inline function get_ndofs(::QuadTrig)
    return 12
end

# Avoids allocation (remove this?)
const QUADTRIG_DOFTYPES = [Du, Dv]
function doftypes(::QuadTrig, ::Int)
    return return QUADTRIG_DOFTYPES
end

function Bmatrix(elem::QuadTrig, gp::GaussPoint2, nodes::Vector{FENode2})
    dNdx = dNdxmatrix(elem.interp, gp.local_coords, elem.vertices, nodes)
    B = elem.storage.B
    for i in 1:6
        B[1, 2*i - 1] = dNdx[i, 1]
        B[2, 2*i]     = dNdx[i, 2]

        B[4, 2*i - 1] = dNdx[i, 2]
        B[4, 2*i]     = dNdx[i, 1]
    end
    return B
end

end # module


text:
import Base: +, -, *, &, |, xor, ==, ~

for op in (:+, :-, :&, :|, :xor, :(==))
    @eval begin
        function ($op)(a::Cenum{T}, b::Cenum{S}) where {T<:Integer,S<:Integer}
            N = promote_type(T, S)
            ($op)(N(a), N(b))
        end
        function ($op)(a::Cenum{T}, b::S) where {T<:Integer,S<:Integer}
            N = promote_type(T, S)
            ($op)(N(a), N(b))
        end
        function ($op)(a::T, b::Cenum{S}) where {T<:Integer,S<:Integer}
            N = promote_type(T, S)
            ($op)(N(a), N(b))
        end
    end
end

function ~(a::Cenum{T}) where {T<:Integer}
    ~(T(a))
end

Base.convert(::Type{T1}, x::Cenum{T2}) where {T1<:Integer,T2<:Integer} = convert(T1, T2(x))

max_stars_repo_name: MiguelMarcelino/PyEnum

max_issues_repo_name: MiguelMarcelino/PyEnum

max_forks_repo_name: MiguelMarcelino/PyEnum


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using TensorFlow
import TensorFlow.API: mul, softmax, arg_max, cast, reduce_sum,
    expand_dims, tile, l2_normalize, cond, truncated_normal, zeros_like

"""
Soft multiplexer (selector) component that can be learned using gradient
descent
"""
type SoftMux <: AbstractMux
    n_muxinput::Int64
    hidden_units::Vector{Int64}
    muxin::Tensor
    muxselect::Tensor
    softness::Tensor
    nn::ReluStack
    weight::Variable
    nnout::Tensor
    softout::Tensor
    harden::Tensor
    hardselect::Tensor
    hardout::Tensor
    override::Int64
end

using Debug
#"""
#harden::Bool that switches mux output to hard
#"""
@debug function SoftMux(n_muxinput::Int64,
    hidden_units::Vector{Int64},
    muxin::Tensor,
    muxselect::Tensor,
    harden::Tensor,
    softness::Tensor=constant(1.0);
    override::Int64=0)

    relustack = ReluStack(muxselect, hidden_units)
    reluout = out(relustack)

    # last layer is softmax
    # softmax select over inputs
    n0 = get_shape(reluout)[end]
    n1 = n_muxinput
    weight = Variable(truncated_normal(constant([n0, n1]), constant(0.0), constant(5e-2)))
    nnout = softmax(mul(softness, reluout * l2_normalize(weight, constant(1))))

    # mux output is the soft selected input
    if override > 0
        hardselect = tile(constant([override-1]), Tensor([size(nnout)[1]]))
    else
        hardselect = arg_max(nnout, constant(1)) #hardened dense selected channel, 0-indexed
    end
    hardselect_1h = one_hot(hardselect, constant(n_muxinput))
    muxin_rank = ndims(muxin)
    if muxin_rank == 1 || muxin_rank == 2 #TODO: avoid switching if possible
        softout = reduce_sum(mul(muxin, nnout), constant(1))
        hardout = reduce_sum(mul(muxin, hardselect_1h), constant(1))
    elseif muxin_rank == 3
        softout = reduce_sum(mul3(muxin, nnout), constant(2))
        hardout = reduce_sum(mul3(muxin, hardselect_1h), constant(2))
    else
        error("Not supported! (rank=$(muxin_rank))")
    end

    SoftMux(n_muxinput, hidden_units, muxin, muxselect, softness, relustack, weight,
        nnout, softout, harden, hardselect, hardout, override)
end

function mul3(X::Tensor, y::Tensor)
    n_time = get_shape(X)[2]
    tmp = expand_dims(y, constant(1))
    Y = tile(tmp, constant(Int32[1, n_time, 1]))
    out = mul(X, Y) #element-wise mul
    out
end

function out(mux::SoftMux)
    cond_(mux.harden[1], ()->hardout(mux), ()->softout(mux))
end

function softselect(mux::SoftMux)
    mux.nnout
end

function softout(mux::SoftMux)
    mux.softout
end

function hardselect(mux::SoftMux)
    mux.hardselect
end

function hardout(mux::SoftMux)
    mux.hardout
end

function get_variables(mux::SoftMux)
    vars = vcat(mux.weight, get_variables(mux.nn))
    vars
end


text:
using MAT
using PyPlot
using LinearAlgebra
using ADCME
using AdFem
using Distributions

sigma = 0.01
q = (1e-1)^2
m = n = 20
h = 0.1

for sigma in [0.0 0.01]

r = (1e-2)^2

d = matread("simdata/nn2-data-$sigma.mat")
y_, hs, H, s, u_est, kappa = d["yobs"], d["hs"], d["H"], d["s"], d["u_est"], d["kappavalue"]

R = r*diagm(0=>ones(length(hs)))
Q = q*diagm(0=>ones(length(s)))
μ, Σ = uqnlin(y_, hs, H, R, s, Q)
σ = diag(Σ)

umax = maximum(u_est)
xc = LinRange(0, umax, 100)|>Array
u2 = @. 0.1 + 1/(1+xc^2) + 100xc^2
figure()
plot(xc, u2, "-", linewidth=2, label="Referencee")

# k = s[1:4:end]
u1 = fem_to_fvm(u_est, m, n, h)
# Is = sortperm(u1)
# plot(u1[Is], k[Is], "+--")

z1 = μ[1:4:end]
σ1 = 2σ[1:4:end]
u1 = u1[Is]
z1 = z1[Is]
σ1 = σ1[Is]
plot(u1, z1, "--", linewidth=2, label="Posterior Mean")
fill_between(u1, z1-σ1, z1+σ1, alpha=0.3, color="orange", label="Uncertainty Region")

xlim(minimum(u1), maximum(u1))

xlabel("\$u\$")
ylabel("\$K(u)\$")
legend()
savefig("nn2-uq$sigma-1.png")

close("all")
plot(u1, σ1, "+")
xlabel("\$u\$")
ylabel("Standard Deviation")
savefig("nn2-uq$sigma-2.png")

end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using Gadfly

plot(layer(x=[1, 3], y=[1, 3], Geom.line),
     layer(x=[1, 2, 3], y=[1, 2, 3], color=[:a, :b, :c], Geom.point))


text:
#PBC for the linked cells/lists

function wrap_cell(n,icell,jcell)
    if icell < 1
      icell = n
    end
      if icell > n
        icell = 1
      end
    if jcell < 1
      jcell = n
    end
      if jcell > n
        jcell = 1
      end
  return(icell,jcell)
  end


text:
is_constraint(::typeof(voltage_magnitude_constraints)) = true

is_linear(polar::PolarForm, ::typeof(voltage_magnitude_constraints)) = true

# We add constraint only on vmag_pq
function voltage_magnitude_constraints(polar::PolarForm, cons, vm, va, pinj, qinj)
    index_pq = polar.indexing.index_pq
    cons .= @view vm[index_pq]
    return
end
function voltage_magnitude_constraints(polar::PolarForm, cons, buffer)
    index_pq = polar.indexing.index_pq
    cons .= @view buffer.vmag[index_pq]
    return
end

function size_constraint(polar::PolarForm, ::typeof(voltage_magnitude_constraints))
    return PS.get(polar.network, PS.NumberOfPQBuses())
end

function bounds(polar::PolarForm, ::typeof(voltage_magnitude_constraints))
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    fr_ = npq + npv + 1
    to_ = 2*npq + npv
    return polar.x_min[fr_:to_], polar.x_max[fr_:to_]
end

function adjoint!(
    polar::PolarForm,
    pbm::AutoDiff.TapeMemory{F, S, I},
    cons, ∂cons,
    vm, ∂vm,
    va, ∂va,
    pinj, ∂pinj,
) where {F<:typeof(voltage_magnitude_constraints), S, I}
    index_pq = polar.indexing.index_pq
    ∂vm[index_pq] .= ∂cons
end

function matpower_jacobian(
    polar::PolarForm,
    X::Union{State,Control},
    cons::typeof(voltage_magnitude_constraints),
    V,
)
    m = size_constraint(polar, cons)
    nᵤ = get(polar, NumberOfControl())
    nₓ = get(polar, NumberOfState())
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    shift = npq + npv

    I = 1:m
    J = (shift+1):(shift+npq)
    V = ones(m)
    if isa(X, State)
        return sparse(I, J, V, m, nₓ)
    elseif isa(X, Control)
        return spzeros(m, nᵤ)
    end
end

function matpower_hessian(polar::PolarForm, ::typeof(voltage_magnitude_constraints), buffer, λ)
    nu = get(polar, NumberOfControl())
    nx = get(polar, NumberOfState())
    return FullSpaceHessian(
        spzeros(nx, nx),
        spzeros(nu, nx),
        spzeros(nu, nu),
    )
end


text:
using Distributed, StatsBase, Distributions, PyPlot

addprocs(4; exeflags="--project")
@everywhere begin

    using TransitionalMCMC

    # Prior Bounds
    lb, ub  = -5, 5

    # Prior Density and sampler
    priorDen(x) = pdf(Uniform(lb, ub), x[1,:]) .* pdf(Uniform(lb, ub), x[2,:])
    priorRnd(Nsamples) = rand(Uniform(lb, ub), Nsamples, 2)

    # Log Likelihood
    function logLik(x)
        sleep(0.005)
        return -1 .* ((x[1,:].^2 .+ x[2,:] .- 11).^2 .+ (x[1,:] .+ x[2,:].^2 .- 7).^2)
    end

end

Nsamples = 200

samps, acc = tmcmc(logLik, priorDen, priorRnd, Nsamples, 5, 2)

# plt.scatter(samps[:,1], samps[:,2])

rmprocs(workers())


text:
module Tables
   const JuliaConnectoR_DummyTables = true
   istable(x) = false
   columnaccess(x) = false
   columns(x) = nothing
   getcolumn(x) = nothing
   columnnames(x) = nothing
   abstract type AbstractColumns end
end


text:
using Plots
using DifferentialEquations
using DiffEqBiological

sirmodel = @reaction_network SIR begin
    β, S + I --> 2I
    γ, I --> R
end β γ

function solveSIR(u0, β, γ, time_span)

    p = [β, γ]

    dprob = DiscreteProblem(sirmodel, u0, time_span, p)
    jprob = JumpProblem(dprob, Direct(), sirmodel)

    sol = solve(jprob, SSAStepper())

    return sol

end

R0=2.0
N=10000
γ=0.1
β=γ*R0/N

sol = solveSIR([N-1,1,0], β, γ, (0, 160.0))

plot_graph = plot(sol, vars=(0,1), title="SIR model/Gillespie", label="Suscepible",
                  color=:blue, seriestype=:line, legend=:left, size=(1200,900))
plot!(plot_graph, sol, vars=(0,2), label="Infectious", color=:red)
plot!(plot_graph, sol, vars=(0,3), label="Removed", color=:green)

png(plot_graph, "sir_gillespie.png")

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
module FormInputs

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template, register_normal_element

export textfield, numberfield, textarea, filefield

register_normal_element("q__input", context = @__MODULE__)
register_normal_element("q__file", context = @__MODULE__)

function textfield( label::String = "",
                    fieldname::Union{Symbol,Nothing} = nothing,
                    args...;
                    content::Union{String,Vector,Function} = "",
                    kwargs...)
  q__input([isa(content, Function) ? content() : join(content)],
            args...;
            kw([:label => label, :fieldname => fieldname, kwargs...])...)
end

textfield(content::Union{Vector,Function},
          fieldname::Union{Symbol,Nothing} = nothing,
          args...;
          label::String = "",
          kwargs...) = textfield(label, fieldname, args...; content, kwargs...)

function numberfield( label::String = "",
                      fieldname::Union{Symbol,Nothing} = nothing,
                      args...;
                      content::Union{String,Vector,Function} = "",
                      kwargs...)
  q__input( [isa(content, Function) ? content() : join(content)],
            args...;
            kw(
              [:label => label, :fieldname => fieldname, kwargs...],
              Dict("fieldname" => "v-model.number"))...
            )
end

function textarea(label::String = "",
                  fieldname::Union{Symbol,Nothing} = nothing,
                  args...;
                  content::Union{String,Vector,Function} = "",
                  kwargs...)
  textfield(label, fieldname, args...; type="textarea", content=content, kwargs...)
end

function filefield( label::String = "",
                    fieldname::Union{Symbol,Nothing} = nothing,
                    args...;
                    kwargs...)

  q__file(args...;
          kw(
            [:label => label, :fieldname => fieldname, kwargs...],
            Dict("maxfiles" => "max-files", "maxfilesize" => "max-file-size",
                 "counterlabel" => "counter-label", "maxtotalsize" => "max-total-size"
            )
          )...)
end

end


text:
# Autogenerated wrapper script for bitwarden_cli_jll for x86_64-linux-gnu-cxx11
export bw

JLLWrappers.@generate_wrapper_header("bitwarden_cli")
JLLWrappers.@declare_executable_product(bw)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        bw,
        "bin/bw",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@propagate_inbounds function Base.to_indices(A, inds::Tuple{<:AbstractIndex, Vararg{Any}}, I::Tuple{Any, Vararg{Any}})
    Base.@_inline_meta
    (to_index(first(inds), first(I)), to_indices(A, maybetail(inds), tail(I))...)
end

@propagate_inbounds function Base.to_indices(A, inds::Tuple{<:AbstractIndex, Vararg{Any}}, I::Tuple{Colon, Vararg{Any}})
    Base.@_inline_meta
    (values(first(inds)), to_indices(A, maybetail(inds), tail(I))...)
end

@propagate_inbounds function Base.to_indices(A, inds::Tuple{<:AbstractIndex, Vararg{Any}}, I::Tuple{CartesianIndex{1}, Vararg{Any}})
    Base.@_inline_meta
    (to_index(first(inds), first(I)), to_indices(A, maybetail(inds), tail(I))...)
end

maybetail(::Tuple{}) = ()
maybetail(t::Tuple) = tail(t)

"""
    indices(x) -> Tuple
"""
function indices(x::AbstractArray)
    if is_static(x)
        return map(i -> SimpleIndex(OneToSRange(i)), size(x))
    elseif is_fixed(x)
        return map(i -> SimpleIndex(OneTo(i)), size(x))
    else
        return map(i -> SimpleIndex(OneToMRange(i)), size(x))
    end
end

function indices(x::AbstractArray{T,N}, dimnames::Tuple{Vararg{Union{Symbol,Nothing},N}}) where {T,N}
    if is_static(x)
        return map((d, s) -> SimpleIndex{d}(OneToSRange(s)), dimnames, size(x))
    elseif is_fixed(x)
        return map((d, s) -> SimpleIndex{d}(OneTo(s)), dimnames, size(x))
    else
        return map((d, s) -> SimpleIndex{d}(OneToMRange(s)), dimnames, size(x))
    end
end


text:
@testset "flow over circle" begin

end