text:
# ------------------------------------------------
# basic math operators like dot mul,dot add,etc ..
# ------------------------------------------------
import Base.+
import Base.-
import Base.*
import Base.^

export dotAdd
export dotMul
export matAddVec
export matMulVec

function Base.:+(x::Variable{T}, constant) where T
    # a matrix add a constant element by element
    C = eltype(ᵛ(x))(constant)
    y = Variable{T}(ᵛ(x) .+ C, x.backprop)
    if x.backprop
        function matAddScalarBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, matAddScalarBackward)
    end
    return y
end

function Base.:+(constant, var::Variable{T}) where T
    return var + constant;
end

function Base.:-(x::Variable{T}, constant) where T
    # a matrix minus a constant element by element
    C = eltype(ᵛ(x))(constant)
    y = Variable{T}(ᵛ(x) .- C, x.backprop)
    if x.backprop
        function matMinusScalarBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, matMinusScalarBackward)
    end
    return y
end

function Base.:-(constant, x::Variable{T}) where T
    # a matrix minus a constant element by element
    C = eltype(ᵛ(x))(constant)
    y = Variable{T}(C .- ᵛ(x), x.backprop)
    if x.backprop
        function scalarMinusMatBackward()
            if need2computeδ!(x)
                δ(x) .-= δ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, scalarMinusMatBackward)
    end
    return y
end

function Base.:*(x::Variable{T}, constant) where T
    # a matrix multiplies a constant element by element
    C = eltype(ᵛ(x))(constant)
    y = Variable{T}(ᵛ(x) .* C, x.backprop)
    if x.backprop
        function matMulScalarBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* constant
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, matMulScalarBackward)
    end
    return y
end

function Base.:*(constant, var::Variable{T}) where T
    return var * constant
end

function Base.:^(x::Variable{T}, n::Int) where T
    # 矩阵、列向量与常数按元素做幂指数运算
    n = eltype(ᵛ(x))(n)
    y = Variable{T}(ᵛ(x) .^ n, x.backprop)
    if x.backprop
        function powerBackward()
            if need2computeδ!(x)
                δ(x) .+= n .* ᵛ(y) ./ ᵛ(x) .* δ(y);
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, powerBackward)
    end
    return y
end

function Base.:+(x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    # a matrix add a matrix element by element: z = x + y
   @assert T1 <: T2 || T1 >: T2
   T = T1 <: T2 ? T1 : T2
   @assert (x.shape == y.shape) "2 inputs shall be the same size"
   backprop = (x.backprop || y.backprop)
   z = Variable{T}(ᵛ(x) + ᵛ(y), backprop)
   if backprop
       function add2varBackward()
           if need2computeδ!(x) δ(x) .+= δ(z) end
           if need2computeδ!(y) δ(y) .+= δ(z) end
           ifNotKeepδThenFreeδ!(z);
       end
       push!(graph.backward, add2varBackward)
   end
   return z
end

function Base.:-(x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    # a matrix minus a matrix element by element : z = x - y
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    @assert (x.shape == y.shape) "2 inputs shall be the same size"
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) - ᵛ(y), backprop)
    if backprop
        function minus2varBackward()
            if need2computeδ!(x) δ(x) .+= δ(z) end
            if need2computeδ!(y) δ(y) .-= δ(z) end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, minus2varBackward)
    end
    return z
end

"""
    dotAdd(var1::Variable{T1}, var2::Variable{T2}) where {T1,T2}
a tensor add a tensor element by element
"""
function dotAdd(x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    # a tensor add a tensor element by element: z = x .+ y
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    @assert (x.shape == y.shape) "2 inputs shall be the same size"
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) .+ ᵛ(y), backprop)
    if backprop
        function dotAddBackward()
            if need2computeδ!(x) δ(x) .+= δ(z) end
            if need2computeδ!(y) δ(y) .+= δ(z) end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, dotAddBackward)
    end
    return z
end

"""
    dotMul(var1::Variable{T1}, var2::Variable{T2}) where {T1,T2}
a tensor multiplies a tensor element by element
"""
function dotMul(x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    # a tensor multiplies a tensor element by element: z = x .* y
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    @assert (x.shape == y.shape) "2 inputs shall be the same size"
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) .* ᵛ(y), backprop)
    if backprop
        function dotMulBackward()
            if need2computeδ!(x) δ(x) .+= δ(z) .* ᵛ(y) end
            if need2computeδ!(y) δ(y) .+= δ(z) .* ᵛ(x) end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, dotMulBackward)
    end
    return z
end

function Base.:*(W::Variable{T1}, X::Variable{T2}) where {T1,T2}
    # matrix W multiplies matrix X
    # 矩阵相乘 Y[i,j] = sum(W[i,k]*X[k,j],k=...)
    # W -- 权重矩阵
    # X -- n个输入列向量组成的矩阵
    # Y -- n个输出列向量组成的矩阵
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    backprop = (W.backprop || X.backprop)
    Y = Variable{T}(ᵛ(W) * ᵛ(X), backprop)
    if backprop
        function matMulBackward()
            if need2computeδ!(W) δ(W) .+= δ(Y)  * ᵛ(X)' end
            if need2computeδ!(X) δ(X) .+= ᵛ(W)' * δ(Y)  end
            ifNotKeepδThenFreeδ!(Y);
        end
        push!(graph.backward, matMulBackward)
    end
    return Y
end

"""
    matAddVec(var1::Variable{T1}, var2::Variable{T2}) where {T1,T2}
a matrix tensor `var1` adds a vector tensor `var2`
"""
function matAddVec(M::Variable{T1}, V::Variable{T2}) where {T1,T2}
    # M -- 充当和节点，非学习的参数
    # V -- 偏置列向量，要学习的参数
    # Z = M .+ V
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    @assert (M.shape[1]==V.shape[1] && V.shape[2]==1)
    backprop = (M.backprop || V.backprop)
    Z = Variable{T}(ᵛ(M) .+ ᵛ(V), backprop)
    if backprop
        function matAddVecBackward()
            if need2computeδ!(M)
                δ(M) .+= δ(Z)
            end
            if need2computeδ!(V)
                δ(V) .+= sum(δ(Z), dims=2)
            end
            ifNotKeepδThenFreeδ!(Z);
        end
        push!(graph.backward, matAddVecBackward)
    end
    return Z
end

"""
    matAddVec(var1::Variable{T1}, var2::Variable{T2}) where {T1,T2}
a matrix tensor `var1` multiplies a vector tensor `var2`
"""
function matMulVec(M::Variable{T1}, V::Variable{T2}) where {T1,T2}
    # M -- 一般充当激活节点，非网络需要学习的参数
    # V -- 列向量，循环权重，是网络需要学习的参数
    # Z = M .* V
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    @assert (M.shape[1]==V.shape[1] && V.shape[2]==1)
    backprop = (M.backprop || V.backprop)
    Z = Variable{T}(ᵛ(M) .* ᵛ(V), backprop)
    if backprop
        function matMulVecBackward()
            if need2computeδ!(M) δ(M) .+=     δ(Z) .* ᵛ(V)          end
            if need2computeδ!(V) δ(V) .+= sum(δ(Z) .* ᵛ(M), dims=2) end
            ifNotKeepδThenFreeδ!(Z);
        end
        push!(graph.backward, matMulVecBackward)
    end
    return Z
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using Revise
using ADSeismic
using ADCME
using PyPlot
using DelimitedFiles
# matplotlib.use("Agg")
close("all")
if has_gpu()
  use_gpu()
  gpu = true
else
  gpu = false
end

output_dir = "data/acoustic"
if !ispath(output_dir)
  mkpath(output_dir)
end
figure_dir = "figure/acoustic/marmousi/"
if !ispath(figure_dir)
  mkpath(figure_dir)
end
result_dir = "result/acoustic/marmousi/"
if !ispath(result_dir)
  mkpath(result_dir)
end
# if isfile(joinpath(result_dir, "loss.txt"))
#   rm(joinpath(result_dir, "loss.txt"))
# end

################### Generate synthetic data #####################
reset_default_graph()

ap_sim = load_acoustic_model("models/marmousi2-model-true.mat", IT_DISPLAY=0)
src = load_acoustic_source("models/marmousi2-model-true.mat")
rcv_sim = load_acoustic_receiver("models/marmousi2-model-true.mat")
## For debug, only run the first source
# src = [src[1]]
# rcv_sim = [rcv_sim[1]]

if gpu
  Rs_ = compute_forward_GPU(ap_sim, src, rcv_sim)
else
  [SimulatedObservation!(ap_sim(src[i]), rcv_sim[i]) for i = 1:length(src)]
  Rs_ = [rcv_sim[i].rcvv for i = 1:length(rcv_sim)]
end

sess = Session(); init(sess)

Rs = run(sess, Rs_)

for i = 1:length(src)
    writedlm(joinpath(output_dir, "marmousi-r$i.txt"), Rs[i])
end

## visualize_wavefield if needed
# u = run(sess, ap_sim(src[div(length(src),2)+1]).u)
# visualize_wavefield(u, ap_sim.param)

################### Inversion using Automatic Differentiation #####################
reset_default_graph()

ap_sim = load_acoustic_model("models/marmousi2-model-smooth.mat"; inv_vp=true, IT_DISPLAY=0)
src = load_acoustic_source("models/marmousi2-model-true.mat")
rcv_sim = load_acoustic_receiver("models/marmousi2-model-true.mat")
## For debug, only run the first source
# src = [src[1]]
# rcv_sim = [rcv_sim[1]]

Rs = Array{Array{Float64,2}}(undef, length(src))
for i = 1:length(src)
    Rs[i] = readdlm(joinpath(output_dir, "marmousi-r$i.txt"))
end

if gpu
  losses, gs = compute_loss_and_grads_GPU(ap_sim, src, rcv_sim, Rs, get_collection()[1])
  g = sum(gs)
  loss = sum(losses)
else
  [SimulatedObservation!(ap_sim(src[i]), rcv_sim[i]) for i = 1:length(src)]
  loss = sum([sum((rcv_sim[i].rcvv-Rs[i])^2) for i = 1:length(rcv_sim)])
end

sess = Session(); init(sess)

vp = (run(sess, ap_sim.vp.contents)).^0.5./1e3
scale = mean(vp)
vmax = maximum(vp)
vmin = minimum(vp)

function callback(vs, iter, loss)
  if iter%10==0
    if ndims(vs) == 1
      vp = Array(reshape(vs[1:(ap_sim.param.NY+2)*(ap_sim.param.NX+2)], ap_sim.param.NY+2, ap_sim.param.NX+2))
    else
      vp = vs[1]'
    end
    clf()
    pcolormesh([0:ap_sim.param.NX+1;]*ap_sim.param.DELTAX/1e3,[0:ap_sim.param.NY+1;]*ap_sim.param.DELTAY/1e3,  vp*scale)
    axis("scaled")
    colorbar(shrink=0.4)
    xlabel("x (km)")
    ylabel("z (km)")
    gca().invert_yaxis()
    title("Iteration = $iter")
    savefig(joinpath(figure_dir, "inv_$(lpad(iter,5,"0")).png"), bbox_inches="tight")
    writedlm(joinpath(result_dir, "inv_$(lpad(iter,5,"0")).txt"), vp*scale)

    open(joinpath(result_dir, "loss.txt"), "a") do io
      writedlm(io, loss)
    end
  else
  end
end

if gpu
  LBFGS!(sess, loss, g, get_collection()[1]; callback=callback)
else
  LBFGS!(sess, loss, vars=[get_collection()[1]], callback=callback)
end

ext: jl


text:
#=
Model: Building (48 variables, 1 input)
=#
using Reachability, MAT, SparseArrays

building(o::Pair{Symbol, <:Any}...) = building(Options(Dict{Symbol, Any}(o)))

function building(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    file = matopen(@relpath "building.mat")
    A = read(file, "A")

    # initial set
    # - x1-x10 are in [0.0002, 0.00025],
    # - x25 is in [-0.0001, 0.0001],
    # - the rest is 0
    X0 = Hyperrectangle([fill(0.000225, 10); zeros(38)],
                        [fill(0.000025, 10); zeros(14); 0.0001; zeros(23)])

    # input set
    B = read(file, "B")
    U = BallInf([0.9], .1)

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # property: x25 < 6e-3
    property = SafeStatesProperty(HalfSpace(sparsevec([25], [1.0], 48), 6e-3))

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:24] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [25],
                                  :partition => partition,
                                  :plot_vars => [0, 25])
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [25],
                                  :partition => partition,
                                  :property => property)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
@testset "Levels & States" begin
    @testset "J values" begin
        @test J_range(T"¹S") == 0:0
        @test J_range(T"²S") == 1/2:1/2
        @test J_range(T"³P") == 0:2
        @test J_range(T"²D") == 3/2:5/2

        @test J_range(half(1)) == 1/2:1/2
        @test J_range(1) == 1:1
    end

    @testset "Construction" begin
        csfs_3s_3p = csfs(c"3s 3p")
        csf_1 = first(csfs_3s_3p)
        @test_throws ArgumentError Level(csf_1, HalfInteger(2))
        @test_throws ArgumentError Level(csf_1, 2)
        @test_throws ArgumentError State(Level(csf_1, HalfInteger(1)), HalfInteger(2))

        @test string(Level(csf_1, HalfInteger(1))) == "|3s(₁²S|²S) 3p(₁²Pᵒ|¹Pᵒ)-, J = 1⟩"
        @test string(Level(csf_1, 1)) == "|3s(₁²S|²S) 3p(₁²Pᵒ|¹Pᵒ)-, J = 1⟩"
        @test string(State(Level(csf_1, 1), HalfInteger(-1))) == "|3s(₁²S|²S) 3p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = -1⟩"
        @test string(State(Level(csf_1, 1), -1)) == "|3s(₁²S|²S) 3p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = -1⟩"

        @test sort([State(Level(csf_1, HalfInteger(1)), M_J) for M_J ∈ reverse(HalfInteger(-1):HalfInteger(1))]) ==
            [State(Level(csf_1, HalfInteger(1)), M_J) for M_J ∈ HalfInteger(-1):HalfInteger(1)]

        @test states.(csfs_3s_3p) == [[[State(Level(csf_1, HalfInteger(1)), M_J) for M_J ∈ HalfInteger(-1):HalfInteger(1)]],
                                      [[State(Level(csfs_3s_3p[2], J), M_J) for M_J ∈ -J:J] for J ∈ HalfInteger(0):HalfInteger(2) ]]

        rcsfs_3s_3p = csfs(rc"3s" ⊗ rcs"3p")

        @test states.(rcsfs_3s_3p)== [[[State(Level(rcsfs_3s_3p[1], HalfInteger(0)), HalfInteger(0))]],
                                      [[State(Level(rcsfs_3s_3p[2], HalfInteger(1)), M_J) for M_J ∈ HalfInteger(-1):HalfInteger(1)]],
                                      [[State(Level(rcsfs_3s_3p[3], HalfInteger(1)), M_J) for M_J ∈ HalfInteger(-1):HalfInteger(1)]],
                                      [[State(Level(rcsfs_3s_3p[4], HalfInteger(2)), M_J) for M_J ∈ HalfInteger(-2):HalfInteger(2)]]]
    end
end


text:
#= ForwardDiff is removed from dependencies in this version.
import MechanicalSketch
import MechanicalSketch: empty_figure, finish
import MechanicalSketch: PALETTE, settext, setline, WI, EM, HE, O
import MechanicalSketch: dimension_aligned, noise, rk4_steps!
import MechanicalSketch: generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: @import_expand, Quantity, @layer, sethue
import MechanicalSketch: place_image, draw_legend, set_scale_sketch, BinLegend
import MechanicalSketch: normalize_datarange, x_y_iterators_at_pixels
import MechanicalSketch: ∙, ∇_rectangle, SVector, draw_streamlines
import Interpolations: interpolate, Gridded, Linear, Flat, extrapolate
import MechanicalSketch: Greys_9
let
empty_figure(filename = joinpath(@__DIR__, "test_26.png"));
include("test_functions_24.jl")
include("test_functions_26.jl")
# Reuse the flow field from test_23.jl, a matrix of complex velocities with an element per pixel.
A = flowfield_23();
set_scale_sketch(PHYSWIDTH_23, round(Int, SCREEN_WIDTH_FRAC_23 * WI))

str = "The grain size or spectrum is important to Line Integral Convolution (LIC) visualization"
settext(str, O + (-WI/2 + EM, -0.5HE + 2EM), markup = true)

# We are going to use the velocity vector field in a lot of calculations,
# and interpolate between the calculated pixel values
xs, ys = x_y_iterators_at_pixels(;physwidth = PHYSWIDTH_23, physheight = PHYSHEIGHT_23)
fxy_inter = interpolate((xs, reverse(ys)), map( cmplx -> (real(cmplx), imag(cmplx)), transpose(A)[ : , end:-1:1]), Gridded(Linear()));
fxy = extrapolate(fxy_inter, Flat());

# Duration to trace the streamlines forward and back
global const DU_26 = 2.0s

# Upper figure centre coordinate
global const OU_26 = O + (0.0, -0.25HE + EM)

NO_2 = noisepic_2(A, PHYSWIDTH_23, PHYSHEIGHT_23)

# Show the noise image at the top figure
legend = BinLegend(;maxlegend = 1.0, noofbins = 128,
                       colorscheme = reverse(Greys_9), name = Symbol("Value{Float64}"))
upleftpoint, lowrightpoint = place_image(OU_26, legend.(NO_2))
draw_legend(lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23), legend)

# And make a function that interpolates noise between pixels:
nxy_inter = interpolate((xs, reverse(ys)), transpose(NO_2)[ : , end:-1:1], Gridded(Linear()));
nxy = extrapolate(nxy_inter, Flat());
nxy(0.1m, 0.02m)

# Centre of the bottom figure
global const OB_26 = O + (0.0, + 0.25HE + 0.5EM )

# Number of sampling points along our streamline in each direction (effectively one less than this)
global const NS_2 = 10

# Step length (duration)
h = DU_26 / (NS_2 -1 )

M = convolute_image_2(xs, reverse(ys), fxy, nxy, h, CUTOFF_23) # 36.019 s (10 allocations: 4.85 MiB)
                                                  # 27.331 s (50803210 allocations: 2.64 GiB)
                                                  # 29.189352 seconds (30.08 M allocations: 642.345 MiB, 0.39% gc time)
                                                  # 30.222620 seconds (24.06 M allocations: 389.484 MiB, 0.19% gc time)
                                                  # 29.965055 seconds (561.78 k allocations: 31.034 MiB, 0.05% gc time)

upleftpoint, lowrightpoint = place_image(OB_26, legend.(M) )
draw_legend(lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23), legend)

@layer begin
    sethue(PALETTE[1])
    setline(8)
    draw_streamlines(OB_26, fxy; h)
    setline(1)
end

@layer begin
    sethue(PALETTE[1])
    dimension_aligned(OB_26 + (-PHYSWIDTH_23  / 2, PHYSHEIGHT_23  / 2), OB_26 + (PHYSWIDTH_23  / 2, PHYSHEIGHT_23  / 2))
    dimension_aligned(OB_26, OB_26 + complex(0.0, 1.0)m)
    dimension_aligned(OB_26 + (-PHYSWIDTH_23  / 2, - PHYSHEIGHT_23  / 2 ),  OB_26 +  (-PHYSWIDTH_23  / 2, PHYSHEIGHT_23  / 2 ))
end
str = "The longest noise wavelength above is $CUTOFF_23 · $DU_26 = $(CUTOFF_23 * DU_26)"
settext(str, O + (-WI/2 + EM, EM) , markup = true)
str = "Shorter wavelength amplitudes vary. We still need to improve:"
settext(str, O + (-WI/2 + EM, 2 * EM) , markup = true)

finish()
set_scale_sketch()
end
=#


text:
# Tools for implementing device functionality

# julia.h: jl_datatype_align
Base.@pure function datatype_align(::Type{T}) where {T}
    # typedef struct {
    #     uint32_t nfields;
    #     uint32_t alignment : 9;
    #     uint32_t haspadding : 1;
    #     uint32_t npointers : 20;
    #     uint32_t fielddesc_type : 2;
    # } jl_datatype_layout_t;
    field = T.layout + sizeof(UInt32)
    unsafe_load(convert(Ptr{UInt16}, field)) & convert(Int16, 2^9-1)
end

# generalization of word-based primitives

## extract bits from a larger value
@inline function extract_word(val, ::Val{i}) where {i}
    extract_value(val, UInt32, Val(32*(i-1)))
end
@generated function extract_value(val, ::Type{sub}, ::Val{offset}) where {sub, offset}
    T_val = convert(LLVMType, val)
    T_sub = convert(LLVMType, sub)

    bytes = Core.sizeof(val)
    T_int = LLVM.IntType(8*bytes, JuliaContext())

    # create function
    llvm_f, _ = create_function(T_sub, [T_val])
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        equiv = bitcast!(builder, parameters(llvm_f)[1], T_int)
        shifted = lshr!(builder, equiv, LLVM.ConstantInt(T_int, offset))
        # extracted = and!(builder, shifted, 2^32-1)
        extracted = trunc!(builder, shifted, T_sub)

        ret!(builder, extracted)
    end

    call_function(llvm_f, UInt32, Tuple{val}, :( (val,) ))
end

## insert bits into a larger value
@inline function insert_word(val, word::UInt32, ::Val{i}) where {i}
    insert_value(val, word, Val(32*(i-1)))
end
@generated function insert_value(val, sub, ::Val{offset}) where {offset}
    T_val = convert(LLVMType, val)
    T_sub = convert(LLVMType, sub)

    bytes = Core.sizeof(val)
    T_out_int = LLVM.IntType(8*bytes, JuliaContext())

    # create function
    llvm_f, _ = create_function(T_val, [T_val, T_sub])
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        equiv = bitcast!(builder, parameters(llvm_f)[1], T_out_int)
        ext = zext!(builder, parameters(llvm_f)[2], T_out_int)
        shifted = shl!(builder, ext, LLVM.ConstantInt(T_out_int, offset))
        inserted = or!(builder, equiv, shifted)
        orig = bitcast!(builder, inserted, T_val)

        ret!(builder, orig)
    end

    call_function(llvm_f, val, Tuple{val, sub}, :( (val, sub) ))
end

# split the invocation of a function `op` on a value `val` with non-struct eltype
# into multiple smaller invocations on byte-sized partial values.
@generated function split_value_invocation(op::Function, val, args...)
    # TODO: control of lower-limit

    ex = quote
        Base.@_inline_meta
    end

    # disassemble into words
    words = Symbol[]
    for i in 1:Core.sizeof(val)÷4
        word = Symbol("word$i")
        push!(ex.args, :( $word = extract_word(val, Val($i)) ))
        push!(words, word)
    end

    # perform the operation
    for word in words
        push!(ex.args, :( $word = op($word, args...)) )
    end

    # reassemble
    push!(ex.args, :( out = zero(val) ))
    for (i,word) in enumerate(words)
        push!(ex.args, :( out = insert_word(out, $word, Val($i)) ))
    end

    push!(ex.args, :( out ))
    return ex
end

# split the invocation of a function `op` on a value `val`
# by invoking the function on each of its fields
@generated function recurse_value_invocation(op::Function, val, args...)
    ex = quote
        Base.@_inline_meta
    end

    fields = fieldnames(val)
    if isempty(fields)
        push!(ex.args, :( split_value_invocation(op, val, args...) ))
    else
        ctor = Expr(:new, val)
        for field in fields
            push!(ctor.args, :(
                recurse_value_invocation(op, getfield(val, $(QuoteNode(field))), args...) ))
        end
        push!(ex.args, ctor)
    end

    return ex
end

# split the invocation of a function `op` on a pointer `ptr` with non-struct eltype
# into multiple smaller invocations on any supported pointer as listed in `supported_ptrs`.
@generated function split_pointer_invocation(op::Function, ptr, ::Type{supported_ptrs},
                                             args...) where {supported_ptrs}
    T = eltype(ptr)
    elsize(x) = Core.sizeof(eltype(x))
    supported_ptrs = reverse(Base.uniontypes(supported_ptrs))

    ex = quote
        Base.@_inline_meta
    end

    # disassemble
    vals = Tuple{Symbol,Int,Type}[]
    offset = 0
    while offset < Core.sizeof(T)
        val = Symbol("value.$(length(vals)+1)")

        # greedy selection of next pointer type
        remaining = Core.sizeof(T)-offset
        valid = filter(ptr->elsize(ptr)<=remaining, supported_ptrs)
        if isempty(valid)
            error("Cannot partition $T into values of $supported_typs")
        end
        ptr = first(sort(collect(valid); by=elsize, rev=true))

        push!(vals, (val, offset, ptr))
        offset += elsize(ptr)
    end

    # perform the operation
    for (val, offset, ptr) in vals
        subptr = :(convert($ptr, ptr+$offset))
        push!(ex.args, :( $val = op($subptr, args...)) )
    end

    # reassemble
    push!(ex.args, :( out = zero($T) ))
    for (val, offset, ptr) in vals
        push!(ex.args, :( out = insert_value(out, $val, Val($offset)) ))
    end

    push!(ex.args, :( out ))
    return ex
end

# split the invocation of a function `op` on a pointer `ptr`
# by invoking the function on a pointer to each of its fields
@generated function recurse_pointer_invocation(op::Function, ptr, ::Type{supported_ptrs},
                                               args...) where {supported_ptrs}
    T = eltype(ptr)

    ex = quote
        Base.@_inline_meta
    end

    fields = fieldnames(T)
    if isempty(fields)
        push!(ex.args, :( split_pointer_invocation(op, ptr, supported_ptrs, args...) ))
    else
        ctor = Expr(:new, T)
        for (i,field) in enumerate(fields)
            field_typ = fieldtype(T, i)
            field_offset = fieldoffset(T, i)
            field_ptr_typ = :($(ptr.name.wrapper){$field_typ})
            # NOTE: this ctor is a leap of faith
            subptr = :(convert($field_ptr_typ, ptr+$field_offset))
            push!(ctor.args, :(
                recurse_pointer_invocation(op, $subptr, supported_ptrs, args...) ))
        end
        push!(ex.args, ctor)
    end

    return ex
end


text:
import Term: Measure
import Term.measure: width, height
import Term: remove_markup, Panel

@testset "\e[34mMeasure - str" begin
    @test Measure("a"^10).w == 10

    strings = (
        "asadasda"^2,
        "~~±±||__ASdaSDvxcvxc\nasdfsiaudfhsdiufhndskjv",
        "dasdadsa\nsdfsd"^7,
        "asdasdsa|dfvcxvashfusdn\nfidsuhfsdf",
    )
    for string in strings
        m = Measure(string)
        @test m.w == lw(string)
        @test m.h == nlines(string)
    end

    strings = (
        "[red]is my color"^2,
        "~~±±||__[red on_green]aasdas[bold]asdsad[/bold]asdas[/red on_green]xc\nasdfsiaudfhsdiufhndskjv",
        "dasda[#ffffff]dsa[#ffffff]\nsdfsd"^7,
        "asdasdsa|dfvcxvashf[bold]usdn\nfid[/bold]suhfsdf",
    )
    for string in strings
        m = Measure(string)
        @test m.w == lw(remove_markup(string))
        @test m.h == nlines(remove_markup(string))
    end
end

# TODO test with renderables

@testset "\r34mMeasure funcs" begin
    @test width("test") == 4
    @test height("test") == 1

    @test width(Panel(;width=5, height=5)) == 5
    @test height(Panel(;width=5, height=5)) == 5
end


text:
using DecisionTree

D = 10
N = 500
NT = 100

x = randn(N, D)
function rastrigin(x)
  D = length(x)
  10 * D + sum( x.^2 ) - 10 * sum( cos( 2 * π * x ) )
end
y = vec(mapslices(rastrigin, x, 2))
xx = randn(NT, D)
yy = vec(mapslices(rastrigin, xx, 2))

model = build_forest(y, x, 3, 1000, 5, 0.7)
ypred2 = apply_forest(model, xx)

medape(ypred, y) = 100.0 * median(abs((ypred .- y) ./ y))
medape(ypred2, yy)


text:
"""
Process simulation results and arrange a training data for Social GAN
"""
using JDL

path_dir_sim_results = "../sim_records"
files = readdir(path_dir_sim_results)

# TODO:

for f in files
    vars = load(join([path_dir_sim_results, f],"/"))
    for n in vars.rec.nframes
    end
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _airtovac{T<:AbstractFloat}(wave_air::T)
    if wave_air >= 2000.0
        wave_vac = wave_air
        for iter= 1:2
            sigma2 = (1e4/wave_vac)^2.0 # Convert to wavenumber squared
            # Computer conversion factor.
            fact = 1.0 + 5.792105e-2/(238.0185e0 - sigma2) +
                1.67917e-3/(57.362e0 - sigma2)
            wave_vac = wave_air*fact # Convert Wavelength
        end
        return wave_vac
    else
        return wave_air
    end
end

"""
    airtovac(wave_air) -> wave_vacuum

### Purpose ###

Converts air wavelengths to vacuum wavelengths.

### Explanation ###

Wavelengths are corrected for the index of refraction of air under standard
conditions.  Wavelength values below \$2000 Å\$ will not be altered.  Uses
relation of Ciddor (1996).

### Arguments ###

* `wave_air`: can be either a scalar or an array of numbers.  Wavelengths are
corrected for the index of refraction of air under standard conditions.
Wavelength values below \$2000 Å\$ will *not* be altered, take care within \$[1
Å, 2000 Å]\$.

### Output ###

Vacuum wavelength in angstroms, same number of elements as `wave_air`.

### Method ###

Uses relation of Ciddor (1996), Applied Optics 62, 958.

### Example ###

If the air wavelength is `w = 6056.125` (a Krypton line), then `airtovac(w)`
yields a vacuum wavelength of `6057.8019`.

### Notes ###

`vactoair` converts vacuum wavelengths to air wavelengths.

Code of this function is based on IDL Astronomy User's Library.
"""
airtovac(wave_air::Real) = _airtovac(float(wave_air))

@vectorize_1arg Real airtovac


text:
include("KinematicModel.jl")

function vars_state_conservative(m::KinematicModel, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        ρq_tot::FT
    end
end

function vars_state_auxiliary(m::KinematicModel, FT)
    @vars begin
        # defined in init_state_auxiliary
        p::FT
        z::FT
        # defined in update_aux
        u::FT
        w::FT
        q_tot::FT
        q_vap::FT
        q_liq::FT
        q_ice::FT
        e_tot::FT
        e_kin::FT
        e_pot::FT
        e_int::FT
        T::FT
        S::FT
        RH::FT
    end
end

function init_kinematic_eddy!(eddy_model, state, aux, (x, y, z), t)
    FT = eltype(state)

    _grav::FT = grav(param_set)

    dc = eddy_model.data_config

    # density
    q_pt_0 = PhasePartition(dc.qt_0)
    R_m, cp_m, cv_m, γ = gas_constants(param_set, q_pt_0)
    T::FT = dc.θ_0 * (aux.p / dc.p_1000)^(R_m / cp_m)
    ρ::FT = aux.p / R_m / T
    state.ρ = ρ

    # moisture
    state.ρq_tot = ρ * dc.qt_0

    # velocity (derivative of streamfunction)
    ρu::FT =
        dc.wmax * dc.xmax / dc.zmax *
        cos(π * z / dc.zmax) *
        cos(2 * π * x / dc.xmax)
    ρw::FT = 2 * dc.wmax * sin(π * z / dc.zmax) * sin(2 * π * x / dc.xmax)
    state.ρu = SVector(ρu, FT(0), ρw)
    u::FT = ρu / ρ
    w::FT = ρw / ρ

    # energy
    e_kin::FT = 1 // 2 * (u^2 + w^2)
    e_pot::FT = _grav * z
    e_int::FT = internal_energy(param_set, T, q_pt_0)
    e_tot::FT = e_kin + e_pot + e_int
    state.ρe = ρ * e_tot

    return nothing
end

function kinematic_model_nodal_update_auxiliary_state!(
    m::KinematicModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    _grav::FT = grav(param_set)

    aux.u = state.ρu[1] / state.ρ
    aux.w = state.ρu[3] / state.ρ

    aux.q_tot = state.ρq_tot / state.ρ

    aux.e_tot = state.ρe / state.ρ
    aux.e_kin = 1 // 2 * (aux.u^2 + aux.w^2)
    aux.e_pot = _grav * aux.z
    aux.e_int = aux.e_tot - aux.e_kin - aux.e_pot

    # saturation adjustment happens here
    ts = PhaseEquil(param_set, aux.e_int, state.ρ, aux.q_tot)
    pp = PhasePartition(ts)

    aux.T = ts.T
    aux.q_vap = aux.q_tot - pp.liq - pp.ice
    aux.q_liq = pp.liq
    aux.q_ice = pp.ice

    q = PhasePartition(aux.q_tot, aux.q_liq, aux.q_ice)
    aux.S =
        max(
            0,
            aux.q_vap / q_vap_saturation(param_set, aux.T, state.ρ, q) - FT(1),
        ) * FT(100)
    aux.RH =
        aux.q_vap / q_vap_saturation(param_set, aux.T, state.ρ, q) * FT(100)
end

function boundary_state!(
    ::RusanovNumericalFlux,
    m::KinematicModel,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
) end

@inline function wavespeed(
    m::KinematicModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
)
    u = state.ρu / state.ρ
    return abs(dot(nM, u))
end

@inline function flux_first_order!(
    m::KinematicModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)

    # advect moisture ...
    flux.ρq_tot = SVector(
        state.ρu[1] * state.ρq_tot / state.ρ,
        FT(0),
        state.ρu[3] * state.ρq_tot / state.ρ,
    )
    # ... energy ...
    flux.ρe = SVector(
        state.ρu[1] / state.ρ * (state.ρe + aux.p),
        FT(0),
        state.ρu[3] / state.ρ * (state.ρe + aux.p),
    )
    # ... and don't advect momentum (kinematic setup)
end

source!(::KinematicModel, _...) = nothing

function main()
    CLIMA.init()

    # Working precision
    FT = Float64
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δx = FT(20)
    Δy = FT(1)
    Δz = FT(20)
    resolution = (Δx, Δy, Δz)
    # Domain extents
    xmax = 1500
    ymax = 10
    zmax = 1500
    # initial configuration
    wmax = FT(0.6)  # max velocity of the eddy  [m/s]
    θ_0 = FT(289) # init. theta value (const) [K]
    p_0 = FT(101500) # surface pressure [Pa]
    p_1000 = FT(100000) # reference pressure in theta definition [Pa]
    qt_0 = FT(7.5 * 1e-3) # init. total water specific humidity (const) [kg/kg]
    z_0 = FT(0) # surface height

    # time stepping
    t_ini = FT(0)
    t_end = FT(60 * 30)
    dt = 40
    output_freq = 9

    driver_config = config_kinematic_eddy(
        FT,
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        wmax,
        θ_0,
        p_0,
        p_1000,
        qt_0,
        z_0,
    )
    solver_config = CLIMA.SolverConfiguration(
        t_ini,
        t_end,
        driver_config;
        ode_dt = dt,
        init_on_cpu = true,
        #Courant_number = CFL,
    )

    mpicomm = MPI.COMM_WORLD

    # output for paraview
    model = driver_config.bl
    step = [0]
    cbvtk =
        GenericCallbacks.EveryXSimulationSteps(output_freq) do (init = false)
            mkpath("vtk/")
            outprefix = @sprintf(
                "vtk/new_ex_1_mpirank%04d_step%04d",
                MPI.Comm_rank(mpicomm),
                step[1]
            )
            @info "doing VTK output" outprefix
            writevtk(
                outprefix,
                solver_config.Q,
                solver_config.dg,
                flattenednames(vars_state_conservative(model, FT)),
                solver_config.dg.state_auxiliary,
                flattenednames(vars_state_auxiliary(model, FT)),
            )
            step[1] += 1
            nothing
        end

    # get aux variables indices for testing
    q_tot_ind = varsindex(vars_state_auxiliary(model, FT), :q_tot)
    q_vap_ind = varsindex(vars_state_auxiliary(model, FT), :q_vap)
    q_liq_ind = varsindex(vars_state_auxiliary(model, FT), :q_liq)
    q_ice_ind = varsindex(vars_state_auxiliary(model, FT), :q_ice)
    S_ind = varsindex(vars_state_auxiliary(model, FT), :S)

    # call solve! function for time-integrator
    result = CLIMA.invoke!(
        solver_config;
        user_callbacks = (cbvtk,),
        check_euclidean_distance = true,
    )

    # no supersaturation
    max_S = maximum(abs.(solver_config.dg.state_auxiliary[:, S_ind, :]))
    @test isequal(max_S, FT(0))

    # qt is conserved
    max_q_tot = maximum(abs.(solver_config.dg.state_auxiliary[:, q_tot_ind, :]))
    min_q_tot = minimum(abs.(solver_config.dg.state_auxiliary[:, q_tot_ind, :]))
    @test isapprox(max_q_tot, qt_0; rtol = 1e-3)
    @test isapprox(min_q_tot, qt_0; rtol = 1e-3)

    # q_vap + q_liq = q_tot
    max_water_diff = maximum(abs.(
        solver_config.dg.state_auxiliary[:, q_tot_ind, :] .-
        solver_config.dg.state_auxiliary[:, q_vap_ind, :] .-
        solver_config.dg.state_auxiliary[:, q_liq_ind, :],
    ))
    @test isequal(max_water_diff, FT(0))

    # no ice
    max_q_ice = maximum(abs.(solver_config.dg.state_auxiliary[:, q_ice_ind, :]))
    @test isequal(max_q_ice, FT(0))

    # q_liq ∈ reference range
    max_q_liq = max(solver_config.dg.state_auxiliary[:, q_liq_ind, :]...)
    min_q_liq = min(solver_config.dg.state_auxiliary[:, q_liq_ind, :]...)
    @test max_q_liq < FT(1e-3)
    @test isequal(min_q_liq, FT(0))
end

main()


text:
module DiffEqSensitivity

using DiffEqBase, ForwardDiff, Tracker, DiffEqDiffTools, Statistics
using DiffEqCallbacks, QuadGK, RecursiveArrayTools, LinearAlgebra
using DataFrames, GLM, RecursiveArrayTools, QuasiMonteCarlo
using Parameters: @unpack, @with_kw
using FFTW, Distributions
import ZygoteRules, ChainRulesCore

abstract type SensitivityFunction end
abstract type GSAMethod end

include("derivative_wrappers.jl")
include("local_sensitivity/sensitivity_algorithms.jl")
include("local_sensitivity/sensitivity_interface.jl")
include("local_sensitivity/forward_sensitivity.jl")
include("local_sensitivity/backsolve_adjoint.jl")
include("local_sensitivity/interpolating_adjoint.jl")
include("local_sensitivity/quadrature_adjoint.jl")
include("local_sensitivity/concrete_solve.jl")
include("global_sensitivity/morris_sensitivity.jl")
include("global_sensitivity/sobol_sensitivity.jl")
include("global_sensitivity/regression_sensitivity.jl")
include("global_sensitivity/DGSM_sensitivity.jl")
include("global_sensitivity/eFAST_sensitivity.jl")

export extract_local_sensitivities

export ODEForwardSensitivityFunction, ODEForwardSensitivityProblem, SensitivityFunction,
       ODEAdjointSensitivityProblem, ODEAdjointProblem, AdjointSensitivityIntegrand,
       adjoint_sensitivities, adjoint_sensitivities_u0, Sobol, Morris, gsa,
       SensitivityAlg, regression_sensitivity, DGSM, eFAST

export BacksolveAdjoint, QuadratureAdjoint, InterpolatingAdjoint, ForwardSensitivity,
       ForwardDiffSensitivity

end # module


text:
"""
PDE_aa(n=31, m=3; Vstore=Float64[], pdata=nothing, beta=1.0, maxit=40)
Solve preconditioned convection-diffusion equation with hardwired left
preconditioner using Anderson acceleration.

If you're putting this in a loop, you should allocate Vstore to
zeros(n*n,2*(mmax+1)). Otherwise I will make a decision for you and
allocate for Vstore with each call to this function. The story on
pdata is the same. If you are calling this several times with the
same value of n, build pdata outside the call.
"""
function PDE_aa(n=31, m=3; Vstore=Float64[], pdata=nothing, beta=1.0, maxit=40)
#
# Process Vstore and pdata
#
(pdata != nothing) || (pdata=pdeinit(n));
(length(Vstore)>0) || (Vstore=zeros(n*n,3*m + 3))
(mv, nv)  = size(Vstore)
dimvtest= ((mv == n*n) && (nv >= 2*m + 4))
dimvtest || error("Vstore too small")
#
# Call aasol and return the results.
#
u0=zeros(n*n);
rtol = 1.e-8;
atol = 1.e-8;
aout = aasol(hardleftFix!, u0, m, Vstore; pdata=pdata, maxit=maxit,
             rtol=rtol, atol=atol, beta=beta);
return aout
end


text:
using Soss
using LaTeXStrings

m = @model x begin
    μ ~ Normal(0,5)
    σ ~ HalfCauchy(3)
    N = length(x)
    x ~ Normal(μ, σ) |> iid(N)
end

f = makeLogdensity(m)
μμ = range(-10,10,length=1000)

normalize(seq) = seq ./ sum(seq)
plot(μμ, [exp(f((μ=μ,σ=1,x=[]))) for μ in μμ] |> normalize, label=L"P(\mu)",legendfontsize=16,dpi=300)
plot!(μμ, [exp(f((μ=μ,σ=1,x=[-5]))) for μ in μμ] |> normalize, label=L"P(\mu|x=-5)")
savefig("/home/chad/prior-posterior.png")

mPrior = prior(m)

rand(mPrior)

particles(mPrior)

nuts(m; x=randn(10))

m2 = @model x,errDist begin
    μ ~ Normal(0,5)
    σ ~ HalfCauchy()
    N = length(x)
    d = LocationScale(μ,σ,errDist)
    x ~ d |> iid(N)
end

nuts(m; x=randn(10), errDist = Normal()).samples[1:3] |> DataFrame

nuts(m; x = randn(10), errDist = TDist(3)).samples[1:3] |> DataFrame

# tdist = @model ν begin
#     w ~ InverseGamma(ν / 2, ν / 2)
#     x ~ Normal(0, w ^ 2)
#     return x
# end

nuts(m; x = randn(10), errDist = tdist, ν=3).samples[1:3] |> DataFrame

p = @model begin
    α ~ Normal(1,1)
    β ~ Normal(α^2,1)
end

q = @model μα,σα,μβ,σβ begin
    α ~ Normal(μα,σα)
    β ~ Normal(μβ,σβ)
end

f = sourceParticleImportance(p,q) |> eval

f = sourceParticleImportance(p,q) |> eval
(ℓ,θ) = f(1000,(μα=1, σα=1, μβ=2, σβ=2.6))
@unpack (α,β) = θ
scatter(α.particles, β.particles
    , alpha=exp(ℓ - maximum(ℓ)).particles
    , legend=false
    , xlim=(-2,3)
    , ylim=(-3,10)
    , xlabel=L"\alpha"
    , ylabel=L"\beta")
savefig("/home/chad/banana.svg")

q = @model μm,μs,σm,σs begin
    μ ~ Normal(μm,μs)
    σ ~ LogNormal(σm,σs)
end

using LaTeXStrings
f = sourceParticleImportance(p,q) |> eval
(ℓ,θ) = f(1000,(μα=1, σα=1, μβ=2, σβ=2.6))
@unpack (α,β) = θ
scatter(α.particles, β.particles
    , alpha=exp(ℓ - maximum(ℓ)).particles
    , legend=false
    , xlim=(-2,3)
    , ylim=(-3,10)
    , xlabel=L"\alpha"
    , ylabel=L"\beta")
savefig("/home/chad/banana.svg")

#############################
# Example problem


text:
using Documenter, MethodAnalysis

using InteractiveUtils
display(versioninfo())

makedocs(;
    modules=[MethodAnalysis],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/timholy/MethodAnalysis.jl/blob/{commit}{path}#L{line}",
    sitename="MethodAnalysis.jl",
    authors="Tim Holy <tim.holy@gmail.com>",
)

deploydocs(;
    repo="github.com/timholy/MethodAnalysis.jl",
    push_preview=true,
)


text:
using StackPointers
using Test

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


text:
printstyled("  findid\n", color=:light_green)
id = "UW.TDH..EHZ"
IDs = ["UW.WWVB..TIM","UW.TCG..TIM","UW.TDH..EHZ","UW.VLM..EHZ"]
@test findid(id, IDs) == findid(codeunits(id), IDs) == findfirst(IDs.==id) == 3
@test findid("aslkasnglknsgf", IDs) == 0


text:
### Without renewables for validation of IEEE test cases
###################
# Run DCOPF #######
###################

DCOPF_mod, P_opt, Θ_opt, P_inj_opt, P_flow_opt =
	DCOPF(time, nodes, lines, pp)

# Optimisation status
JuMP.termination_status(DCOPF_mod)
JuMP.primal_status(DCOPF_mod)
JuMP.dual_status(DCOPF_mod)

# Results
JuMP.objective_value(DCOPF_mod)
P_opt
Θ_opt
P_inj_opt
P_flow_opt

###################
# Run ED ##########
###################

ED_mod, P_opt, price = EconomicDispatch(time, nodes, pp)

# Optimisation status
JuMP.termination_status(ED_mod)
JuMP.primal_status(ED_mod)
JuMP.dual_status(ED_mod)

# Results
JuMP.objective_value(ED_mod)
P_opt
price

# Dispatch: Unit commitment and power generation
output_EconomicDispatch = DataFrame(time=Int64[],unit=Int64[],output=Float64[])
for t in time
	for g in pp.unit
		val = P_opt[t, g]
		push!(output_EconomicDispatch, [t g val])
	end
end

CSV.write("output/EconomicDispatch.csv", output_EconomicDispatch)
electrondisplay(output_EconomicDispatch)
electrondisplay(powerplants_df)

###################
# Run CM ##########
###################

CM_mod, ΔP_up_opt, ΔP_dn_opt, Θ_opt, P_inj_opt, P_flow_opt =
	CongestionManagement(time, nodes, lines, price, P_opt, pp)

# Optimisation status
JuMP.termination_status(CM_mod)
JuMP.primal_status(CM_mod)
JuMP.dual_status(CM_mod)

# Results
JuMP.objective_value(CM_mod)
ΔP_up_opt
ΔP_dn_opt
Θ_opt
P_inj_opt
P_flow_opt

# Testing
P_opt.data .+ ΔP_up_opt.data .- ΔP_dn_opt.data
P_opt.data


text:
#! /usr/bin/julia

# Rosetta Code, Arithmetic/Rational

function isperfect{T<:Integer}(n::T)
    !isprime(n) || return false
    tal = 1//n
    hi = isqrt(n)
    if hi^2 == n
        tal += 1//hi
        hi -= 1
    end
    for i in 2:hi
        (d, r) = divrem(n, i)
        if r == 0
            tal += (1//i + 1//d)
        end
    end
    return tal == 1//1
end

lo = 2
hi = 2^19
println("Searching for perfect numbers from ", lo, " to ", hi, ".")
for i in 2:2^19
    isperfect(i) || continue
    println(@sprintf("%8d", i))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
function arithmetic(n)
    sum = 0
    for i = 1:n
        sum += i
    end
    return sum
end

function main()
    n::Int64 = 12
    arithmetic(n)
end

# println("arithmetic, Int64:")
# code_native(arithmetic, (Int64,), syntax = :intel)

# println("arithmetic, Float64:")
# code_native(arithmetic, (Float64,), syntax = :intel)

# println("main:")
# code_native(main, (), syntax = :intel)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


text:
@testset "show" begin
    x = CSTParser.parse("a + (b*c) - d")
    @test sprint(show, x) ===
    """
      1:13  BinaryOpCall
      1:10   BinaryOpCall
      1:2     a
      3:4     OP: PLUS
      5:10    InvisBrackets
      5:5      (
      6:8      BinaryOpCall
      6:6       b
      7:7       OP: STAR
      8:8       c
      9:10     )
     11:12   OP: MINUS
     13:13   d
    """
end


text:
#################################################################################################
# AIM 1: PR
#################################################################################################
## Deps
using Distributed
using CodecZlib
using StatsPlots
using Distributions
using DelimitedFiles
using Plots.PlotMeasures

## Constants
const noise_levels = ["$(s)" for s in collect(2.0:0.5:3.5)]
const data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-1"

# Thresholds
const nano_thrs = vcat(collect(-500.0:50.0:-20.0), collect(-20.0:0.1:20.0), collect(20.0:50.0:500.0))

# Colors
const call_thrs = Dict("Nanopolish" => nano_thrs)
const cllr_color_code = Dict("Nanopolish" => "#D55E00")
const blnd_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]
const s_color_cde = Dict("0.5" => blnd_col[1], "1.0" => blnd_col[2], "1.5" => blnd_col[3], "2.0" => blnd_col[4],
                         "2.5" => blnd_col[5], "3.0" => blnd_col[6], "3.5" => blnd_col[7], "4.0" => blnd_col[8])

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

# Function to threshold calls
@everywhere function thresh_calls(conf_x, thrs)

    # Init output vectors
    out_conf_x = []

    # Compare confidence score to threshold to determine call
    @inbounds for i = 1:length(conf_x)

        if conf_x[i] == "n/a"
            # Fails to detect
            push!(out_conf_x, 0)
        elseif conf_x[i] < thrs
            # Negative
            push!(out_conf_x, -1)
        elseif conf_x[i] >= thrs
            # Positive
            push!(out_conf_x, 1)
        else
            println("Something went wrong when thresholding calls ...")
        end

    end

    return out_conf_x

end

# Function to compute sensitivity (recall, true positive rate)
@everywhere function comp_sens(truth, pred)

    # Get positives
    pos_ind = pred .== 1

    # Get true positives
    true_pos = sum(pos_ind) > 0 ? sum(truth[pos_ind] .== pred[pos_ind]) : 0.0

    # Return sensitivity
    return true_pos / sum(truth .== 1)

end

# Function to compute precision
@everywhere function comp_prec(truth, pred)

    # Get positives
    pos_ind = pred .== 1

    # Get true positives
    true_pos = sum(pos_ind) > 0 ? sum(truth[pos_ind] .== pred[pos_ind]) : 0.0

    # Return precision
    return true_pos / sum(pos_ind)

end

@everywhere function pmap_thresh(true_x, conf_x, t)

    # Threshold calls
    thresh_x = thresh_calls(conf_x, t)

    # Return
    return comp_prec(true_x, thresh_x), comp_sens(true_x, thresh_x)

end

function get_prec_recall(caller, s)

    # Read in data
    stream = GzipDecompressorStream(open("$(data_dir)/$(caller)/gm12878_chr22_sigma_$(s)_$(caller)_tuples_sample.tsv.gz"))
    in_data = readdlm(stream)
    close(stream)

    # Remove NaN
    kp_ind = .! map(isequal("n/a"), in_data[:,2])
    in_data = in_data[kp_ind,:]

    # Get true calls and confidence scores for predictions
    true_x = in_data[:,1]
    conf_x = in_data[:,2]

    # Get tpr and fpr
    pmap_out = pmap(t -> pmap_thresh(true_x, conf_x, t), call_thrs[caller])

    # Disentangle data
    prec_vec = [round(x[1], digits=3) for x in pmap_out]
    recall_vec = [round(x[2], digits=3) for x in pmap_out]

    # Remove NaN
    kp_ind = isfinite.(recall_vec) .& isfinite.(prec_vec)
    recall_vec = recall_vec[kp_ind]
    prec_vec = prec_vec[kp_ind]

    # Clean data
    kp_ind = []
    sort_ind = sortperm(recall_vec)
    recall_vec = recall_vec[sort_ind]
    prec_vec = prec_vec[sort_ind]
    for i = 2:length(prec_vec)
        if prec_vec[i] < prec_vec[i - 1]
            push!(kp_ind, i)
        end
    end
    recall_vec = recall_vec[kp_ind]
    prec_vec = prec_vec[kp_ind]

    # Return tuple
    return recall_vec, prec_vec

end

function get_auc(recall_vec, prec_vec)

    # Integrate ROC curve
    auc = 0.0
    @inbounds for i in 1:(length(prec_vec) - 1)
        x = recall_vec[i + 1] - recall_vec[i]
        y = prec_vec[i]
        auc += x * y
    end

    # Return AUC
    return auc

end

#################################################################################################
# PRC & AUC
#################################################################################################

# Init plots
p = plot(xlabel="sensitivity (recall)", ylabel="precision", title="PRC - Nanopolish", xlim=(0.0, 1.0), ylim=(0.5, 1.0));

# Init AUROC array
auc_nano = []

# Loop over noise levels
for s in noise_levels

    # Print noise level
    println("Working on sigma=$(s)")

    # Get precision & recall
    recall_vec, prec_vec = get_prec_recall("Nanopolish", s)

    # Plot curve
    col = s_color_cde[s]
    plot!(p, recall_vec, prec_vec, seriestype=:line, alpha=0.5, color=col, label="sd=$(s)")

    # Calculate AUC
    push!(auc_nano, get_auc(recall_vec, prec_vec))

end

# Plot ROC ROC_Curves
plot(p, size=(700, 600))
savefig("$(data_dir)/PRC-Curves.pdf")

# Plot AUROC
col = cllr_color_code["Nanopolish"]
p4 = plot(xlabel="signal noise level (sd)", ylabel="AUC-PRC", ylim=(0.5, 1.0));
plot!(p4, noise_levels, auc_nano, seriestype=:scatter, markershape=:circle, color=col, label="")
plot!(p4, noise_levels, auc_nano, seriestype=:line, alpha=0.5, color=col, label="")
plot(p4, size=(700, 600), legend=:bottomleft)
savefig("$(data_dir)/AUC-PRC.pdf")


text:
@testset "83.remove-duplicates-from-sorted-list.jl" begin
    @test delete_duplicates!(ListNode{Int}([1, 1, 2])) == ListNode{Int}([1, 2])
    @test delete_duplicates!(ListNode{Int}([1, 1, 2, 3, 3])) == ListNode{Int}([1, 2, 3])
    @test delete_duplicates!(ListNode{Int}([1, 1, 1])) == ListNode(1)
end


text:
struct GPSL1{C <: AbstractMatrix} <: AbstractGNSS{C}
    codes::C
end

get_system_string(s::GPSL1) = "GPSL1"

function read_gpsl1_codes()
    read_in_codes(
        Int8,
        joinpath(dirname(pathof(GNSSSignals)), "..", "data", "codes_gps_l1.bin"),
        37,
        1023
    )
end

function GPSL1(;use_gpu = Val(false))
    _GPSL1(use_gpu)
end

function _GPSL1(::Val{false})
    GPSL1(Int16.(read_gpsl1_codes()))
end
function _GPSL1(::Val{true})
    GPSL1(CuMatrix{Float32}(read_gpsl1_codes()))
end

"""
$(SIGNATURES)

Get code length of GNSS system GPSL1.
```julia-repl
julia> get_code_length(gpsl1)
```
"""
@inline function get_code_length(gpsl1::GPSL1)
    1023
end

"""
$(SIGNATURES)

Get secondary code length of GNSS system GPSL1.
```julia-repl
julia> get_secondary_code_length(gpsl1)
```
"""
@inline function get_secondary_code_length(gpsl1::GPSL1)
    1
end

"""
$(SIGNATURES)

Get center frequency of GNSS system GPSL1.
```julia-repl
julia> get_center_frequency(gpsl1)
```
"""
@inline function get_center_frequency(gpsl1::GPSL1)
    1_575_420_000Hz
end

"""
$(SIGNATURES)

Get code frequency of GNSS system GPSL1.
```julia-repl
julia> get_code_frequency(gpsl1)
```
"""
@inline function get_code_frequency(gpsl1::GPSL1)
    1_023_000Hz
end
"""
$(SIGNATURES)

Get data frequency of GNSS system GPSL1.
```julia-repl
julia> get_data_frequency(gpsl1)
```
"""
@inline function get_data_frequency(gpsl1::GPSL1)
    50Hz
end

"""
$(SIGNATURES)
Get the spectral power of the GPSL1 CA code
"""
function get_code_spectrum(s::GPSL1, f)
    get_code_spectrum_BPSK(get_code_frequency(s), f)
end


text:
using ACE1pack, JuLIP, LazyArtifacts, Test
using JuLIP.Testing: print_tf

test_train_set = joinpath(artifact"TiAl_tiny_dataset", "TiAl_tiny.xyz")
expected_errors_json = joinpath(artifact"ACE1pack_test_files", "expected_fit_errors.json")

@info("test full fit from script")

species = [:Ti, :Al]
r0 = 2.88

data = data_params(fname = test_train_set,
    energy_key = "energy",
    force_key = "force",
    virial_key = "virial")

ace_basis = basis_params(
    type = "ace",
    species = species,
    N = 3,
    maxdeg = 6,
    r0 = r0,
    radial = basis_params(
        type = "radial",
        rcut = 5.0,
        rin = 1.44,
        pin = 2))

pair_basis = basis_params(
    type = "pair",
    species = species,
    maxdeg = 6,
    r0 = r0,
    rcut = 5.0,
    rin = 0.0,
    pcut = 2, # TODO: check if it should be 1 or 2?
    pin = 0)

basis = Dict(
    "ace" => ace_basis,
    "pair" => pair_basis
)

solver = solver_params(type = :lsqr)

# symbols for species (e.g. :Ti) would work as well
e0 = Dict("Ti" => -1586.0195, "Al" => -105.5954)

weights = Dict(
    "default" => Dict("E" => 5.0, "F" => 1.0, "V" => 1.0),
    "FLD_TiAl" => Dict("E" => 5.0, "F" => 1.0, "V" => 1.0),
    "TiAl_T5000" => Dict("E" => 30.0, "F" => 1.0, "V" => 1.0))

P = regularizer_params(type = "laplacian", rlap_scal = 3.0)

params = fit_params(
    data = data,
    basis = basis,
    solver = solver,
    e0 = e0,
    weights = weights,
    P = P,
    ACE_fname = "")

IP, lsqinfo = ACE1pack.fit_ace(params)

errors = lsqinfo["errors"]

expected_errors = load_dict(expected_errors_json)

for error_type in keys(errors),
        config_type in keys(errors[error_type]),
            property in keys(errors[error_type][config_type])
    print_tf(@test errors[error_type][config_type][property] <= 2 * expected_errors[error_type][config_type][property])
end
println()


text:
module TypedMissings

export TypedMissing, MissingKind, MissingKinds
module MissingKinds
    export MissingKind

    @enum MissingKind::UInt8 begin
        NI=0    # No information
        INV=1   # Invalid
        OTH=2   # Other
        NINF=3  # Negative infinity
        PINF=4  # Positive infinity
        UNC=5   # Unencoded
        DER=6   # Derived
        UNK=7   # Unknown
        ASKU=8  # Asked but unknown
        NAV=9   # Temporarily not available
        NAVU=10 # Not available
        NASK=11 # Not asked
        QS=12   # Sufficient quantity
        TRC=13  # Trace
        MSK=14  # Masked
        NA=15   # Not applicable

        a=101
        b=102
        c=103
        d=104
        e=105
        f=106
        g=107
        h=108
        i=109
        j=110
        k=111
        l=112
        m=113
        n=114
        o=115
        p=116
        q=117
        r=118
        s=119
        t=120
        u=121
        v=122
        w=123
        x=124
        y=125
        z=126
    end
end
using .MissingKinds

"""
    TypedMissing(kind::MissingKind=MissingKinds.NI)

A type similar to `Missing` that allows representing multiple kinds
of missing values. The default kind is `MissingKinds.NI` ("No information"),
which is equivalent to `missing`.
Values of type `TypedMissing` behave identically to `missing`, except that
`isequal(TypedMissing(kind), missing)` is `true` if and only if
`kind == MissingKinds.NI`.

`TypedMissing` values propagate their `kind` across operations
which involve only missing values of the same `kind` or non-missing
values. Operations mixing `TypedMissing` values of different `kind`s
fall back to the `MissingKinds.NI` kind.

If provided, `kind` must be an instance of the `MissingKind` enum.
Supported `kind`s are `NI`, lowercase letters from a to z, and "null flavour"
values [defined by FHIR HL7 v3](https://www.hl7.org/fhir/v3/NullFlavor/cs.html)/
ISO 21090:
- NI: No information
- INV: Invalid
- OTH: Other
- NINF: Negative infinity
- PINF: Positive infinity
- UNC: Unencoded
- DER: Derived
- UNK: Unknown
- ASKU: Asked but unknown
- NAV: Temporarily not available
- NAVU: Not available
- QS: Sufficient quantity
- NASK:Not asked
- TRC: Trace
- MSK: Masked
- NA: Not applicable

# Examples
```jldoctest
julia> TypedMissing()
TypedMissing()

julia> TypedMissing(MissingKinds.NI)
TypedMissing()

julia> TypedMissing() + 1
TypedMissing()

julia> TypedMissing(MissingKinds.a) + 1
TypedMissing(MissingKinds.a)

julia> TypedMissing(MissingKinds.NASK) + 1
TypedMissing(MissingKinds.NASK)

julia> TypedMissing(MissingKinds.NASK) + TypedMissing(MissingKinds.INV)
TypedMissing()

julia> TypedMissing(MissingKinds.NASK) + missing
TypedMissing()

julia> TypedMissing(MissingKinds.NASK) + 1
TypedMissing(MissingKinds.NASK)

julia> isequal(TypedMissing(MissingKinds.NASK), TypedMissing(MissingKinds.NASK))
true

julia> isequal(TypedMissing(MissingKinds.NASK), TypedMissing(MissingKinds.INV))
false

julia> isequal(TypedMissing(MissingKinds.NASK), missing)
false

julia> isequal(TypedMissing(MissingKinds.NI), missing)
true
```
"""
struct TypedMissing
    kind::MissingKind
end
TypedMissing() = TypedMissing(MissingKinds.NI)

function Base.show(io::IO, x::TypedMissing)
    show(io, TypedMissing)
    print(io, '(')
    if x.kind != MissingKinds.NI
        print(io, "MissingKinds.", x.kind)
    end
    print(io, ')')
end

function Base.show(io::IO, ::MIME"text/plain", x::TypedMissing)
    show(io, TypedMissing)
    print(io, '(')
    if x.kind != MissingKinds.NI
        get(io, :compact, false) || print(io, "MissingKinds.")
        print(io, x.kind)
    end
    print(io, ')')
end

Base.ismissing(::TypedMissing) = true
Base.coalesce(x::TypedMissing, y...) = coalesce(y...)

# Semi type piracy
Base.nonmissingtype(::Type{T}) where {T >: TypedMissing} =
    Base.typesplit(T, Union{Missing, TypedMissing})

Base.broadcastable(x::TypedMissing) = Ref(x)

Base.promote_rule(T::Type{TypedMissing}, S::Type) = Union{S, TypedMissing}
Base.promote_rule(T::Type{TypedMissing}, S::Type{Missing}) = TypedMissing
Base.promote_rule(T::Type{Missing}, S::Type{TypedMissing}) = TypedMissing
Base.promote_rule(T::Type{Union{Nothing, TypedMissing}}, S::Type) = Union{S, Nothing, TypedMissing}

function _promote_nothing(T::Type, S::Type)
    R = Base.nonnothingtype(T)
    R >: T && return Any
    T = R
    R = nonmissingtype(T)
    R >: T && return Any
    T = R
    R = promote_type(T, S)
    return Union{R, Nothing, TypedMissing}
end
Base.promote_rule(T::Type{>:Union{Nothing, TypedMissing}}, S::Type) = _promote_nothing(T, S)
Base.promote_rule(T::Type{>:Union{Nothing, Missing, TypedMissing}}, S::Type) =
    _promote_nothing(T, S)
Base.promote_rule(T::Type{>:Union{Nothing, Missing, TypedMissing}}, S::Type{>:TypedMissing}) =
    _promote_nothing(T, S)
Base.promote_rule(T::Type{>:Union{Nothing, Missing}}, S::Type{>:TypedMissing}) =
    _promote_nothing(T, S)

function _promote_missing(T::Type, S::Type)
    R = nonmissingtype(T)
    R >: T && return Any
    T = R
    R = promote_type(T, nonmissingtype(S))
    return Union{R, TypedMissing}
end
Base.promote_rule(T::Type{>:TypedMissing}, S::Type) =
    _promote_missing(T, S)
Base.promote_rule(T::Type{>:Union{Missing, TypedMissing}}, S::Type) =
    _promote_missing(T, S)
Base.promote_rule(T::Type{>:Missing}, S::Type{>:TypedMissing}) =
    _promote_missing(T, S)
Base.promote_rule(T::Type{>:Union{Missing, TypedMissing}}, S::Type{>:TypedMissing}) =
    _promote_missing(T, S)
Base.promote_rule(T::Type{Missing}, S::Type{>:TypedMissing}) =
    _promote_missing(T, S)

Base.convert(::Type{T}, x::T) where {T>:TypedMissing} = x
Base.convert(::Type{T}, x::T) where {T>:Union{TypedMissing, Nothing}} = x
Base.convert(::Type{T}, x::T) where {T>:Union{TypedMissing, Missing}} = x
Base.convert(::Type{T}, x) where {T>:TypedMissing} =
    convert(Base.nonmissingtype_checked(T), x)
Base.convert(::Type{T}, x::Missing) where {T>:TypedMissing} =
    TypedMissing()
Base.convert(::Type{T}, x) where {T>:Union{TypedMissing, Nothing}} =
    convert(Base.nonmissingtype_checked(Base.nonnothingtype_checked(T)), x)
Base.convert(::Type{T}, x) where {T>:Union{TypedMissing, Missing}} =
    convert(Base.nonmissingtype_checked(T), x)
Base.convert(::Type{T}, x) where {T>:Union{TypedMissing, Missing, Nothing}} =
    convert(Base.nonmissingtype_checked(Base.nonnothingtype_checked(T)), x)
Base.convert(::Type{T}, x::Missing) where {T>:Union{TypedMissing, Nothing}} =
    TypedMissing()
Base.convert(::Type{T}, x::Missing) where {T>:Union{TypedMissing, Missing}} =
    missing
Base.convert(::Type{T}, x::Missing) where {T>:Union{TypedMissing, Missing, Nothing}} =
    missing
Base.convert(::Type{Any}, x::Missing) = missing
Base.convert(::Type{T}, x::TypedMissing) where {T>:Union{TypedMissing, Missing, Nothing}} =
    x
Base.convert(::Type{Any}, x::TypedMissing) = x
Base.convert(::Type{T}, x::TypedMissing) where {T>:Missing} = missing
Base.convert(::Type{T}, x::TypedMissing) where {T>:Union{Missing, Nothing}} = missing

Array{T,N}(x::TypedMissing, d...) where {T,N} = fill!(Array{T,N}(undef, d...), x)
Array{T}(x::TypedMissing, d...) where {T} = fill!(Array{T}(undef, d...), x)

# Comparison operators
Base.:(==)(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.:(==)(x::TypedMissing, ::Any) = x
Base.:(==)(::Any, y::TypedMissing) = y
# To fix ambiguity
Base.:(==)(x::TypedMissing, ::WeakRef) = x
Base.:(==)(::WeakRef, y::TypedMissing) = y
Base.:(==)(x::TypedMissing, ::Missing) = TypedMissing()
Base.:(==)(::Missing, y::TypedMissing) = TypedMissing()

Base.isequal(x::TypedMissing, y::TypedMissing) = x === y
Base.isequal(x::TypedMissing, y::Missing) = x.kind == MissingKinds.NI
Base.isequal(x::Missing, y::TypedMissing) = y.kind == MissingKinds.NI
Base.isequal(x::TypedMissing, y::Any) = false
Base.isequal(x::Any, y::TypedMissing) = false

Base.:(<)(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.:(<)(x::TypedMissing, ::Any) = x
Base.:(<)(::Any, y::TypedMissing) = y
Base.:(<)(x::TypedMissing, ::Missing) = TypedMissing()
Base.:(<)(::Missing, y::TypedMissing) = TypedMissing()

Base.isless(x::TypedMissing, y::TypedMissing) =
    isless(x.kind, y.kind)
Base.isless(x::TypedMissing, y::Missing) = false
Base.isless(x::Missing, y::TypedMissing) = y.kind != MissingKinds.NI
Base.isless(x::TypedMissing, ::Any) = false
Base.isless(::Any, y::TypedMissing) = true

Base.isapprox(x::TypedMissing, y::TypedMissing; kwargs...) =
    x === y ? x : TypedMissing()
Base.isapprox(x::TypedMissing, ::Any; kwargs...) = x
Base.isapprox(::Any, y::TypedMissing; kwargs...) = y
Base.isapprox(x::TypedMissing, ::Missing; kwargs...) = TypedMissing()
Base.isapprox(::Missing, y::TypedMissing; kwargs...) = TypedMissing()

# Unary operators/functions
for f in (:(!), :(~), :(+), :(-), :(*), :(&), :(|), :(xor),
          :(zero), :(one), :(oneunit),
          :(isfinite), :(isinf), :(isodd),
          :(isinteger), :(isreal), :(isnan),
          :(iszero), :(transpose), :(adjoint), :(float), :(complex), :(conj),
          :(abs), :(abs2), :(iseven), :(ispow2),
          :(real), :(imag), :(sign), :(inv))
    @eval (Base.$f)(x::TypedMissing) = x
end
for f in (:(Base.zero), :(Base.one), :(Base.oneunit))
    @eval ($f)(x::Type{TypedMissing}) = TypedMissing()
    @eval function $(f)(::Type{Union{T, TypedMissing}}) where T
        T === Any && throw(MethodError($f, (Any,)))  # To prevent StackOverflowError
        $f(T)
    end
end
for f in (:(Base.float), :(Base.complex))
    @eval $f(::Type{TypedMissing}) = TypedMissing
    @eval function $f(::Type{Union{T, TypedMissing}}) where T
        T === Any && throw(MethodError($f, (Any,)))  # To prevent StackOverflowError
        Union{$f(T), TypedMissing}
    end
end

# Binary operators/functions
for f in (:(+), :(-), :(*), :(/), :(^), :(mod), :(rem))
    @eval begin
        # Scalar with missing
        (Base.$f)(x::TypedMissing, y::TypedMissing) =
            x === y ? x : TypedMissing()
        (Base.$f)(x::TypedMissing, ::Number) = x
        (Base.$f)(::Number, y::TypedMissing) = y
        (Base.$f)(x::TypedMissing, ::Missing) = TypedMissing()
        (Base.$f)(::Missing, y::TypedMissing) = TypedMissing()
    end
end

Base.div(x::TypedMissing, y::TypedMissing, r::RoundingMode) =
    x === y ? x : TypedMissing()
Base.div(x::TypedMissing, ::Number, r::RoundingMode) = x
Base.div(::Number, y::TypedMissing, r::RoundingMode) = y
Base.div(x::TypedMissing, ::Missing, r::RoundingMode) = TypedMissing()
Base.div(::Missing, y::TypedMissing, r::RoundingMode) = TypedMissing()

Base.min(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.min(x::TypedMissing, ::Any) = x
Base.min(::Any, y::TypedMissing) = y
Base.min(x::TypedMissing, ::Missing) = TypedMissing()
Base.min(::Missing, y::TypedMissing) = TypedMissing()
Base.max(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.max(x::TypedMissing, ::Any)     = x
Base.max(::Any,     y::TypedMissing) = y
Base.max(x::TypedMissing, ::Missing) = TypedMissing()
Base.max(::Missing, y::TypedMissing) = TypedMissing()

# Math functions
for f in (:sin, :cos, :tan, :asin, :atan, :acos,
          :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
          :exp, :exp2, :exp10, :expm1, :log, :log2, :log10, :log1p,
          :exponent, :sqrt, :cbrt)
    @eval Base.$(f)(x::TypedMissing) = x
end

for f in (:atan, :hypot, :log)
    @eval Base.$(f)(x::TypedMissing, y::Missing) = TypedMissing()
    @eval Base.$(f)(x::Missing, y::TypedMissing) = TypedMissing()
    @eval Base.$(f)(x::TypedMissing, y::TypedMissing) =
        x === y ? x : TypedMissing()
    @eval Base.$(f)(x::Number, y::TypedMissing) = y
    @eval Base.$(f)(x::TypedMissing, y::Number) = x
end

Base.clamp(x::TypedMissing, lo, hi) = x

# Rounding and related functions
Base.round(x::TypedMissing, ::RoundingMode=RoundNearest;
           sigdigits::Integer=0, digits::Integer=0, base::Integer=0) = x
Base.round(::Type{>:TypedMissing}, x::TypedMissing, ::RoundingMode=RoundNearest) = x
Base.round(::Type{T}, ::TypedMissing, ::RoundingMode=RoundNearest) where {T} =
    throw(MissingException("cannot convert a typed missing value to type $T: use Union{$T, TypedMissing} instead"))
Base.round(::Type{T}, x::Any, r::RoundingMode=RoundNearest) where {T>:TypedMissing} =
    round(Base.nonmissingtype_checked(T), x, r)
# to fix ambiguities
Base.round(::Type{T}, x::Rational{Tr}, r::RoundingMode=RoundNearest) where {T>:TypedMissing,Tr} =
    round(Base.nonmissingtype_checked(T), x, r)
Base.round(::Type{T}, x::Rational{Bool}, r::RoundingMode=RoundNearest) where {T>:TypedMissing} =
    round(Base.nonmissingtype_checked(T), x, r)

# Handle ceil, floor, and trunc separately as they have no RoundingMode argument
for f in (:(ceil), :(floor), :(trunc))
    @eval begin
        (Base.$f)(x::TypedMissing; sigdigits::Integer=0, digits::Integer=0, base::Integer=0) = x
        (Base.$f)(::Type{>:TypedMissing}, x::TypedMissing) = x
        (Base.$f)(::Type{T}, ::TypedMissing) where {T} =
            throw(MissingException("cannot convert a typed missing value to type $T: use Union{$T, TypedMissing} instead"))
        (Base.$f)(::Type{T}, x::Any) where {T>:TypedMissing} =
            $f(Base.nonmissingtype_checked(T), x)
        # to fix ambiguities
        (Base.$f)(::Type{T}, x::Rational) where {T>:TypedMissing} =
            $f(Base.nonmissingtype_checked(T), x)
    end
end

# to avoid ambiguity warnings
Base.:(^)(x::TypedMissing, ::Integer) = x

# Bit operators
Base.:(&)(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.:(&)(x::TypedMissing, y::Missing) = TypedMissing()
Base.:(&)(x::Missing, y::TypedMissing) = TypedMissing()
Base.:(&)(x::TypedMissing, y::Bool) = ifelse(y, x, false)
Base.:(&)(x::Bool, y::TypedMissing) = ifelse(x, y, false)
Base.:(&)(x::TypedMissing, ::Integer) = x
Base.:(&)(::Integer, y::TypedMissing) = y
Base.:(|)(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.:(|)(x::TypedMissing, y::Missing) = TypedMissing()
Base.:(|)(x::Missing, y::TypedMissing) = TypedMissing()
Base.:(|)(x::TypedMissing, y::Bool) = ifelse(y, true, x)
Base.:(|)(x::Bool, y::TypedMissing) = ifelse(x, true, y)
Base.:(|)(x::TypedMissing, ::Integer) = x
Base.:(|)(::Integer, y::TypedMissing) = y
Base.xor(x::TypedMissing, y::TypedMissing) =
    x === y ? x : TypedMissing()
Base.xor(x::TypedMissing, y::Missing) = TypedMissing()
Base.xor(x::Missing, y::TypedMissing) = TypedMissing()
Base.xor(x::TypedMissing, y::Bool) = x
Base.xor(x::Bool, y::TypedMissing) = y
Base.xor(x::TypedMissing, ::Integer) = x
Base.xor(::Integer, y::TypedMissing) = y

Base.:(*)(d::TypedMissing, x::Union{AbstractString,AbstractChar}) = d
Base.:(*)(d::Union{AbstractString,AbstractChar}, x::TypedMissing) = x

function Base.float(A::AbstractArray{Union{T, TypedMissing}}) where {T}
    U = typeof(float(zero(T)))
    convert(AbstractArray{Union{U, TypedMissing}}, A)
end
Base.float(A::AbstractArray{TypedMissing}) = A

## Redefinition of Base methods
## (temporary type piracy to use `ismissing` instead of `=== missing`,
## needed until Base is fixed to do the same)

# This method is the only one which cannot be submitted to Base as-is,
# an extensible mechanism would be needed instead (JuliaLang/julia#38241)
Base._promote_typesubtract(@nospecialize(a)) =
    Base.typesplit(a, Union{Nothing, Missing, TypedMissing})

function Base.iterate(itr::Base.SkipMissing, state...)
    y = iterate(itr.x, state...)
    y === nothing && return nothing
    item, state = y
    while ismissing(item)
        y = iterate(itr.x, state)
        y === nothing && return nothing
        item, state = y
    end
    item, state
end

Base.eachindex(itr::Base.SkipMissing) =
    Iterators.filter(i -> !ismissing(@inbounds(itr.x[i])), eachindex(itr.x))
Base.keys(itr::Base.SkipMissing) =
    Iterators.filter(i -> !ismissing(@inbounds(itr.x[i])), keys(itr.x))
Base.@propagate_inbounds function Base.getindex(itr::Base.SkipMissing, I...)
    v = itr.x[I...]
    ismissing(v) && throw(MissingException("the value at index $I is missing"))
    v
end

Base.mapreduce(f, op, itr::Base.SkipMissing{<:AbstractArray}) =
    Base._mapreduce(f, op, IndexStyle(itr.x),
                    eltype(itr.x) >: Missing || eltype(itr.x) >: TypedMissing ? itr : itr.x)

function Base._mapreduce(f, op, ::IndexLinear, itr::Base.SkipMissing{<:AbstractArray})
    A = itr.x
    # TODO: use TypedMissing() when applicable?
    ai = missing
    inds = LinearIndices(A)
    i = first(inds)
    ilast = last(inds)
    for outer i in i:ilast
        @inbounds ai = A[i]
        !ismissing(ai) && break
    end
    ismissing(ai) && return Base.mapreduce_empty(f, op, eltype(itr))
    a1::eltype(itr) = ai
    i == typemax(typeof(i)) && return Base.mapreduce_first(f, op, a1)
    i += 1
    ai = missing
    for outer i in i:ilast
        @inbounds ai = A[i]
        !ismissing(ai) && break
    end
    ismissing(ai) && return Base.mapreduce_first(f, op, a1)
    # We know A contains at least two non-missing entries: the result cannot be nothing
    something(Base.mapreduce_impl(f, op, itr, first(inds), last(inds)))
end

# Returns nothing when the input contains only missing values, and Some(x) otherwise
@noinline function Base.mapreduce_impl(f, op, itr::Base.SkipMissing{<:AbstractArray},
                                       ifirst::Integer, ilast::Integer, blksize::Int)
    A = itr.x
    if ifirst > ilast
        return nothing
    elseif ifirst == ilast
        @inbounds a1 = A[ifirst]
        if ismissing(ai)
            return nothing
        else
            return Some(Base.mapreduce_first(f, op, a1))
        end
    elseif ilast - ifirst < blksize
        # sequential portion
        # TODO: use TypedMissing() when applicable?
        ai = missing
        i = ifirst
        for outer i in i:ilast
            @inbounds ai = A[i]
            !ismissing(ai) && break
        end
        ismissing(ai) && return nothing
        a1 = ai::eltype(itr)
        i == typemax(typeof(i)) && return Some(Base.mapreduce_first(f, op, a1))
        i += 1
        ai = missing
        for outer i in i:ilast
            @inbounds ai = A[i]
            !ismissing(ai) && break
        end
        ismissing(ai) && return Some(Base.mapreduce_first(f, op, a1))
        a2 = ai::eltype(itr)
        i == typemax(typeof(i)) && return Some(op(f(a1), f(a2)))
        i += 1
        v = op(f(a1), f(a2))
        @simd for i = i:ilast
            @inbounds ai = A[i]
            if !ismissing(ai)
                v = op(v, f(ai))
            end
        end
        return Some(v)
    else
        # pairwise portion
        imid = ifirst + (ilast - ifirst) >> 1
        v1 = Base.mapreduce_impl(f, op, itr, ifirst, imid, blksize)
        v2 = Base.mapreduce_impl(f, op, itr, imid+1, ilast, blksize)
        if v1 === nothing && v2 === nothing
            return nothing
        elseif v1 === nothing
            return v2
        elseif v2 === nothing
            return v1
        else
            return Some(op(something(v1), something(v2)))
        end
    end
end

function Base.filter(f, itr::Base.SkipMissing{<:AbstractArray})
    y = similar(itr.x, eltype(itr), 0)
    for xi in itr.x
        if !ismissing(xi) && f(xi)
            push!(y, xi)
        end
    end
    y
end

if VERSION >= v"1.7.0"
    import Base: @coalesce
    macro coalesce(args...)
        expr = :(missing)
        for arg in reverse(args)
            expr = :((val = $arg); !ismissing(val) ? val : $expr)
        end
        return esc(:(let val; $expr; end))
    end
end

end # module TypedMissings


text:
include("WeightGeneration.jl")

using .WeightGeneration

using Formatting

in_filename  = "b.e11.B1850C5CN.f09_g16.005.pop.h.SST.100001-109912.nc"
out_filename = "test_SST.nc"
wgt_filename = "wgt_gx1v6_to_gx3v7.nc"

#=
WeightGeneration.convertFile(
    in_filename,
    out_filename,
    wgt_filename,
    varnames=("SHF",);
    xdim = "nlon",
    ydim = "nlat",
    zdim = "z_w",
)
=#
WeightGeneration.convertFile(
    in_filename,
    out_filename,
    wgt_filename,
    varnames=("SST",);
    xdim = "nlon",
    ydim = "nlat",
    zdim = "z_t",
)

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
centralpolymat_m_ks = [
    [
    (1, 4), # compile run
    (1, 15),
    (1, 25),
    (1, 50),
    (1, 75),
    (1, 100),
    (1, 125),
    (1, 150),
    (1, 175),
    ],
    [
    (3, 2), # compile run
    (4, 2),
    (4, 3),
    (4, 4),
    (4, 5),
    (4, 6),
    ],
    ]
centralpolymat_insts(ext::MatSpecExt) = [
    [(m, k, ext) for (m, k) in mks]
    for mks in centralpolymat_m_ks
    ]

insts = OrderedDict()
insts["nat"] = (nothing, vcat(centralpolymat_insts.([
    MatNegEntropy(),
    MatNegEntropyConj(),
    MatPower12(1.5),
    MatPower12Conj(1.5),
    ])...))
insts["ext"] = (nothing, vcat(centralpolymat_insts.([
    MatNegEntropyEigOrd(),
    MatNegEntropyConjEigOrd(),
    MatPower12EigOrd(1.5),
    MatPower12ConjEigOrd(1.5),
    ])...))
return (CentralPolyMatJuMP, insts)


text:
type SFWParameters{G <: AbstractFloat}
    a::G
    d::G
    e::G
    Ec::G
    rlim::G
    b::G
    q::G
    Jb::G
    ρs::G
    rs::G
    α::G
    β::G
    γ::G

    #secondary parameters
    Φs::G
    xlim::G
    Φ0::G
    Φlim::G

    adjEc::G
    adjJb::G
end

function SFWParameters(a::G, d::G, e::G, Ec::G, rlim::G, b::G, q::G, Jb::G, ρs::G, rs::G, α::G, β::G, γ::G; x0::G = 10.0^-8) where G <: AbstractFloat

    Φs = ρs * rs ^ 2
    xlim = rlim / rs
    adjEc = Φs * Ec
    adjJb = Jb * rs * sqrt(Φs)

    p = SFWParameters(a, d, e, Ec, rlim, b, q, Jb, ρs, rs, α, β, γ, Φs, xlim, 0.0, 0.0, adjEc, adjJb)

    p.Φ0 = gravitational_potential0(x0, p)
    p.Φlim = gravitational_potential(xlim, p)
    return p
end


text:
Dict(
     "T1" => "T1 - empty",
     "T2" => "T2 - empty",
     "T3" => "T3 - empty",
     "T4" => "T4 - empty",
     "T5" => "T5 - empty",
    )


text:
struct HTML <: Format end

# Parsing

function parse(io::IO, format::HTML)
    return parse(read(io, String), format)

end

function parse(s::String, ::HTML)
    return convert(XTree, Gumbo.parsehtml(s).root)
end

extensionformat(::Val{:html}) = HTML()

function Base.convert(::Type{XTree}, htmlnode::Gumbo.HTMLElement{S}) where S
    attrs = Dict{Symbol, Any}(Symbol(key) => val for (key, val) in htmlnode.attributes)
    return XNode(S, attrs, XTree[convert(XTree, c) for c in htmlnode.children])
end

function Base.convert(::Type{XTree}, htmltext::Gumbo.HTMLText)
    return XLeaf(htmltext.text)
end

# Rendering

function render!(io, x::XNode, format::HTML, ::Val)
    print(io, "<", x.tag)
    if !isempty(x.attributes)
        print(io, " ")
        for (name, attr) in x.attributes
            print(io, string(name), "=", attr, ' ')
        end
    end
    print(io, ">")
    foreach(child -> render!(io, child, format), children(x))
    print(io, "</", x.tag, ">")
end

const HTML_MIMES = [
    MIME"image/png"(),
    MIME"image/jpeg"(),
    MIME"image/svg+xml"(),
    #MIME"text/markdown"(),
    MIME"text/html"(),
    #MIME"text/latex"(),
    MIME"text/plain"(),
]

function render!(io, x::XLeaf, ::HTML)
    val = x[]

    for m in HTML_MIMES
        try
            if IJulia._showable(m, val)
                s = htmlstr(m, val)
                print(io, s)
                return
            end
        catch
            if m == MIME("text/plain")
                rethrow() # text/plain is required
            end
        end
    end
end

function htmlstr(mime::MIME, x)
    s = IJulia.limitstringmime(mime, x)
    return adapthtmlstr(mime, s)
end

adapthtmlstr(::MIME, s) = s

function adapthtmlstr(::MIME{Symbol("image/png")}, s)
    return """"""
end

function adapthtmlstr(::MIME{Symbol("image/jpeg")}, s)
    return """"""
end

function render!(io, x::XLeaf{<:AbstractString}, ::HTML)
    print(io, ansistringtohtml(x[]))
end

function ansistringtohtml(s)
    buf = IOBuffer()
    printer = HTMLPrinter(IOBuffer(s), root_tag="span")
    ANSIColoredPrinters.show_body(buf, printer)
    #show(buf, MIME"text/html"(), )
    return String(take!(buf))
end

IMAGEMIMES = [
    "image/jpeg",
    "image/svg+xml",
    "image/png",
    "image/gif",
]

formatextension(::HTML) = "html"

# Constants

const HTMLTAGS = [
    :DOCTYPE,
    :a,
    :abbr,
    :acronym,
    :address,
    :applet,
    :area,
    :article,
    :aside,
    :audio,
    :b,
    :base,
    :basefont,
    :bdi,
    :bdo,
    :big,
    :blockquote,
    :body,
    :br,
    :button,
    :canvas,
    :caption,
    :center,
    :cite,
    :code,
    :col,
    :colgroup,
    :data,
    :datalist,
    :dd,
    :del,
    :details,
    :dfn,
    :dialog,
    :dir,
    :div,
    :dl,
    :dt,
    :em,
    :embed,
    :fieldset,
    :figcaption,
    :figure,
    :font,
    :footer,
    :form,
    :frame,
    :frameset,
    :h1,
    :h2,
    :h3,
    :h4,
    :h5,
    :h6,
    :head,
    :header,
    :hr,
    :html,
    :HTML,
    :i,
    :iframe,
    :img,
    :input,
    :ins,
    :kbd,
    :label,
    :legend,
    :li,
    :link,
    :main,
    :map,
    :mark,
    :meta,
    :meter,
    :nav,
    :noframes,
    :noscript,
    :object,
    :ol,
    :optgroup,
    :option,
    :output,
    :p,
    :param,
    :picture,
    :pre,
    :progress,
    :q,
    :rp,
    :rt,
    :ruby,
    :s,
    :samp,
    :script,
    :section,
    :select,
    :small,
    :source,
    :span,
    :strike,
    :strong,
    :style,
    :sub,
    :summary,
    :sup,
    :svg,
    :table,
    :tbody,
    :td,
    :template,
    :textarea,
    :tfoot,
    :th,
    :thead,
    :time,
    :title,
    :tr,
    :track,
    :tt,
    :u,
    :ul,
    :var,
    :video,
    :wbr
]


text:
export VectorialTrajectory

const VectorialTrajectory = Trajectory{
    names,
    types,
    NamedTuple{names,trace_types},
} where {names,types,trace_types<:Tuple{Vararg{<:Vector}}}

"""
    VectorialTrajectory(;trace_name=trace_type ...)

Use `Vector` to store the traces.

# Example

```julia-repl
julia> t = VectorialTrajectory(;a=Int, b=Symbol)
0-element Trajectory{(:a, :b),Tuple{Int64,Symbol},NamedTuple{(:a, :b),Tuple{Array{Int64,1},Array{Symbol,1}}}}

julia> push!(t, a=0, b=:x)

julia> push!(t, a=1, b=:y)

julia> t
2-element Trajectory{(:a, :b),Tuple{Int64,Symbol},NamedTuple{(:a, :b),Tuple{Array{Int64,1},Array{Symbol,1}}}}:
 (a = 0, b = :x)
 (a = 1, b = :y)

julia> get_trace(t, :b)
2-element Array{Symbol,1}:
 :x
 :y

julia> pop!(t)

julia> t
1-element Trajectory{(:a, :b),Tuple{Int64,Symbol},NamedTuple{(:a, :b),Tuple{Array{Int64,1},Array{Symbol,1}}}}:
 (a = 0, b = :x)
```
"""
function VectorialTrajectory(; kwargs...)
    names = keys(kwargs.data)
    types = values(kwargs.data)
    trajectories =
        merge(NamedTuple(), (name, Vector{type}()) for (name, type) in zip(names, types))
    VectorialTrajectory{names,Tuple{types...},typeof(values(trajectories))}(trajectories)
end


text:
# indicator of the ball of matrices with (at most) a given rank

using LinearAlgebra
using TSVD

export IndBallRank

"""
**Indicator of rank ball**

    IndBallRank(r=1)

Returns the indicator function of the set of matrices of rank at most `r`:
```math
S = \\{ X : \\mathrm{rank}(X) \\leq r \\},
```
Parameter `r` must be a positive integer.
"""
struct IndBallRank{I <: Integer} <: ProximableFunction
    r::I
    function IndBallRank{I}(r::I) where {I <: Integer}
        if r <= 0
            error("parameter r must be a positive integer")
        else
            new(r)
        end
    end
end

is_set(f::IndBallRank) = true
is_prox_accurate(f::IndBallRank) = false

IndBallRank(r::I=1) where {I <: Integer} = IndBallRank{I}(r)

function (f::IndBallRank)(x::AbstractArray{T, 2}) where {R <: Real, T <: RealOrComplex{R}}
    maxr = minimum(size(x))
    if maxr <= f.r return R(0) end
    U, S, V = tsvd(x, f.r+1)
    # the tolerance in the following line should be customizable
    if S[end]/S[1] <= 1e-7
        return R(0)
    end
    return R(Inf)
end

function prox!(y::AbstractMatrix{T}, f::IndBallRank, x::AbstractMatrix{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    maxr = minimum(size(x))
    if maxr <= f.r
        y .= x
        return R(0)
    end
    U, S, V = tsvd(x, f.r)
    # TODO: the order of the following matrix products should depend on the shape of x
    M = S .* V'
    mul!(y, U, M)
    return R(0)
end

fun_name(f::IndBallRank) = "indicator of the set of rank-r matrices"
fun_dom(f::IndBallRank) = "AbstractArray{Real,2}, AbstractArray{Complex,2}"
fun_expr(f::IndBallRank) = "x ↦ 0 if rank(x) ⩽ r, +∞ otherwise"
fun_params(f::IndBallRank) = "r = $(f.r)"

function prox_naive(f::IndBallRank, x::AbstractMatrix{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    maxr = minimum(size(x))
    if maxr <= f.r
        y = x
        return y, R(0)
    end
    F = svd(x)
    y = F.U[:,1:f.r]*(Diagonal(F.S[1:f.r])*F.V[:,1:f.r]')
    return y, R(0)
end


text:
"""
    ramsey()
``ramsey'' is an MINLPLib test instance with binary and continuous variables
(see http://www.minlplib.org/ramsey.html). The best known primal feasible
objective value for this problem is 2.48746864. This example contains 22 nonlinear
terms, 11 are x**0.25 and remaining are log(x)
"""

function post_ramsey_linear_constraints!(m)
@variable(m, objvar)
    x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]
    @variable(m, x[x_Idx])
    setlowerbound(x[1], 3.0)
    setupperbound(x[1], 3.0)
    setlowerbound(x[2], 3.0)
    setlowerbound(x[3], 3.0)
    setlowerbound(x[4], 3.0)
    setlowerbound(x[5], 3.0)
    setlowerbound(x[6], 3.0)
    setlowerbound(x[7], 3.0)
    setlowerbound(x[8], 3.0)
    setlowerbound(x[9], 3.0)
    setlowerbound(x[10], 3.0)
    setlowerbound(x[11], 3.0)
    setlowerbound(x[12], 0.95)
    setlowerbound(x[13], 0.95)
    setlowerbound(x[14], 0.95)
    setlowerbound(x[15], 0.95)
    setlowerbound(x[16], 0.95)
    setlowerbound(x[17], 0.95)
    setlowerbound(x[18], 0.95)
    setlowerbound(x[19], 0.95)
    setlowerbound(x[20], 0.95)
    setlowerbound(x[21], 0.95)
    setlowerbound(x[22], 0.95)
    setlowerbound(x[23], 0.05)
    setupperbound(x[23], 0.05)
    setlowerbound(x[24], 0.05)
    setupperbound(x[24], 0.0575)
    setlowerbound(x[25], 0.05)
    setupperbound(x[25], 0.066125)
    setlowerbound(x[26], 0.05)
    setupperbound(x[26], 0.07604375)
    setlowerbound(x[27], 0.05)
    setupperbound(x[27], 0.0874503125)
    setlowerbound(x[28], 0.05)
    setupperbound(x[28], 0.100567859375)
    setlowerbound(x[29], 0.05)
    setupperbound(x[29], 0.11565303828125)
    setlowerbound(x[30], 0.05)
    setupperbound(x[30], 0.133000994023437)
    setlowerbound(x[31], 0.05)
    setupperbound(x[31], 0.152951143126953)
    setlowerbound(x[32], 0.05)
    setupperbound(x[32], 0.175893814595996)
    setlowerbound(x[33], 0.05)
    setupperbound(x[33], 0.202277886785395)

    # ----- Constraints ----- #
    @constraint(m, e12, -x[1]+x[2]-x[23] == 0.0)
    @constraint(m, e13, -x[2]+x[3]-x[24] == 0.0)
    @constraint(m, e14, -x[3]+x[4]-x[25] == 0.0)
    @constraint(m, e15, -x[4]+x[5]-x[26] == 0.0)
    @constraint(m, e16, -x[5]+x[6]-x[27] == 0.0)
    @constraint(m, e17, -x[6]+x[7]-x[28] == 0.0)
    @constraint(m, e18, -x[7]+x[8]-x[29] == 0.0)
    @constraint(m, e19, -x[8]+x[9]-x[30] == 0.0)
    @constraint(m, e20, -x[9]+x[10]-x[31] == 0.0)
    @constraint(m, e21, -x[10]+x[11]-x[32] == 0.0)
    @constraint(m, e22, 0.03*x[11]-x[33] <= 0.0)

    println("creating lifted variables and constraints")
    @variable(m, y[1:22])
    @constraint(m, y[1] == 3^0.25)
    @constraint(m, e1, 0.759835685651593*y[1]-x[12]-x[23] == 0.0)
    @constraint(m, e2, 0.77686866556676*y[2]-x[13]-x[24] == 0.0)
    @constraint(m, e3, 0.794283468039448*y[3]-x[14]-x[25] == 0.0)
    @constraint(m, e4, 0.812088652256959*y[4]-x[15]-x[26] == 0.0)
    @constraint(m, e5, 0.830292969275008*y[5]-x[16]-x[27] == 0.0)
    @constraint(m, e6, 0.848905366318769*y[6]-x[17]-x[28] == 0.0)
    @constraint(m, e7, 0.867934991180342*y[7]-x[18]-x[29] == 0.0)
    @constraint(m, e8, 0.88739119671479*y[8]-x[19]-x[30] == 0.0)
    @constraint(m, e9, 0.907283545436972*y[9]-x[20]-x[31] == 0.0)
    @constraint(m, e10, 0.92762181422141*y[10]-x[21]-x[32] == 0.0)
    @constraint(m, e11, 0.948415999107521*y[11]-x[22]-x[33] == 0.0)
    @constraint(m, e23, -(0.95*y[12]+0.9025*y[13]+0.857375*y[14]+0.81450625*y[15]+0.7737809375*y[16]+0.735091890625*y[17]+0.69833729609375*y[18]+0.663420431289062*y[19]+0.630249409724609*y[20]+0.598736939238379*y[21]+11.3760018455292*y[22])-objvar == 0.0)

    # ----- Objective ----- #
    @objective(m, Min, objvar)

    y_indexes_1 = collect(2:11)
    y_indexes_2 = collect(12:22)
    x_indexes_1 = collect(2:11)
    x_indexes_2 = collect(12:22)
    partitions_1 = []
    partitions_2 = []
    for index in x_indexes_1
        ub = 1e2
        (has_upper_bound(x[index])) && (ub = upper_bound(x[index]))
        push!(partitions_1, [lower_bound(x[index]), ub])
    end
    for index in x_indexes_2
        ub = 1e2
        (has_upper_bound(x[index])) && (ub = upper_bound(x[index]))
        push!(partitions_2, [lower_bound(x[index]), ub])
    end

    return y_indexes_1, x_indexes_1, partitions_1, y_indexes_2, x_indexes_2, partitions_2
end

function ramsey()
    best_known_objective = -2.48746864
    cplex_optimizer =
        JuMP.optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0)
    error_tolerances = [NaN64, 1e-1, 1e-2]
    # base partition holds the inflection points and the end-points for the trigonometric functions
    milp_statistics = zeros(length(error_tolerances), 5)
    lp_statistics = zeros(length(error_tolerances), 5)
    func_1 = x -> x^0.25
    func_2 = x -> log(x)
    # solve MILP for various error tolerances
    for i = 1:length(error_tolerances)
        error_tolerance = error_tolerances[i]
        milp = Model(cplex_optimizer)
        y_indexes_1, x_indexes_1, base_partition_1,
        y_indexes_2, x_indexes_2, base_partition_2 = post_ramsey_linear_constraints!(milp)
        # construct MILP relaxations
        for j = 1:length(y_indexes_1)
            var, x_index, y_index = add_milp_relaxation(
                func_1,
                base_partition_1[j],
                error_tolerance,
                milp,
                name = "1_" * string(j),
            )
            @constraint(milp, var[x_index] == milp[:x][x_indexes_1[j]])
            @constraint(milp, var[y_index] == milp[:y][y_indexes_1[j]])
        end
        for j = 1:length(y_indexes_2)
            var, x_index, y_index = add_milp_relaxation(
                func_2,
                base_partition_2[j],
                error_tolerance,
                milp,
                name = "2_" * string(j),
            )
            @constraint(milp, var[x_index] == milp[:x][x_indexes_2[j]])
            @constraint(milp, var[y_index] == milp[:y][y_indexes_2[j]])
        end
        _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(milp)
        relaxation_objective = objective_value(milp)
        @assert relaxation_objective ≤ best_known_objective
        relative_gap =
            abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
        milp_statistics[i, :] = [
            error_tolerance,
            best_known_objective,
            relaxation_objective,
            relative_gap,
            round(solve_time, digits = 2),
        ]
    end

    for i = 1:length(error_tolerances)
        error_tolerance = error_tolerances[i]
        lp = Model(cplex_optimizer)
        y_indexes_1, x_indexes_1, base_partition_1,
        y_indexes_2, x_indexes_2, base_partition_2 = post_ramsey_linear_constraints!(lp)
        # construct MILP relaxations
        for j = 1:length(y_indexes_1)
            var, x_index, y_index = add_lp_relaxation(
                func_1,
                base_partition_1[j],
                error_tolerance,
                lp,
                name = "1_" * string(j),
            )
            @constraint(lp, var[x_index] == lp[:x][x_indexes_1[j]])
            @constraint(lp, var[y_index] == lp[:y][y_indexes_1[j]])
        end
        for j = 1:length(y_indexes_2)
            var, x_index, y_index = add_lp_relaxation(
                func_2,
                base_partition_2[j],
                error_tolerance,
                lp,
                name = "2_" * string(j),
            )
            @constraint(lp, var[x_index] == lp[:x][x_indexes_2[j]])
            @constraint(lp, var[y_index] == lp[:y][y_indexes_2[j]])
        end
        _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(lp)
        relaxation_objective = objective_value(lp)
        @assert relaxation_objective ≤ best_known_objective
        relative_gap =
            abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
        lp_statistics[i, :] = [
            error_tolerance,
            best_known_objective,
            relaxation_objective,
            relative_gap,
            round(solve_time, digits = 2),
        ]
    end

    open("csv/ramsey_milp.csv", "w") do io
        writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
        writedlm(io, milp_statistics, ",")
    end

    open("csv/ramsey_lp.csv", "w") do io
        writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
        writedlm(io, lp_statistics, ",")
    end

end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
"""
    Models.GaussianShells()

2-D Gaussian shells centered at `[-3.5, 0]` and `[3.5, 0]` with a radius of 2 and a shell width of 0.1

# Examples
```jldoctest
julia> model, lnZ = Models.GaussianShells();

julia> lnZ
-1.75
```
"""
function GaussianShells()
    μ1 = [-3.5, 0]
    μ2 = [3.5, 0]

    prior(X) = 12 .* X .- 6
    loglike(X) = logaddexp(logshell(X, μ1), logshell(X, μ2))

    lnZ = -1.75
    return NestedModel(loglike, prior), lnZ
end

function logshell(X, μ, radius=2, width=0.1)
    d = LinearAlgebra.norm(X - μ)
    norm = -log(sqrt(2 * π * width^2))
    return norm - (d - radius)^2 / (2 * width^2)
end


text:
function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigenvectors::Array{Complex{T}}, psource::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}};
        basis_order::Int = 2,
        kws...
    ) where T

    if size(eigenvectors)[end] > 1
        @warn "The effective wavenumber: $k_eff has more than one eigenvector. For plane-waves this case has not been fully implemented"
    end

    scale_number_density = one(T) - one(T) / material.numberofparticles

    direction_eff = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, material.shape.normal)

    θ_eff = transmission_angle(direction_eff, material.shape.normal)
    θin  = transmission_angle(psource, material)

    kcos_in = (ω / psource.medium.c) * dot(- conj(material.shape.normal), psource.direction)
    kcos_eff = k_eff * dot(- conj(material.shape.normal), direction_eff)

    extinction_matrix = T(2) .* transpose(vec(
        [
            exp(im*n*(θin - θ_eff)) * scale_number_density * number_density(s)
        for n = -basis_order:basis_order, s in material.species]
    ))

    forcing = [im * field(psource,zeros(T,2),ω) * kcos_in * (kcos_eff - kcos_in)]

    # where extinction_matrix * eigenvectors * α = forcing
    α = (extinction_matrix * reshape(eigenvectors,(:,size(eigenvectors)[end]))) \ forcing

    return α

end

function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Complex{T}}, psource::PlaneSource{T,3,1,Acoustic{T,3}}, material::Material{3,Halfspace{T,3}};
        basis_order::Int = 2,
        kws...
    ) where T

    if size(eigvectors)[end] > 1
        @warn "The effective wavenumber: $k_eff has more than one eigenvector. For plane-waves this case has not been fully implemented"
    end

    scale_number_density = one(T) - one(T) / material.numberofparticles

    # First we calculate the outward pointing normal
    n = material.shape.normal;
    n = - n .* sign(real(dot(n,psource.direction)));

    # calculate the distances from the origin to the face of the halfspace
    Z0 = dot(-n, material.shape.origin)

    # next the components of the wavenumbers in the direction of the inward normal
    direction = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, n)

    k = (ω / psource.medium.c)
    kz = k * dot(-conj(n), psource.direction)
    k_effz = k_eff * dot(-conj(n), direction)

    rθφ = cartesian_to_radial_coordinates(psource.direction);

    Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
    lm_to_n = lm_to_spherical_harmonic_index

    I1(F::Array{Complex{T}}, k_effz::Complex{T}) = sum(
        - F[lm_to_n(dl,dm),j,p] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dl *
        exp(im * (k_effz - kz)*(Z0 + outer_radius(material.species[j]))) *
        scale_number_density * number_density(material.species[j]) / (k*kz * (kz - k_effz))
    for p = 1:size(F,3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species))

    forcing = - field(psource,zeros(T,3),ω)

    α = I1(eigvectors,k_effz) \ forcing

    return α

end

function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{Complex{T}}, eigvectors2::Array{Complex{T}}, psource::PlaneSource{T,3,1,Acoustic{T,3}}, material::Material{3,Plate{T,3}};
        basis_order::Int = 2,
        kws...
    ) where T

    if size(eigvectors1)[end] > 1 || size(eigvectors2)[end] > 1
        @warn "The effective wavenumber: $k_eff has more than one eigenvector. For plane-waves this case has not been fully implemented"
    end

    scale_number_density = one(T) - one(T) / material.numberofparticles

    # First we calculate the outward pointing normal
    n = material.shape.normal;
    n = - n .* sign(real(dot(n,psource.direction)));

    # calculate the distances from the origin to the two faces of the plate
    Z0 = dot(-n, material.shape.origin)
    Z1 = Z0 - material.shape.width / 2
    Z2 = Z0 + material.shape.width / 2

    # next the components of the wavenumbers in the direction of the inward normal
    k = (ω / psource.medium.c)
    k1 = k_eff;
    k2 = -k_eff;

    direction1 = transmission_direction(k1, k * psource.direction, n)
    direction2 = direction1

    kz = k * dot(-conj(n), psource.direction)
    k1_effz = k1 * dot(-conj(n), direction1)
    k2_effz = k2 * dot(-conj(n), direction2)

    rθφ = cartesian_to_radial_coordinates(psource.direction);

    Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
    lm_to_n = lm_to_spherical_harmonic_index

    I1(F::Array{Complex{T}}, k_effz::Complex{T}) = sum(
        - F[lm_to_n(dl,dm),j,p] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dl *
        exp(im * (k_effz - kz)*(Z1 + outer_radius(material.species[j]))) *
        scale_number_density * number_density(material.species[j]) / (k*kz * (kz - k_effz))
    for p = 1:size(F,3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species))

    I2(F::Array{Complex{T}}, k_effz::Complex{T}) = sum(
        - F[lm_to_n(dl,dm),j,p] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dm *
        exp(im * (k_effz + kz)*(Z2 - outer_radius(material.species[j]))) *
        scale_number_density * number_density(material.species[j]) / (k*kz * (kz + k_effz))
    for p = 1:size(F,3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species))

    MIs = [I1(eigvectors1,k1_effz) I1(eigvectors2,k2_effz);
           I2(eigvectors1,k1_effz) I2(eigvectors2,k2_effz)]

    forcing = [ - field(psource,zeros(T,3),ω), T(0)]

    α = MIs \ forcing

    return α

end


text:
module TestDebugEval
using Debug
import Debug.Node, Debug.Scope, Debug.@instrument, Debug.debug_eval

firstline = -1
function trap(node::Node, scope::Scope)
    global firstline = (firstline == -1) ? node.loc.line : firstline
    line = node.loc.line - firstline + 1
    print(line, ":")

    if (line == 2); debug_eval(scope, :(x += 1)) end
    if (line == 7); debug_eval(scope, :((y,z) = (z,y))) end

    if (line >  1); print("\tx = ", debug_eval(scope, :x)) end
    if (line == 3); print("\tk = ", debug_eval(scope, :k)) end
    if (line >  5)
        print("\ty = ", debug_eval(scope, :y))
        print("\tz = ", debug_eval(scope, :z))
    end
    println()
end

@instrument trap function f(n)
    x = 0       # line 1
    for k=1:n   # line 2
        x += k  # line 3
    end         # line 4
    y, z = 1:2  # line 5
    x = x*x     # line 6
    x           # line 7
end

f(3)

end


text:
using OpenStreetMapX

pth = joinpath(dirname(pathof(OpenStreetMapX)),"..","test","data","reno_east3.osm")
map_data =  OpenStreetMapX.get_map_data(pth,use_cache = false);

using Random
Random.seed!(0);

origin = generate_point_in_bounds(map_data);
destination = generate_point_in_bounds(map_data);

pointA = point_to_nodes(origin, map_data)
pointB = point_to_nodes(destination, map_data)

shortest_route1, shortest_distance1, shortest_time1 = shortest_route(map_data, pointA, pointB)
fastest_route1, fastest_distance1, fastest_time1 = fastest_route(map_data, pointA, pointB)

println("shortest_route nodes: ",shortest_route1)
println("fastest route nodes: ",fastest_route1)

### Create this file if you want to test routing with Google API
### The file should only contain your Google API key
google_api_file = joinpath("some_folder","googleapi.key")

if isfile(google_api_file)
    google_api_key = readlines(google_api_file)[1]
    google_route = get_google_route(pointA, pointB,map_data,google_api_key)[1]
    println("Google API route nodes : ",google_route)
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, WildBootTests

makedocs(sitename="WildBootTests.jl",
         authors="David Roodman",
         pages = [
                  "Overview" => "index.md",
                  "OLS examples" => "OLSexamples.md",
                  "IV/2SLS examples" => "IVexamples.md",
                  "Public functions and types" => "exported.md"
                 ]
        )

deploydocs(
  deploy_config = Documenter.GitHubActions(),
  repo = "github.com/droodman/WildBootTests.jl.git",
)