text:
# the demo was given making use of julia 0.6.4

import ArchGDAL
using PerceptualColourMaps
using GeoInterface
using ColorTypes

GDAL.allregister()

#--- From the lake rough bbox, get image

jsonstring = """{"type":"Polygon","coordinates":[[[-3.01379,54.52648],[-3.00396,54.52648],[-3.00396,54.53042],[-3.01379,54.53042],[-3.01379,54.52648]]]}"""

geom = ArchGDAL.importEPSG(4326) do wgs84
    ArchGDAL.importEPSG(27700) do osgb1936
        ArchGDAL.createcoordtrans(wgs84, osgb1936) do transform
            geom = ArchGDAL.fromJSON(jsonstring)
            ArchGDAL.transform!(geom, transform)
end end end

ArchGDAL.issimple(geom)

geotransform = ArchGDAL.read("lake-district/ortho-20cm.tif") do dataset
    ArchGDAL.getgeotransform(dataset)
end

function rasterindex(centerpoint, geotransform)
    x_min = geotransform[1]
    x_res = geotransform[2]
    y_max = geotransform[4]
    y_res = geotransform[6]
    x, y = centerpoint
    row = Int(fld(y - y_max, y_res)) + 1
    col = Int(fld(x - x_min, x_res)) + 1
    row, col
end

function rasterindex(bbox::GDAL.OGREnvelope, geotransform)
    topleft = (bbox.MinX, bbox.MaxY)
    bottomright = (bbox.MaxX, bbox.MinY)
    r1, c1 = rasterindex(topleft, geotransform)
    r2, c2 = rasterindex(bottomright, geotransform)
    r1:r2, c1:c2
end

bbox = ArchGDAL.getenvelope(geom)
rows, cols = rasterindex(bbox, geotransform)

subortho = ArchGDAL.read("lake-district/ortho-20cm.tif") do dataset
    ArchGDAL.imread(RGB, dataset, (1,2,3), rows, cols)
end

collect(subortho)

#--- Get the three faced peak from the DTM

# indices of peak
rows = 3001:4500
cols = 4001:7000

# read DTM of this area
subdem = ArchGDAL.read("lake-district/dtm-1m.tif") do dataset
    PermutedDimsArray(ArchGDAL.read(dataset, 1, rows, cols), (2,1))
end

# rolling down the hill
function roll!(canvas, dem, i, j)
    count = 1
    while checkbounds(Bool, dem, i-1:i+1, j-1:j+1) && count < 100_000
        surroundings = dem[i-1:i+1, j-1:j+1]
        _, linidx = findmin(surroundings)
        local_i, local_j = ind2sub(surroundings, linidx)
        i += local_i - 2
        j += local_j - 2
        count += 1
        canvas[i-1:i+1, j-1:j+1] = 0.0
    end
end

# hillshade as background map
hillshade = Gray.(relief(collect(subdem), 45.0, 45.0, 1.0))
canvas = copy(hillshade)
roll!(canvas, subdem, 1215, 755); canvas
roll!(canvas, subdem, 1175, 855); canvas
roll!(canvas, subdem, 1240, 800); canvas

#---


text:
using VectorizationBase: vnoaliasstore!

@inline vstoreadditivereduce!(args...) = vnoaliasstore!(args...)
@inline vstoremultiplicativevereduce!(args...) = vnoaliasstore!(args...)
@inline function vstoreadditivereduce!(ptr::VectorizationBase.AbstractStridedPointer, v::VectorizationBase.SVec, i::Tuple{Vararg{Union{Integer,Static}}})
    vnoaliasstore!(ptr, SIMDPirates.vsum(v), i)
end
@inline function vstoreadditivereduce!(ptr::VectorizationBase.AbstractStridedPointer, v::VectorizationBase.SVec, i::Tuple{Vararg{Union{Integer,Static}}}, m::VectorizationBase.Mask)
    vnoaliasstore!(ptr, SIMDPirates.vsum(v), i, m)
end
@inline function vstoremultiplicativevereduce!(ptr::VectorizationBase.AbstractStridedPointer, v::VectorizationBase.SVec, i::Tuple{Vararg{Union{Integer,Static}}})
    vnoaliasstore!(ptr, SIMDPirates.vprod(v), i)
end
@inline function vstoremultiplicativevereduce!(ptr::VectorizationBase.AbstractStridedPointer, v::VectorizationBase.SVec, i::Tuple{Vararg{Union{Integer,Static}}}, m::VectorizationBase.Mask)
    vnoaliasstore!(ptr, SIMDPirates.vprod(v), i, m)
end

function storeinstr(op::Operation, vectorized::Symbol)
    # defaultstoreop = :vstore!
    defaultstoreop = :vnoaliasstore!
    vectorized ∉ reduceddependencies(op) && return lv(defaultstoreop)
    vectorized ∈ loopdependencies(op) && return lv(defaultstoreop)
    # vectorized is not a loopdep, but is a reduced dep
    opp = first(parents(op))
    while vectorized ∉ loopdependencies(opp)
        oppold = opp
        for oppp ∈ parents(opp)
            if vectorized ∈ reduceddependencies(oppp)
                @assert opp !== oppp "More than one parent is a reduction over the vectorized variable."
                opp = oppp
            end
        end
        @assert opp !== oppold "Failed to find any parents "
    end
    instr_class = reduction_instruction_class(instruction(opp))
    instr = if instr_class === ADDITIVE_IN_REDUCTIONS
        :vstoreadditivereduce!
    elseif instr_class === MULTIPLICATIVE_IN_REDUCTIONS
        :vstoremultiplicativevereduce!
    else #FIXME
        defaultstoreop
    end
    lv(instr)
end

# const STOREOP = :vstore!
# variable_name(op::Operation, ::Nothing) = mangledvar(op)
# variable_name(op::Operation, suffix) = Symbol(mangledvar(op), suffix, :_)
# # variable_name(op::Operation, suffix, u::Int) = (n = variable_name(op, suffix); u < 0 ? n : Symbol(n, u))
function reduce_range!(q::Expr, toreduct::Symbol, instr::Instruction, Uh::Int, Uh2::Int)
    if 2Uh == Uh2
        for u ∈ 0:2:Uh2-1
            instrexpr = callexpr(instr)
            push!(instrexpr.args, Symbol(toreduct, u))
            push!(instrexpr.args, Symbol(toreduct, u + 1))
            push!(q.args, Expr(:(=), Symbol(toreduct, (u>>>1)), instrexpr))
        end
    else
        for u ∈ Uh:Uh2-2
            tru = Symbol(toreduct, u - Uh)
            instrexpr = callexpr(instr)
            push!(instrexpr.args, tru)
            push!(instrexpr.args, Symbol(toreduct, u))
            push!(q.args, Expr(:(=), tru, instrexpr))
        end
        for u ∈ 2Uh:Uh2-1
            tru = Symbol(toreduct, u - 2Uh)
            instrexpr = callexpr(instr)
            push!(instrexpr.args, tru)
            push!(instrexpr.args, Symbol(toreduct, u))
            push!(q.args, Expr(:(=), tru, instrexpr))
        end
    end
end
function reduce_range!(q::Expr, ls::LoopSet, Ulow::Int, Uhigh::Int)
    for or ∈ ls.outer_reductions
        op = ls.operations[or]
        var = mangledvar(op)
        instr = Instruction(reduction_to_single_vector(op.instruction))
        reduce_range!(q, var, instr, Ulow, Uhigh)
    end
end

function reduce_expr!(q::Expr, toreduct::Symbol, instr::Instruction, U::Int)
    U == 1 && return nothing
    @assert U > 1 "U = $U somehow < 1"
    instr = Instruction(reduction_to_single_vector(instr))
    Uh2 = U
    iter = 0
    while true # combine vectors
        Uh = Uh2 >> 1
        reduce_range!(q, toreduct, instr, Uh, Uh2)
        Uh == 1 && break
        Uh2 = Uh
        iter += 1; iter > 4 && throw("Oops! This seems to be excessive unrolling.")
    end
    nothing
end

# pvariable_name(op::Operation, ::Nothing) = mangledvar(first(parents(op)))
# pvariable_name(op::Operation, ::Nothing, ::Symbol) = mangledvar(first(parents(op)))
# pvariable_name(op::Operation, suffix) = Symbol(pvariable_name(op, nothing), suffix, :_)
# function pvariable_name(op::Operation, suffix, tiled::Symbol)
#     parent = first(parents(op))
#     mname = mangledvar(parent)
#     tiled ∈ loopdependencies(parent) ? Symbol(mname, suffix, :_) : mname
# end

function lower_conditionalstore_scalar!(
    q::Expr, op::Operation, ua::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned}, inds_calc_by_ptr_offset::Vector{Bool}
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    mvar, opu₁, opu₂ = variable_name_and_unrolled(first(parents(op)), u₁loopsym, u₂loopsym, suffix)
    # var = pvariable_name(op, suffix, tiled)
    cond = last(parents(op))
    condvar, condu₁ = condvarname_and_unroll(cond, u₁loopsym, u₂loopsym, suffix, opu₂)
    loopdeps = loopdependencies(op)
    opu₁ = u₁loopsym ∈ loopdeps
    ptr = vptr(op)
    for u ∈ 0:u₁-1
        varname = varassignname(mvar, u, opu₁)
        condvarname = varassignname(condvar, u, condu₁)
        td = UnrollArgs(ua, u)
        push!(q.args, Expr(:&&, condvarname, Expr(:call, storeinstr(op, vectorized), ptr, varname, mem_offset_u(op, td, inds_calc_by_ptr_offset))))
    end
    nothing
end
function lower_conditionalstore_vectorized!(
    q::Expr, op::Operation, ua::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned}, isunrolled::Bool, inds_calc_by_ptr_offset::Vector{Bool}
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    loopdeps = loopdependencies(op)
    @assert vectorized ∈ loopdeps
    mvar, opu₁, opu₂ = variable_name_and_unrolled(first(parents(op)), u₁loopsym, u₂loopsym, suffix)
    # var = pvariable_name(op, suffix, tiled)
    if isunrolled
        umin = 0
        U = u₁
    else
        umin = -1
        U = 0
    end
    ptr = vptr(op)
    vecnotunrolled = vectorized !== u₁loopsym
    cond = last(parents(op))
    condvar, condu₁ = condvarname_and_unroll(cond, u₁loopsym, u₂loopsym, suffix, opu₂)
    # @show parents(op) cond condvar
    for u ∈ 0:U-1
        td = UnrollArgs(ua, u)
        name, mo = name_memoffset(mvar, op, td, opu₁, inds_calc_by_ptr_offset)
        condvarname = varassignname(condvar, u, condu₁)
        instrcall = Expr(:call, storeinstr(op, vectorized), ptr, name, mo)
        if mask !== nothing && (vecnotunrolled || u == U - 1)
            push!(instrcall.args, Expr(:call, :&, condvarname, mask))
        else
            push!(instrcall.args, condvarname)
        end
        push!(q.args, instrcall)
    end
end

function lower_store_scalar!(
    q::Expr, op::Operation, ua::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned}, inds_calc_by_ptr_offset::Vector{Bool}
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    mvar, opu₁, opu₂ = variable_name_and_unrolled(first(parents(op)), u₁loopsym, u₂loopsym, suffix)
    ptr = vptr(op)
    # var = pvariable_name(op, suffix, tiled)
    # parentisunrolled = unrolled ∈ loopdependencies(first(parents(op)))
    for u ∈ 0:u₁-1
        varname = varassignname(mvar, u, opu₁)
        td = UnrollArgs(ua, u)
        push!(q.args, Expr(:call, storeinstr(op, vectorized), ptr, varname, mem_offset_u(op, td, inds_calc_by_ptr_offset)))
    end
    nothing
end
function lower_store_vectorized!(
    q::Expr, op::Operation, ua::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned}, isunrolled::Bool, inds_calc_by_ptr_offset::Vector{Bool}
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    loopdeps = loopdependencies(op)
    @assert vectorized ∈ loopdeps
    mvar, opu₁, opu₂ = variable_name_and_unrolled(first(parents(op)), u₁loopsym, u₂loopsym, suffix)
    ptr = vptr(op)
    # var = pvariable_name(op, suffix, tiled)
    # parentisunrolled = unrolled ∈ loopdependencies(first(parents(op)))
    if isunrolled
        umin = 0
        U = u₁
    else
        umin = -1
        U = 0
    end
    vecnotunrolled = vectorized !== u₁loopsym
    for u ∈ umin:U-1
        td = UnrollArgs(ua, u)
        name, mo = name_memoffset(mvar, op, td, opu₁, inds_calc_by_ptr_offset)
        instrcall = Expr(:call, storeinstr(op, vectorized), ptr, name, mo)
        if mask !== nothing && (vecnotunrolled || u == U - 1)
            push!(instrcall.args, mask)
        end
        push!(q.args, instrcall)
    end
end
function lower_store!(
    q::Expr, ls::LoopSet, op::Operation, ua::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned} = nothing
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    isunrolled₁ = isu₁unrolled(op) #u₁loopsym ∈ loopdependencies(op)
    # isunrolled₂ = isu₂unrolled(op)
    inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)
    # if isunrolled₁ & ((!isnothing(suffix)) && isu₂unrolled(op))
    #     u₁ind = findfirst(isequal(u₁loopsym), loopdependencies(op))::Int
    #     u₂ind = findfirst(isequal(u₁loopsym), loopdependencies(op))::Int
    #     if inds_calc_by_ptr_offset[u₁ind] && inds_calc_by_ptr_offset[u₂ind]
    #         if u₁ind < u₂ind
    #             zero_offset = u₂ind
    #             gespind = gensym("gesp")
    #             # gesp by suffix
    #         else
    #             zero_offset = u₁ind
    #             if iszero(suffix) # gesp for each

    #             end
    #         end
    #     end
    # else
    # end
    ua = UnrollArgs(ua, isunrolled₁ ? u₁ : 1)
    if instruction(op).instr !== :conditionalstore!
        if vectorized ∈ loopdependencies(op)
            lower_store_vectorized!(q, op, ua, mask, isunrolled₁, inds_calc_by_ptr_offset)
        else
            lower_store_scalar!(q, op, ua, mask, inds_calc_by_ptr_offset)
        end
    else
        if vectorized ∈ loopdependencies(op)
            lower_conditionalstore_vectorized!(q, op, ua, mask, isunrolled₁, inds_calc_by_ptr_offset)
        else
            lower_conditionalstore_scalar!(q, op, ua, mask, inds_calc_by_ptr_offset)
        end
    end
end

ext: jl


text:
include("../../module/CustomStencil.jl")
include("../../misc/misc.jl")
include("../../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4
st_insts = []
for i in [1 2 4 8 16]
    global st_insts
    coefs = round.([1/j for j = 1:(i+1)],digits=4)
    star_stencil = def_stencil_expression(:(c[0]D[x,y,z] + @sum(i, 1,$i, c[i]*(
                D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
                D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))))

    st_def = CreateStencilDefinition(star_stencil, coefs)
    bdim = 32
    if i > 8
        bdim = 16
    end
    st_inst = NewStencilInstance(st_def, m_step=false, bdim=bdim)
    push!(st_insts, st_inst)

end

## Warm Up
function warmup(st_insts)
    nx = 6
    ny = 6
    nz = 6
    dx = 1<<(nx)
    dy = 1<<(ny)
    dz = 1<<(nz)

    data = CreateData(dx,dy,dz)

    t_steps = 1

    gpu_out = ApplyFFTstencil(st_insts[1], data, t_steps)
end
warmup(st_insts)

## Benchmark
function bench(st_insts)
    nx = 8
    ny = 8
    nz = 8
    dx = 1<<(nx)
    dy = 1<<(ny)
    dz = 1<<(nz)

    data = CreateData(dx,dy,dz)

    t_steps = 1

    CUDA.cuProfilerStart()
    for i in st_insts
        ##global data, t_steps
        NVTX.@range "r$(i.max_radius)" begin
            gpu_out = ApplyFFTstencil(i, data, t_steps)
        end
    end
end
bench(st_insts)
exit()


text:
"""
find upper and lower bounds on range_theta based on MLE estimates for extremal range_lambda values

INPUTS:

"""
function mle(btg, glz, range_lambda, cholesky)
    N =

    g = btg.g
    lmin = range_lambda[1]
    lmax = range_lambda[2]

    glzs1 = copy(glz)
    for i = 1:length(glz)
        glzs1[i] = g.(glz[i], lmin)
    end
    glzs2 = copy(glz)
    for i = 1:length(glz)
        glzs2[i] = g.(glz[i], lmax)
    end

    function choltheta

    L1 = 0.5*logdet(cholesky) + 0.5 * glzs1

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Autogenerated wrapper script for Nettle_jll for x86_64-apple-darwin
export libhogweed, libnettle, nettle_hash

using GMP_jll
JLLWrappers.@generate_wrapper_header("Nettle")
JLLWrappers.@declare_library_product(libhogweed, "@rpath/libhogweed.6.dylib")
JLLWrappers.@declare_library_product(libnettle, "@rpath/libnettle.8.dylib")
JLLWrappers.@declare_executable_product(nettle_hash)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libhogweed,
        "lib/libhogweed.6.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnettle,
        "lib/libnettle.8.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nettle_hash,
        "bin/nettle-hash",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# https://github.com/FluxML/Flux.jl/pull/524

Base.prod(xs::TrackedArray; dims=:) = track(prod, xs; dims=dims)

@grad prod(xs; dims=:) = _prod(xs.data, prod(xs.data, dims=dims), dims)
_prod(xd, p, ::Colon) = p, Δ -> (nobacksies(:prod, ∇prod(xd, p, data(Δ)) ),)
_prod(xd, p, dims) = count(iszero, p) == 0 ?
  (p, Δ -> (nobacksies(:prod, p ./ xd .* data(Δ) ),)) :
  (p, Δ -> (nobacksies(:prod, mapslices(∇prod, xd; dims=dims) .* data(Δ)),))

function ∇prod(x, p=prod(x), Δ=1)
  numzero = count(iszero, x)
  if numzero == 0
    ∇ = p ./ x .* Δ
  elseif numzero > 1
    ∇ = zero(x)
  else
    ∇ = ∇prod_one(x, Δ)
  end
end
function ∇prod_one(x::Array, Δ)
  zloc = findfirst(iszero, x)
  ∇ = copy(x)
  ∇[zloc] = 1
  nonzero = prod(∇) * Δ
  ∇ .= 0
  ∇[zloc] = nonzero
  ∇
end
∇prod_one(x::AbstractArray, Δ) = ForwardDiff.gradient(y -> prod(y) * Δ, x)

Base.cumsum(xs::TrackedArray; dims=1) = track(cumsum, xs; dims=dims)

@grad cumsum(xs; dims=1) = _cumsum(xs.data, dims)
_cumsum(xd::Array, d) = cumsum(xd; dims=d), Δ -> (reverse(cumsum(reverse(Δ,dims=d),dims=d),dims=d),)
_cumsum(xd::AbstractArray, d) = cumsum(xd; dims=d), Δ -> (mapslices(reverse∘cumsum∘reverse,Δ, dims=d),)

Base.cumprod(xs::TrackedArray; dims=nothing) = track(cumprod, xs; dims=dims)

@grad cumprod(xs; dims=nothing) = _cumprod(xs.data, dims)
_cumprod(xd, ::Nothing, p = cumprod(xd)) = p, Δ -> (nobacksies(:cumprod, ∇cumprod(xd, p, data(Δ)) ),)
_cumprod(xd, d, p = cumprod(xd, dims=d)) = p, Δ -> (nobacksies(:cumprod, ∇cumprod_d(xd, Val(d), p, data(Δ)) ),)

function ∇cumprod(x::Vector, p, Δ)
  len = length(x)
  z = something(findfirst(iszero, x), len+1)
  ∇ = zero(x)
  @inbounds for i=1:z-1
    ixi = 1/x[i]
    for k=i:z-1
      ∇[i] += p[k] * Δ[k] * ixi
    end
  end
  @inbounds if z != len+1
    pk = z==1 ? one(p[1]) : p[z-1] # will be prod(x[j] for j=1:k if j!=z)
    ∇[z] += pk * Δ[z]
    for k=(z+1):len
      pk *= x[k]
      ∇[z] += pk * Δ[k]
    end
  end
  ∇
end
∇cumprod(x::AbstractVector, p, Δ) = vec(Δ' * ForwardDiff.jacobian(cumprod, x))
@noinline function ∇cumprod_d(x::AbstractArray{T,N}, ::Val{d}, p, Δ) where {T,N,d}
  ∇ = similar(x)
  for i in Iterators.product(ntuple(k -> k==d ? Ref(:) : axes(x,k), Val(N))...)
    copyto!(view(∇,i...), ∇cumprod(x[i...], p[i...], Δ[i...]))
  end
  ∇ # roughly mapslices(∇cumprod, x,p,Δ; dims=d) if that existed
end


text:
# Julia wrapper for header: /usr/include/libavcodec/jni.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_jni_set_java_vm,
    av_jni_get_java_vm

function av_jni_set_java_vm(vm, log_ctx)
    ccall((:av_jni_set_java_vm, libavcodec), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), vm, log_ctx)
end

function av_jni_get_java_vm(log_ctx)
    ccall((:av_jni_get_java_vm, libavcodec), Ptr{Cvoid}, (Ptr{Cvoid},), log_ctx)
end


text:
### Calculate the number of passengers on each flight
reload("nmr")

nmr.runjob(nmr.NMR(3, ["AComp_Passenger_data.csv"], nmr.mapper_parserecordacomp, "acomp.csv"))
nmr.runjob(nmr.NMR(3, ["acomp.csv"], nmr.reducer_numberofpassengers, "acomp_reduced.json"))
c = nmr.runcombiner("acomp_reduced.json", nmr.combiner_parsejson)

@test typeof(c) == Vector{String}
@test length(c) == 17

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
let
    n = mean_motion(R_EARTH + 500e3, use_degrees=false)
    @test isapprox(n, 0.0011067836148773837, atol=1e-12)

    n = mean_motion(R_EARTH + 500e3, use_degrees=true)
    @test isapprox(n, 0.0634140299667068, atol=1e-12)
end

let
    a = semimajor_axis(0.0011067836148773837, use_degrees=false)
    @test isapprox(a, R_EARTH + 500e3, atol=1e-6)

    a = semimajor_axis(0.0634140299667068, use_degrees=true)
    @test isapprox(a, R_EARTH + 500e3, atol=1e-6)
end

let
    T = orbit_period(R_EARTH + 500e3)
    @test isapprox(T, 5676.977164028288, atol=1e-9)
end

let
    iss = sun_sync_inclination(R_EARTH + 574e3, 0.0, use_degrees=false)
    @test isapprox(iss, 97.685*pi/180, atol=1.0e-3)

    iss = sun_sync_inclination(R_EARTH + 574e3, 0.0, use_degrees=true)
    @test isapprox(iss, 97.685, atol=1.0e-3)
end

let
    # 0
    M = anomaly_eccentric_to_mean(0.0, 0.0, use_degrees=false)
    @test M == 0

    M = anomaly_eccentric_to_mean(0.0, 0.0, use_degrees=true)
    @test M == 0

    # 180
    M = anomaly_eccentric_to_mean(pi/2, 0.1, use_degrees=false)
    @test isapprox(M, 1.4707, atol=1e-3)

    M = anomaly_eccentric_to_mean(90.0, 0.1, use_degrees=true)
    @test isapprox(M, 84.270, atol=1e-3)

    # 180
    M = anomaly_eccentric_to_mean(pi, 0.0, use_degrees=false)
    @test isapprox(M, pi, atol=1e-12)

    M = anomaly_eccentric_to_mean(180.0, 0.0, use_degrees=true)
    @test M == 180.0
end

let
    # 0
    E = anomaly_mean_to_eccentric(0.0, 0.0, use_degrees=false)
    @test E == 0

    E = anomaly_mean_to_eccentric(0.0, 0.0, use_degrees=true)
    @test E == 0

    # 180
    E = anomaly_mean_to_eccentric(1.4707, 0.1, use_degrees=false)
    @test isapprox(E, pi/2, atol=1e-3)

    E = anomaly_mean_to_eccentric(84.270, 0.1, use_degrees=true)
    @test isapprox(E, 90.0, atol=1e-3)

    # 180
    E = anomaly_mean_to_eccentric(pi, 0.0, use_degrees=false)
    @test isapprox(E, pi, atol=1e-12)

    E = anomaly_mean_to_eccentric(180.0, 0.0, use_degrees=true)
    @test E == 180.0

    # Large Eccentricities
    E = anomaly_mean_to_eccentric(180.0, 0.9, use_degrees=true)
    @test E == 180.0

    # Max Iteration exit on hyperbolic trajectory
    # @test_throws ErrorException anomaly_mean_to_eccentric(pi, 2.0)
end

let
    oe  = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    eci = sOSCtoCART(oe, use_degrees=true)

    tol = 1e-6
    @test isapprox(eci[1], R_EARTH + 500e3, atol=tol)
    @test isapprox(eci[2], 0.0, atol=tol)
    @test isapprox(eci[3], 0.0, atol=tol)
    @test isapprox(eci[4], 0.0, atol=tol)
    @test isapprox(eci[5], 0.0, atol=tol)
    @test isapprox(eci[6], sqrt(GM_EARTH/(R_EARTH + 500e3)), atol=tol)
end

let
    # Using radians
    oe   = [R_EARTH + 500e3, 0, pi/2.0, 0, 0, 0]
    eci  = sOSCtoCART(oe, use_degrees=false)
    eci2 = sOSCtoCART(oe, use_degrees=false)

    tol = 1e-6
    @test isapprox(eci[1], eci2[1], atol=tol)
    @test isapprox(eci[2], eci2[2], atol=tol)
    @test isapprox(eci[3], eci2[3], atol=tol)
    @test isapprox(eci[4], eci2[4], atol=tol)
    @test isapprox(eci[5], eci2[5], atol=tol)
    @test isapprox(eci[6], eci2[6], atol=tol)

    # Using degrees
    oe   = [R_EARTH + 500e3, 0, 90.0, 0, 0, 0]
    eci  = sOSCtoCART(oe, use_degrees=true)
    eci2 = sOSCtoCART(oe, use_degrees=true)

    tol = 1e-6
    @test isapprox(eci[1], eci2[1], atol=tol)
    @test isapprox(eci[2], eci2[2], atol=tol)
    @test isapprox(eci[3], eci2[3], atol=tol)
    @test isapprox(eci[4], eci2[4], atol=tol)
    @test isapprox(eci[5], eci2[5], atol=tol)
    @test isapprox(eci[6], eci2[6], atol=tol)
end

let
    eci   = [R_EARTH + 500e3, 100e3, 575e3, 0, 0, 7300]
    eci2  = sOSCtoCART(sCARTtoOSC(eci))

    tol = 1e-6
    @test isapprox(eci[1], eci2[1], atol=tol)
    @test isapprox(eci[2], eci2[2], atol=tol)
    @test isapprox(eci[3], eci2[3], atol=tol)
    @test isapprox(eci[4], eci2[4], atol=tol)
    @test isapprox(eci[5], eci2[5], atol=tol)
    @test isapprox(eci[6], eci2[6], atol=tol)

    # Equatorial circulator
    a   = R_EARTH + 1000e3
    e   = 0.0
    eci = [a, 0, 0, 0, 0, sqrt(GM_EARTH/a)]
    oe  = sCARTtoOSC(eci, use_degrees=true)

    tol = 1e-6
    @test isapprox(oe[1], a,     atol=tol)
    @test isapprox(oe[2], e,     atol=tol)
    @test isapprox(oe[3], 90.0,  atol=tol)
end

let
    # Test near-circular conversions
    a   = R_EARTH + 500e3
    eci = [a, 0.0, 0.0, 0.0, 0.0, sqrt(GM_EARTH/a)]
    oe  = sCARTtoOSC(eci, use_degrees=true)

    tol = 1e-6
    @test isapprox(oe[1], a, atol=tol)
    @test isapprox(oe[2], 0.0, atol=tol)
    @test isapprox(oe[3], 90.0, atol=tol)
    @test isapprox(oe[4], 0.0, atol=tol)
    @test isapprox(oe[5], 0.0, atol=tol)
    @test isapprox(oe[6], 0.0, atol=tol)
end


text:
######
###### `inv`
######

function frule((_, Δx), ::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    return Ω, -Ω * Δx * Ω
end

function rrule(::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    function inv_pullback(ΔΩ)
        return NO_FIELDS, -Ω' * ΔΩ * Ω'
    end
    return Ω, inv_pullback
end

#####
##### `*`
#####

function rrule(
    ::typeof(*),
    A::AbstractVecOrMat{<:CommutativeMulNumber},
    B::AbstractVecOrMat{<:CommutativeMulNumber},
)
    function times_pullback(Ȳ)
        return (
            NO_FIELDS,
            InplaceableThunk(
                @thunk(Ȳ * B'),
                X̄ -> mul!(X̄, Ȳ, B', true, true)
            ),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, A', Ȳ, true, true)
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
    ::typeof(*),
    A::AbstractVector{<:CommutativeMulNumber},
    B::AbstractMatrix{<:CommutativeMulNumber},
)
    function times_pullback(Ȳ)
        @assert size(B, 1) === 1   # otherwise primal would have failed.
        return (
            NO_FIELDS,
            InplaceableThunk(
                @thunk(Ȳ * vec(B')),
                X̄ -> mul!(X̄, Ȳ, vec(B'), true, true)
            ),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, A', Ȳ, true, true)
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
   ::typeof(*), A::CommutativeMulNumber, B::AbstractArray{<:CommutativeMulNumber}
)
    function times_pullback(Ȳ)
        return (
            NO_FIELDS,
            @thunk(dot(Ȳ, B)'),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, conj(A), Ȳ, true, true)
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
    ::typeof(*), B::AbstractArray{<:CommutativeMulNumber}, A::CommutativeMulNumber
)
    function times_pullback(Ȳ)
        return (
            NO_FIELDS,
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, conj(A), Ȳ, true, true)
            ),
            @thunk(dot(Ȳ, B)'),
        )
    end
    return A * B, times_pullback
end

#####
##### `/`
#####

function rrule(::typeof(/), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Aᵀ, dA_pb = rrule(adjoint, A)
    Bᵀ, dB_pb = rrule(adjoint, B)
    Cᵀ, dS_pb = rrule(\, Bᵀ, Aᵀ)
    C, dC_pb = rrule(adjoint, Cᵀ)
    function slash_pullback(Ȳ)
        # Optimization note: dAᵀ, dBᵀ, dC are calculated no matter which partial you want
        _, dC = dC_pb(Ȳ)
        _, dBᵀ, dAᵀ = dS_pb(unthunk(dC))

        ∂A = last(dA_pb(unthunk(dAᵀ)))
        ∂B = last(dA_pb(unthunk(dBᵀ)))

        (NO_FIELDS, ∂A, ∂B)
    end
    return C, slash_pullback
end

#####
##### `\`
#####

function rrule(::typeof(\), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Y = A \ B
    function backslash_pullback(Ȳ)
        ∂A = @thunk begin
            B̄ = A' \ Ȳ
            Ā = -B̄ * Y'
            Ā = add!!(Ā, (B - A * Y) * B̄' / A')
            Ā = add!!(Ā, A' \ Y * (Ȳ' - B̄'A))
            Ā
        end
        ∂B = @thunk A' \ Ȳ
        return NO_FIELDS, ∂A, ∂B
    end
    return Y, backslash_pullback

end

#####
##### `\`, `/` matrix-scalar_rule
#####

function rrule(::typeof(/), A::AbstractArray{<:Real}, b::Real)
    Y = A/b
    function slash_pullback(Ȳ)
        return (NO_FIELDS, @thunk(Ȳ/b), @thunk(-dot(Ȳ, Y)/b))
    end
    return Y, slash_pullback
end

function rrule(::typeof(\), b::Real, A::AbstractArray{<:Real})
    Y = b\A
    function backslash_pullback(Ȳ)
        return (NO_FIELDS, @thunk(-dot(Ȳ, Y)/b), @thunk(Ȳ/b))
    end
    return Y, backslash_pullback
end

#####
##### Negation (Unary -)
#####

function rrule(::typeof(-), x::AbstractArray)
    function negation_pullback(ȳ)
        return NO_FIELDS, InplaceableThunk(@thunk(-ȳ), ā -> ā .-= ȳ)
    end
    return -x, negation_pullback
end


text:
# Interface definition for Soapy SDR devices.

# General design rules about the API:
# The caller must free non-const array results.

# Forward declaration of device handle
mutable struct SoapySDRDevice
end

# Forward declaration of stream handle
mutable struct SoapySDRStream
end

# Get the last status code after a Device API call.
# The status code is cleared on entry to each Device call.
# When an device API call throws, the C bindings catch
# the exception, and set a non-zero last status code.
# Use lastStatus() to determine success/failure for
# Device calls without integer status return codes.
function SoapySDRDevice_lastStatus()
    ccall((:SoapySDRDevice_lastStatus, lib), Cint, ())
end

# Get the last error message after a device call fails.
# When an device API call throws, the C bindings catch
# the exception, store its message in thread-safe storage,
# and return a non-zero status code to indicate failure.
# Use lastError() to access the exception's error message.
function SoapySDRDevice_lastError()
    ccall((:SoapySDRDevice_lastError, lib), Cstring, ())
end

# Enumerate a list of available devices on the system.
# param args device construction key/value argument filters
# param [out] length the number of elements in the result.
# return a list of arguments strings, each unique to a device
function SoapySDRDevice_enumerate()
    sz = Ref{Csize_t}()
    kwargs = ccall((:SoapySDRDevice_enumerate, lib), Ptr{SoapySDRKwargs}, (Ptr{Nothing}, Ref{Csize_t}), C_NULL, sz)
    return (kwargs, sz)
end

# Enumerate a list of available devices on the system.
# Markup format for args: "keyA=valA, keyB=valB".
# param args a markup string of key/value argument filters
# param [out] length the number of elements in the result.
# return a list of arguments strings, each unique to a device
function SoapySDRDevice_enumerateStrArgs(args)
    sz = Ref{Csize_t}()
    kwargs = ccall((:SoapySDRDevice_enumerateStrArgs, lib), Ptr{SoapySDRKwargs}, (Cstring, Ref{Csize_t}), args, sz)
    return (kwargs, sz)
end

# Make a new Device object given device construction args.
# The device pointer will be stored in a table so subsequent calls
# with the same arguments will produce the same device.
# For every call to make, there should be a matched call to unmake.
# param args device construction key/value argument map
# return a pointer to a new Device object
function SoapySDRDevice_make(args) # have not checked
    r = ccall((:SoapySDRDevice_make, lib), Ptr{SoapySDRDevice}, (Ref{SoapySDRKwargs},), args)
end

# Make a new Device object given device construction args.
# The device pointer will be stored in a table so subsequent calls
# with the same arguments will produce the same device.
# For every call to make, there should be a matched call to unmake.
#
# param args a markup string of key/value arguments
# return a pointer to a new Device object or null for error
function SoapySDRDevice_makeStrArgs(args) # have not checked
    ccall((:SoapySDRDevice_makeStrArgs, lib), Ptr{SoapySDRDevice}, (Cstring,), args)
end

# Unmake or release a device object handle.
# param device a pointer to a device object
# return 0 for success or error code on failure
function SoapySDRDevice_unmake(device) # have not checked
    ccall((:SoapySDRDevice_unmake, lib), Cint, (Ptr{SoapySDRDevice},), device)
end

######################
## PARALLEL SUPPORT ##      # have not checked any of these
######################

# Create a list of devices from a list of construction arguments.
# This is a convenience call to parallelize device construction,
# and is fundamentally a parallel for loop of make(Kwargs).
#
# \param argsList a list of device arguments per each device
# \param length the length of the argsList array
# \return a list of device pointers per each specified argument
function SoapySDRDevice_make_list(argsList, length::Cint)
    ccall((:SoapySDRDevice_make_list, lib), Ptr{Ptr{SoapySDRDevice}}, (Ptr{Cint}, Cint), argsList, length)
end

# Unmake or release a list of device handles
# and free the devices array memory as well.
# This is a convenience call to parallelize device destruction,
# and is fundamentally a parallel for loop of unmake(Device *).
# param devices a list of pointers to device objects
# param length the length of the devices array
# return 0 for success or error code on failure
function SoapySDRDevice_unmake_list(devices, length::Cint)
    ccall((:SoapySDRDevice_unmake_list, lib), Cint, (Ptr{Ptr{SoapySDRDevice}}, Cint), devices, length)
end

########################
## Identificatoin API ##    # have not checked any of these
########################

# A key that uniquely identifies the device driver.
# This key identifies the underlying implementation.
# Serveral variants of a product may share a driver.
# param device a pointer to a device instance
function SoapySDRDevice_getDriverKey(device)
    ccall((:SoapySDRDevice_getDriverKey, lib), Cstring, (Ptr{SoapySDRDevice},), device)
end

# A key that uniquely identifies the hardware.
# This key should be meaningful to the user
# to optimize for the underlying hardware.
# \param device a pointer to a device instance
function SoapySDRDevice_getHardwareKey(device)
    ccall((:SoapySDRDevice_getHardwareKey, lib), Cstring, (Ptr{SoapySDRDevice},), device)
end

# Query a dictionary of available device information.
# This dictionary can any number of values like
# vendor name, product name, revisions, serials...
# This information can be displayed to the user
# to help identify the instantiated device.
# \param device a pointer to a device instance
function SoapySDRDevice_getHardwareInfo()
    ccall((:SoapySDRDevice_getHardwareInfo, lib), Cint, ())
end

#################
## ANTENNA API ##
#################

# Get a list of available antennas to select on a given chain.
# param device a pointer to a device instance
# param direction the channel direction RX or TX
# param channel an available channel on the device
# param [out] length the number of antenna names
# return a list of available antenna names
function SoapySDRDevice_listAntennas(device, direction, channel)
    leng = Ref{Csize_t}()
    names = ccall((:SoapySDRDevice_listAntennas, lib), Ptr{Cstring}, (Ptr{SoapySDRDevice}, Cint, Csize_t, Ref{Csize_t}), device, direction, channel, leng)
    return (names, leng)
end

##############
## GAIN API ##
##############

# List available amplification elements.
# Elements should be in order RF to baseband.
# param device a pointer to a device instance
# param direction the channel direction RX or TX
# param channel an available channel
# param [out] length the number of gain names
# return a list of gain string names
#function SoapySDRDevice_listGains(device, direction, channel, length)
function SoapySDRDevice_listGains(device, direction, channel)
    leng = Ref{Csize_t}()
    names = ccall((:SoapySDRDevice_listGains, "libSoapySDR.so"), Ptr{Cstring}, (Ptr{SoapySDRDevice}, Cint, Csize_t, Ref{Csize_t}), device, direction, channel, leng)
    return (names, leng)
end

###################
## FREQUENCY API ##
###################

# Get the range of overall frequency values.
# param device a pointer to a device instance
# param direction the channel direction RX or TX
# param channel an available channel on the device
# param [out] length the number of ranges
# return a list of frequency ranges in Hz
function SoapySDRDevice_getFrequencyRange(device, direction, channel)
    leng = Ref{Csize_t}()
    ranges = ccall((:SoapySDRDevice_getFrequencyRange, lib), Ptr{SoapySDRRange}, (Ptr{SoapySDRDevice}, Cint, Csize_t, Ref{Csize_t}), device, direction, channel, leng)
    return (ranges, leng)
end

#####################
## SAMPLE RATE API ##
#####################

# Set the baseband sample rate of the chain.
# param device a pointer to a device instance
# param direction the channel direction RX or TX
# param channel an available channel on the device
# param rate the sample rate in samples per
function SoapySDRDevice_setSampleRate(device, direction, channel, rate)
    ccall((:SoapySDRDevice_setSampleRate, lib), Cint, (Ref{SoapySDRDevice}, Cint, Csize_t, Cdouble), device, direction, channel, rate)
end

###################
## FREQUENCY API ##
###################

# Set the center frequency of the chain.
#  - For RX, this specifies the down-conversion frequency.
#  - For TX, this specifies the up-conversion frequency.
# The default implementation of setFrequency() will tune the "RF"
# component as close as possible to the requested center frequency.
# Tuning inaccuracies will be compensated for with the "BB" component.
# The args can be used to augment the tuning algorithm.
#  - Use "OFFSET" to specify an "RF" tuning offset,
#    usually with the intention of moving the LO out of the passband.
#    The offset will be compensated for using the "BB" component.
#  - Use the name of a component for the key and a frequency in Hz
#    as the value (any format) to enforce a specific frequency.
#    The other components will be tuned with compensation
#    to achieve the specified overall frequency.
#  - Use the name of a component for the key and the value "IGNORE"
#    so that the tuning algorithm will avoid altering the component.
#  - Vendor specific implementations can also use the same args to augment
#    tuning in other ways such as specifying fractional vs integer N tuning.
# param device a pointer to a device instance
# param direction the channel direction RX or TX
# param channel an available channel on the device
# param frequency the center frequency in Hz
# param args optional tuner arguments
# return an error code or 0 for success
function SoapySDRDevice_setFrequency(device, direction, channel, frequency)
    ccall((:SoapySDRDevice_setFrequency, lib), Int, (Ref{SoapySDRDevice}, Cint, Csize_t, Cdouble, Ptr{Nothing}), device, direction, channel, frequency, C_NULL)
end

################
## STREAM API ##
################

# Initialize a stream given a list of channels and stream arguments.
# The implementation may change switches or power-up components.
# All stream API calls should be usable with the new stream object
# after setupStream() is complete, regardless of the activity state.
#
# The API allows any number of simultaneous TX and RX streams, but many dual-channel
# devices are limited to one stream in each direction, using either one or both channels.
# This call will return an error if an unsupported combination is requested,
# or if a requested channel in this direction is already in use by another stream.
#
# When multiple channels are added to a stream, they are typically expected to have
# the same sample rate. See SoapySDRDevice_setSampleRate().
#
# param device a pointer to a device instance
# param [out] stream the opaque pointer to a stream handle.
#
# The returned stream is not required to have internal locking, and may not be used
# concurrently from multiple threads.
#
# param direction the channel direction (`SOAPY_SDR_RX` or `SOAPY_SDR_TX`)
# param format A string representing the desired buffer format in read/writeStream()
# parblock
#
# The first character selects the number type:
#   - "C" means complex
#   - "F" means floating point
#   - "S" means signed integer
#   - "U" means unsigned integer
#
# The type character is followed by the number of bits per number (complex is 2x this size per sample)
#
#  Example format strings:
#   - "CF32" -  complex float32 (8 bytes per element)
#   - "CS16" -  complex int16 (4 bytes per element)
#   - "CS12" -  complex int12 (3 bytes per element)
#   - "CS4" -  complex int4 (1 byte per element)
#   - "S32" -  int32 (4 bytes per element)
#   - "U8" -  uint8 (1 byte per element)
#
# endparblock
# param channels a list of channels or empty for automatic
# param numChans the number of elements in the channels array
# param args stream args or empty for defaults
# parblock
#
#   Recommended keys to use in the args dictionary:
#    - "WIRE" - format of the samples between device and host
# endparblock
# return 0 for success or error code on failure
function SoapySDRDevice_setupStream(device, stream, direction, format, channels, numChans)
    ccall((:SoapySDRDevice_setupStream, lib), Int, (Ref{SoapySDRDevice}, Ref{Ref{SoapySDRStream}}, Cint, Cstring, Ref{Csize_t}, Csize_t, Ptr{Nothing}), device, stream, direction, format, channels, numChans, C_NULL)
end

# Activate a stream.
# Call activate to prepare a stream before using read/write().
# The implementation control switches or stimulate data flow.
#
# The timeNs is only valid when the flags have SOAPY_SDR_HAS_TIME.
# The numElems count can be used to request a finite burst size.
# The SOAPY_SDR_END_BURST flag can signal end on the finite burst.
# Not all implementations will support the full range of options.
# In this case, the implementation returns SOAPY_SDR_NOT_SUPPORTED.
#
# param device a pointer to a device instance
# param stream the opaque pointer to a stream handle
# param flags optional flag indicators about the stream
# param timeNs optional activation time in nanoseconds
# param numElems optional element count for burst control
# return 0 for success or error code on failure
function SoapySDRDevice_activateStream(device, stream, flags, timeNs, numElems)
    ccall((:SoapySDRDevice_activateStream, lib), Cint, (Ref{SoapySDRDevice}, Ref{SoapySDRStream}, Cint, Clonglong, Csize_t), device, stream, flags, timeNs, numElems)
end

# Read elements from a stream for reception.
# This is a multi-channel call, and buffs should be an array of void *,
# where each pointer will be filled with data from a different channel.
#
# **Client code compatibility:**
# The readStream() call should be well defined at all times,
# including prior to activation and after deactivation.
# When inactive, readStream() should implement the timeout
# specified by the caller and return SOAPY_SDR_TIMEOUT.
#
# param device a pointer to a device instance
# param stream the opaque pointer to a stream handle
# param buffs an array of void* buffers num chans in size
# param numElems the number of elements in each buffer
# param [out] flags optional flag indicators about the result
# param [out] timeNs the buffer's timestamp in nanoseconds
# param timeoutUs the timeout in microseconds
# return the number of elements read per buffer or error code
function SoapySDRDevice_readStream(device, stream, buffs, numElems, flags, timeNs, timeoutUs)
    ccall((:SoapySDRDevice_readStream, lib), Cint, (Ref{SoapySDRDevice}, Ref{SoapySDRStream}, Ptr{Ptr{Cvoid}}, Csize_t, Ptr{Cint}, Ptr{Clonglong}, Clong), device, stream, map(pointer, buffs), numElems, flags, timeNs, timeoutUs)
end

# Deactivate a stream.
# Call deactivate when not using using read/write().
# The implementation control switches or halt data flow.
#
# The timeNs is only valid when the flags have SOAPY_SDR_HAS_TIME.
# Not all implementations will support the full range of options.
# In this case, the implementation returns SOAPY_SDR_NOT_SUPPORTED.
#
# param device a pointer to a device instance
# param stream the opaque pointer to a stream handle
# param flags optional flag indicators about the stream
# param timeNs optional deactivation time in nanoseconds
# return 0 for success or error code on failure
function SoapySDRDevice_deactivateStream(device, stream, flags, timeNs)
    ccall((:SoapySDRDevice_deactivateStream, lib), Cint, (Ref{SoapySDRDevice}, Ref{SoapySDRStream}, Cint, Clonglong), device, stream, flags, timeNs)
end

# Close an open stream created by setupStream
# param device a pointer to a device instance
# param stream the opaque pointer to a stream handle
# return 0 for success or error code on failure
function SoapySDRDevice_closeStream(device, stream)
    ccall((:SoapySDRDevice_closeStream, lib), Cint, (Ref{SoapySDRDevice}, Ref{SoapySDRStream}), device, stream)
end


text:
export Coulomb

mutable struct Coulomb <: MoleculeDescribe
    number::Int64
    value::Array{Float64, 1}

    Coulomb() = new(0, [])
    function Coulomb(molecule::Atoms)
        coulombMatrix = Array{Float64, 2}(undef, Atoms.number, Atoms.number)
        for i in 1:(Atoms.number-1)
            chargeI = getCharge(Atoms.atom[i].symbol)
            coulombMatrix[i, i] = chargeI^2.4

            for j in (i+1):(Atoms.number-1)
                chargeJ = getCharge(Atoms.atom[j].symbol)
                distance = dist(Atoms.atom[i], Atoms.atom[j])
                coulombMatrix[i, j] = chargeI * chargeJ / distance
                coulombMatrix[j, i] = coulombMatrix[i, j]
            end
        end

        eigenValues = eigvals(coulombMatrix)
        sort!(eigenValues, rev=true)

        new(Atoms.number, eigenValues)
    end
end


text:
const simplegaussstandata = [
  Dict(
  "N" => 2,
  "x" => [1.5, 2.0]
  )
]


text:
"""
Continuous SIMP algorithm, see [TarekRay2020](@cite).
"""
@params mutable struct ContinuationSIMP <: AbstractSIMP
    simp::Any
    result::Any
    options::Any
    callback::Any
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::ContinuationSIMP)
    return println("TopOpt continuation SIMP algorithm")
end

@params struct CSIMPOptions
    p_cont::Any
    option_cont::Any
    reuse::Bool
    steps::Integer
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::CSIMPOptions)
    return println("TopOpt continuation SIMP options")
end

function CSIMPOptions(
    ::Type{T}=Float64;
    steps=40,
    p_gen=nothing,
    initial_options=MMAOptions(),
    pstart=T(1),
    pfinish=T(5),
    reuse=false,
    options_gen=nothing,
) where {T}
    if p_gen == nothing
        p_cont = PowerContinuation{T}(;
            b=T(1), start=pstart, steps=steps + 1, finish=pfinish
        )
    else
        @assert steps == p_gen.length - 1
        p_cont = p_gen
    end

    if options_gen == nothing
        options_cont = MMAOptionsGen(; steps=steps, initial_options=initial_options)
    else
        options_cont = options_gen
    end

    return CSIMPOptions(p_cont, options_cont, reuse, steps)
end

function ContinuationSIMP(
    simp::SIMP{T},
    steps::Int=40,
    options::CSIMPOptions=CSIMPOptions(
        T; steps=steps, initial_options=deepcopy(simp.optimizer.options)
    ),
    callback=(i) -> (),
) where {T}
    @assert steps == options.steps
    ncells = getncells(simp.solver.problem)
    result = NewSIMPResult(T, simp.optimizer, ncells)
    return ContinuationSIMP(simp, result, options, callback)
end

# Unused functions
default_ftol(solver) = xmin / 10
function default_grtol(solver)
    @unpack xmin, problem = solver
    return xmin / 1000 * sum(problem.elementinfo.cellvolumes) * YoungsModulus(problem)
end

function MMAOptionsGen(;
    steps::Int=40,
    initial_options=MMAOptions(),
    ftol_gen=nothing,
    xtol_gen=nothing,
    grtol_gen=nothing,
    kkttol_gen=nothing,
    maxiter_gen=nothing,
    outer_maxiter_gen=nothing,
    s_init_gen=nothing,
    s_incr_gen=nothing,
    s_decr_gen=nothing,
    store_trace_gen=nothing,
    show_trace_gen=nothing,
    auto_scale_gen=nothing,
    dual_options_gen=nothing,
)
    if maxiter_gen == nothing
        maxiter_cont = FixedContinuation(initial_options.maxiter, steps + 1)
    else
        @assert steps == maxiter_gen.length - 1
        maxiter_cont = maxiter_gen
    end

    if outer_maxiter_gen == nothing
        outer_maxiter_cont = FixedContinuation(initial_options.outer_maxiter, steps + 1)
    else
        @assert steps == outer_maxiter_gen.length - 1
        outer_maxiter_cont = outer_maxiter_gen
    end

    if xtol_gen == nothing
        xtol_cont = FixedContinuation(initial_options.tol.x, steps + 1)
    else
        @assert steps == xtol_gen.length - 1
        xtol_cont = xtol_gen
    end

    if ftol_gen == nothing
        ftol_cont = FixedContinuation(initial_options.tol.frel, steps + 1)
    else
        @assert steps == ftol_gen.length - 1
        ftol_cont = ftol_gen
    end

    if kkttol_gen == nothing
        kkttol_cont = FixedContinuation(initial_options.tol.kkt, steps + 1)
    else
        @assert steps == kkttol_gen.length - 1
        kkttol_cont = kkttol_gen
    end

    tol_cont = Nonconvex.Tolerance(; x=xtol_cont, f=ftol_cont, kkt=kkttol_cont)

    if s_init_gen == nothing
        s_init_cont = FixedContinuation(initial_options.s_init, steps + 1)
    else
        @assert steps == s_init_gen.length - 1
        s_init_cont = s_init_gen
    end

    if s_incr_gen == nothing
        s_incr_cont = FixedContinuation(initial_options.s_incr, steps + 1)
    else
        @assert steps == s_incr_gen.length - 1
        s_incr_cont = s_incr_gen
    end
    if s_decr_gen == nothing
        s_decr_cont = FixedContinuation(initial_options.s_decr, steps + 1)
    else
        @assert steps == s_decr_gen.length - 1
        s_decr_cont = s_decr_gen
    end

    if store_trace_gen == nothing
        store_trace_cont = FixedContinuation(initial_options.store_trace, steps + 1)
    else
        @assert steps == store_trace_gen.length - 1
        store_trace_cont = store_trace_gen
    end

    if show_trace_gen == nothing
        show_trace_cont = FixedContinuation(initial_options.show_trace, steps + 1)
    else
        @assert steps == show_trace_gen.length - 1
        show_trace_cont = show_trace_gen
    end

    if auto_scale_gen == nothing
        auto_scale_cont = FixedContinuation(initial_options.auto_scale, steps + 1)
    else
        @assert steps == auto_scale_gen.length - 1
        auto_scale_cont = auto_scale_gen
    end

    if dual_options_gen == nothing
        dual_options_cont = FixedContinuation(initial_options.dual_options, steps + 1)
    else
        @assert steps == dual_options_gen.length - 1
        dual_options_cont = dual_options_gen
    end

    return MMAOptionsGen(
        maxiter_cont,
        outer_maxiter_cont,
        tol_cont,
        s_init_cont,
        s_incr_cont,
        s_decr_cont,
        store_trace_cont,
        show_trace_cont,
        auto_scale_cont,
        dual_options_cont,
    )
end

function (c_simp::ContinuationSIMP)(x0=copy(c_simp.simp.solver.vars))
    @unpack optimizer = c_simp.simp
    @unpack workspace = optimizer
    @unpack options = workspace

    update!(c_simp, 1)
    # Does the first function evaluation
    Nonconvex.NonconvexCore.reset!(workspace, x0)
    c_simp.callback(0)
    r = c_simp.simp(workspace)

    f_x_previous = workspace.solution.prevf
    g_previous = copy(c_simp.simp.optimizer.workspace.solution.g)

    original_maxiter = options.maxiter
    for i in 1:(c_simp.options.steps)
        c_simp.callback(i)
        update!(c_simp, i + 1)
        Nonconvex.NonconvexCore.reset!(workspace)

        maxiter = 1
        workspace.options.maxiter = maxiter
        r = c_simp.simp(workspace)

        if !workspace.solution.convstate.converged
            while !workspace.solution.convstate.converged
                maxiter += 1
                workspace.options.maxiter = maxiter
                c_simp.simp(workspace)
            end
        end
        f_x_previous = workspace.solution.prevf
        g_previous = copy(c_simp.simp.optimizer.workspace.solution.g)

        workspace.options.maxiter = original_maxiter
    end
    return c_simp.result = c_simp.simp.result
end

function update!(c_simp::ContinuationSIMP, i)
    p = c_simp.options.p_cont(i)
    setpenalty!(c_simp.simp, p)
    options = c_simp.options.option_cont(c_simp.simp.optimizer.workspace.options, i)
    c_simp.simp.optimizer.workspace.options = options
    return c_simp
end
frac(x) = 2 * min(abs(x), abs(x - 1))

function undo_values!(workspace, f_x_previous, g_previous)
    workspace.solution.f = workspace.solution.prevf
    workspace.solution.prevf = f_x_previous
    workspace.solution.g .= g_previous
    workspace.solution.x .= workspace.solution.prevx
    workspace.solution.prevx .= workspace.tempx
    return workspace
end


text:
"""
    Constructor for the Bayesian Parameters (weights on private signal and prior).
        Tested
"""
function BPars(n, par::ParsGrowth)
    @assert n >= 0
    @unpack_ParsGrowth par
    V_n = 1 / (σ²_θ^-1 + n * σ²_ϵ^-1)
    BParsGrowth(n = n,
        V_n = V_n,
        w_pn = σ²_θ^-1 * V_n,
        w_an = n * σ²_ϵ^-1 * V_n)
end

BPars(n; par::ParsGrowth) = BParsGrowth(n, par)

function BPars_dep(n, par::ParsGrowth)
    @assert n >= 0
    @unpack_ParsGrowth par
    V_n = par.uncertainty_grid[n+1]
    BParsGrowth(n = n,
        V_n = V_n,
        w_pn = σ²_θ^-1 * V_n,
        w_an = n * σ²_ϵ^-1 * V_n)
end

BPars_dep(n; par::ParsGrowth) = BPars_dep(n, par)

"""
    The expected mean given n signals with mean a_bar.
"""
μ_n(a_bar, par::ParsGrowth, bpar::BParsGrowth) = bpar.w_pn * par.θ_bar + bpar.w_an * a_bar
μ_n(a_bar; par::ParsGrowth, bpar::BParsGrowth) = μ_n(a_bar, par, bpar)
# Now lets write a function to compute the weights on each node

function Transition_Probability(par::ParsGrowth, bpar::BParsGrowth; normalize::Symbol = :yes)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    Π = zeros(N, N)
    μ_ns = μ_n.(a_grid, par = par, bpar = bpar)

    # Rows: Today's state i
    # Columns: Tomorrow's state j
    for i = 1:N, j = 1:N
        if j == 1
            d = (n + 1) * (a_grid[j+1] - a_grid[j])
            Π[i, j] = normcdf(((n + 1) * a_grid[j] - n * a_grid[i] - μ_ns[i] + d / 2) / √(V_n + σ²_ϵ))
        elseif j == N
            d = (n + 1) * (a_grid[j] - a_grid[j-1])
            Π[i, j] = 1 - normcdf(((n + 1) * a_grid[j] - n * a_grid[i] - μ_ns[i] - d / 2) / √(V_n + σ²_ϵ))
        else
            d_up = (n + 1) * (a_grid[j+1] - a_grid[j])
            d_down = (n + 1) * (a_grid[j] - a_grid[j-1])
            Π[i, j] = normcdf(((n + 1) * a_grid[j] - n * a_grid[i] - μ_ns[i] + d_up / 2) / √(V_n + σ²_ϵ)) -
                      normcdf(((n + 1) * a_grid[j] - n * a_grid[i] - μ_ns[i] - d_down / 2) / √(V_n + σ²_ϵ))
        end
    end
    # Renormalize the transition matrix
    if normalize == :yes
        for i = 1:N
            Π[i, :] = Π[i, :] ./ sum(Π[i, :])
        end
    end
    return Π
end

Transition_Probability(n::Integer, par::ParsGrowth; normalize::Symbol = :yes) = Transition_Probability(par::ParsGrowth, BPars(n, par); normalize)

function find_limit(par::ParsGrowth; tol = 1e-8)
    diff = 1e8
    n = 5
    while diff > tol
        Π_first = Transition_Probability(n, par)
        Π_next = Transition_Probability(n + 1, par)

        diff = sum((Π_first - Π_next) .^ 2)
        if diff > tol
            n += 1
        end
    end
    return n
end

"""
    Expectation of next exp(a/σ) having received n signals with mean a_bar.
        Seems correct.
"""
function EA(a_bar, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    exp(μ_n(a_bar, par, bpar) / σ + ϵ_bar / σ + 1 / (2 * σ^2) * (V_n + σ²_ϵ))
end

function EA_dep(μ_n, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    exp(μ_n / σ + ϵ_bar / σ + 1 / (2 * σ^2) * (V_n + σ²_ϵ))
end

"""
    Expected Profits given the optimal labor demand decision.
        Tested
"""
function Eπ(a_bar, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    ea = EA(a_bar, par, bpar)
    con = (σ - 1)^(σ - 1) / (σ^σ)
    u^(σ - 1) * ea^σ * con - w * f
end

Eπ(a_bar; par::ParsGrowth, bpar::BParsGrowth) = Eπ(a_bar, par, bpar)

function Eπ_dep(μ_n, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    ea  = EA_dep(μ_n, par, bpar)
    con = (σ - 1)^(σ - 1) / (σ^σ)
    u^(σ - 1) * ea^σ * con - w * f
end

Eπ_dep(μ_n; par::ParsGrowth, bpar::BParsGrowth) = Eπ_dep(μ_n, par, bpar)

"""
    Expected Revenue given the optimal labor demand decision.
        Tested
"""
function Er(a_bar, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    ea = EA(a_bar, par, bpar)
    con = ((σ - 1) / σ)^(σ - 1)
    u^(σ - 1) * ea^σ * con
end

Er(a_bar; par::ParsGrowth, bpar::BParsGrowth) = Er(a_bar, par, bpar)

function Er_dep(μ_n, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    ea = EA_dep(a_bar, par, bpar)
    con = ((σ - 1) / σ)^(σ - 1)
    u^(σ - 1) * ea^σ * con
end

Er_dep(μ_n; par::ParsGrowth, bpar::BParsGrowth) = Er_dep(μ_n, par, bpar)

"""
    Labor costs resulting from the optimal labor demand decision.
"""
function labor_costs(a_bar, par::ParsGrowth, bpar::BParsGrowth)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    ea = EA(a_bar, par, bpar)
    con = ((σ - 1) / σ)^σ
    u^(σ - 1) * ea^σ * con
end

labor_costs(a_bar; par::ParsGrowth, bpar::BParsGrowth) = labor_costs(a_bar, par, bpar)

"""
    For some n_max, compute the steady state value function.
"""
function find_V_n_max(n_max, par::ParsGrowth; bpar::BParsGrowth = BPars(n_max, par), tol = 1e-16)
    @unpack_ParsGrowth par
    Π = Transition_Probability(n_max, par)
    V_guess = zeros(N)
    V_new = zeros(N)
    diff = 1e8
    while diff > tol
        V_new = max.(Eπ.(a_grid; par, bpar) + β * (1 - δ) * Π * V_guess, 0)
        diff = sum((V_new - V_guess) .^ 2)
        diff > tol ? V_guess = copy(V_new) : nothing
    end
    return V_new
end

function find_V_n(n, V_n_plus_1, par::ParsGrowth; bpar::BParsGrowth = BPars(n, par))
    @unpack_ParsGrowth par
    Π = Transition_Probability(n, par)
    max.(Eπ.(a_grid, par = par, bpar = bpar) + β * (1 - δ) * Π * V_n_plus_1, 0)
end

"""
    Value function of an entrant with forced continuation. Necessary to find the cutoff for u.
"""
function find_V_0(V_1, par::ParsGrowth; bpar::BParsGrowth = BPars(0, par))
    @unpack_ParsGrowth par
    Π = Transition_Probability(0, par)
    Eπ.(a_grid, par = par, bpar = bpar) + β * (1 - δ) * Π * V_1
end

function find_all_n(n_max, V_n_max, par::ParsGrowth)
    V = zeros(par.N, n_max + 1)
    V[:, end] = V_n_max
    for n = n_max:-1:1
        V[:, n] = find_V_n(n - 1, V[:, n+1], par)
    end
    return V
end

function find_u_star(par::ParsGrowth)
    n_max = find_limit(par, tol = 1e-16)
    f(u) = begin
        @set! par.u = u
        V_n_max = find_V_n_max(n_max, par)
        V = find_all_n(n_max, V_n_max, par)
        return find_V_0(V[:, 2], par)[1] # This is the zeroth value function with forced continuation.
    end
    find_zero(f, 0.15)
end

function find_all_V_n(n_max, par::ParsGrowth)
    V_n_large = zeros(par.N, n_max + 1, par.N_z)
    for i = 1:par.N_z
        @set! par.u = par.u_grid[i]
        V_n_max = find_V_n_max(n_max, par)
        V = find_all_n(n_max, V_n_max, par)
        V_n_large[:, :, i] = V
    end
    return V_n_large
end

function u_to_z(par::ParsGrowth)
    @unpack_ParsGrowth par
    1 / (σ - 1) * log(u^(σ - 1) * w^(σ - 1) / (P^(σ - 1) * Y))
end

function z_to_u(z, par::ParsGrowth)
    @unpack_ParsGrowth par
    (exp(z)^(σ - 1) * P^(σ - 1) * Y / (w^(σ - 1)))^(1 / (σ - 1))
end

# function Firm_Distribution(par::ParsGrowth; normalize::Symbol = :yes)
#     @unpack_ParsGrowth par
#     Π    = zeros(N_z)
#     d    = z_grid[2] - z_grid[1] # equidistant grid
#     # Rows: Today's state i
#     # Columns: Tomorrow's state j
#     for i in 1:N_z
#         if i == 1
#             Π[i] = normcdf((z_grid[i] - z_bar + d/2)/√(σ²_z))
#         elseif i == N
#             Π[i] = 1 - normcdf((z_grid[i] - z_bar - d/2)/√(σ²_z))
#         else
#             Π[i] = normcdf((z_grid[i] - z_bar + d/2)/√(σ²_z)) - normcdf((z_grid[i] - z_bar - d/2)/√(σ²_z))
#         end
#     end
#     # Renormalize the transition matrix
#     if normalize == :yes
#             Π = Π./sum(Π)
#     end
#     return Π
# end

function Firm_Distribution(par::ParsGrowth; normalize::Symbol = :yes)
    @unpack_ParsGrowth par
    Π = zeros(N_z)
    d = (u_grid[2:end] - u_grid[1:end-1]) ./ 2
    # Rows: Today's state i
    # Columns: Tomorrow's state j
    for i = 1:N_z
        if i == 1
            Π[i] = cdf(Pareto(ξ, u_star), u_grid[i] + d[i])
        elseif i == N_z
            Π[i] = 1 - cdf(Pareto(ξ, u_star), u_grid[i] - d[i-1])
        else
            Π[i] = cdf(Pareto(ξ, u_star), u_grid[i] + d[i]) - cdf(Pareto(ξ, u_star), u_grid[i] - d[i-1])
        end
    end
    # Renormalize the transition matrix
    if normalize == :yes
        Π = Π ./ sum(Π)
    end
    return Π
end

"""
    How the population of firms evolves from one period to the next.
"""
function iterate_population(n_max, Π_set, V_set, par::ParsGrowth)
    # Initialize array for the distribution of firms
    pop = zeros(par.N, n_max)
    # Which firms survive / decide not to exit each period
    V_set[V_set.>0] .= 1
    # Set the first column of the firm distribution. All of these firms survive and see their first sales and learn from them.
    # A fraction δ dies after production in their first period.
    pop[:, 1] = (1 - par.δ) .* Π_set[1][1, :] # Distribution over different a_bar after the initial period

    # Then start from the second period (n=2) on and go to the last period.
    for n = 2:1:n_max
        # A fraction δ dies.
        # Then take the population from last period and take away all firms that decide to exit due to too negative demand realizations.
        # Mismatch in timing: V_set starts at period 0, pop starts in period 1. Therefore n-1 for pop and n for V_set.
        # Then study which firms go to another state.
        # Finally, the transition matrix is appropriately transposed.
        pop[:, n] = (1 - par.δ) .* Π_set[n]' * (pop[:, n-1] .* V_set[:, n])
    end
    return pop
end

function Equilibrium_Distribution(n_max, par::ParsGrowth)
    pop = zeros(par.N_z, par.N, n_max)
    Π_u = Firm_Distribution(par)
    Π_set = [Transition_Probability(par, BPars(n, par)) for n = 0:n_max] # find all Transition Matrices
    for i = 1:length(par.u_grid)
        @set! par.u = par.u_grid[i]
        V_set = find_all_n(n_max, find_V_n_max(n_max, par), par)# find all V for a given u
        pop[i, :, :] = Π_u[i] .* iterate_population(n_max, Π_set, V_set, par)
    end
    return pop
end

function Eπ_Distribution(n_max, par::ParsGrowth)
    Eπs = zeros(par.N_z, par.N, n_max)
    for i = 1:length(par.u_grid), j = 1:length(par.a_grid), n = 1:n_max
        @set! par.u = par.u_grid[i]
        bpar = BPars(n, par)
        Eπs[i, j, n] = Eπ(par.a_grid[j], par, bpar)
    end
    return Eπs
end

function Er_Distribution(n_max, par::ParsGrowth)
    Ers = zeros(par.N_z, par.N, n_max)
    for i = 1:length(par.u_grid), j = 1:length(par.a_grid), n = 1:n_max
        @set! par.u = par.u_grid[i]
        bpar = BPars(n, par)
        Ers[i, j, n] = Er(par.a_grid[j], par, bpar)
    end
    return Ers
end

"""
    Find the steady state and compute the amss of firms and the per-period entry.
"""
function mass_and_entry(par::ParsGrowth; n_max = 35)
    u_star = find_u_star(par)
    @show u_star

    @set! par.u = u_star
    @set! par.u_star = u_star

    u_max = quantile(Pareto(par.ξ, par.u_star), 0.9999)
    @set! par.u_grid = range(1.001 * par.u_star, u_max, length = par.N_z)

    #! Lets start with the firms on the first grid point of productivity
    pop = Equilibrium_Distribution(n_max, par)

    Eπs = Eπ_Distribution(n_max, par)
    Ers = Er_Distribution(n_max, par)

    total_pop = sum(pop) + 1 # the one is for the initial set of firms that did not observe any signals yet!

    init_π = Eπ(0, par, BPars(0, par))
    init_r = Er(0, par, BPars(0, par))

    avg_r = (sum(pop .* Ers) + init_r) / total_pop
    avg_π = (sum(pop .* Eπs) + init_π) / total_pop

    M = par.L / (avg_r - avg_π)

    z_under_bar = log((par.J * (sum(pop) + total_pop) / M)^(1 / par.ξ) * exp(par.z_min))

    entry = 1 - cdf(Pareto(par.ξ, par.z_min), z_under_bar)

    Y = par.L + M * avg_π
    P = par.u_star / (Y^(1 / (par.σ - 1)) * exp(z_under_bar))

    return StSt(M = M, E = entry, Y = Y, P = P)
end

function mean_preserving_noise(par, σ²_ϵ)
    @set! par.σ²_ϵ = σ²_ϵ
    @set! par.ϵ_bar = -par.σ²_ϵ / (par.σ * 2)
    return par
end


text:
gmsh_root = nothing

if haskey(ENV,"GMSHROOT") && !isempty(ENV["GMSHROOT"])
  gmsh_root = ENV["GMSHROOT"]
  @info """
  Using the gmsh installation found via the GMSHROOT environment variable.
  GMSHROOT=$gmsh_root"""
else
  gmsh = Sys.which("gmsh")
  if gmsh === nothing
    s = """
    The automatic discovery of a gmsh installation failed.
    Try one of the following:

    - Install gmsh and make sure that the gmsh executable
      can be found via the PATH environment variable.

    - Set an environemnt variable called GMSHROOT with the path
    to the location of a gmsh installation. Make sure that:
      GMSHROOT/bin/gmsh is the path of the gmsh binary and
      GMSHROOT/lib/gmsh.jl is the path of the gmsh Julia API.
    """
    @warn s
  else
    gmsh_root = dirname(gmsh)[1:end-4]
    @info """
    Using the gmsh installation automatically found in $gmsh_root"""
  end
end

GMSH_FOUND = false
if gmsh_root != nothing
  gmsh_bin = joinpath(gmsh_root,"bin","gmsh")
  if isfile(gmsh_bin)
    GMSH_FOUND = true
    @info "gmsh binary found in $gmsh_bin"
  else
    @warn """
    gmsh binary not found in $gmsh_bin
    Make sure that GMSHROOT points to a correct gmsh installation"""
  end
  gmsh_jl = joinpath(gmsh_root,"lib","gmsh.jl")
  if isfile(gmsh_jl)
    GMSH_FOUND = true
    @info "Using the gmsh Julia API found in $gmsh_jl"
  else
    @warn """
    gmsh Julia API not found in $gmsh_jl
    Make sure that you have installed the full gmsh SDK"""
  end
end

open("deps.jl","w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const GMSH_FOUND = $GMSH_FOUND))
  println(f, :(const gmsh_jl    = $gmsh_jl))
end


text:
@testset "divergence" begin

N = [4,5]
M = prod(N)
r = reshape(collect(1:2M), M, 2)'[:]  # index mapping from block matrix to narrowly banded matrix
Z = spzeros(M,M)

F = rand(ComplexF64, N..., 2)  # input array of vector field
g = zeros(ComplexF64, N...)  # output array of scalar
gvec = zeros(ComplexF64, M)  # column vector representation of g

@testset "create_divg and apply_divg!" begin
    for ci = CartesianIndices((false:true,false:true))
        # Construct Divg for a uniform grid and Bloch boundaries.
        isfwd = Vector{Bool}([ci.I...])
        Divg = create_divg(isfwd, N, order_cmpfirst=false)

        # Test the overall coefficients.
        @test size(Divg) == (M,2M)
        @test all(any(Divg.≠0, dims=1))  # no zero columns
        @test all(any(Divg.≠0, dims=2))  # no zero rows
        @test all(sum(Divg, dims=1) .== 0)  # all column sums are zero, because each input field to Divg is used twice, once multiplied with +1 and once with -1
        @test all(sum(Divg, dims=2) .== 0)  # all row sums are zero, because Divg * ones(M) = 0
        @test all(sum(abs.(Divg), dims=1) .== 2)  # each column of Divg has two nonzero entries, which are ±1's
        @test all(sum(abs.(Divg), dims=2) .== 4)  # each row of Divg has six nonzero entries, which are ±1's

        ∂x = (nw = 1; create_∂(nw, isfwd[nw], N))
        ∂y = (nw = 2; create_∂(nw, isfwd[nw], N))
        @test Divg == [∂x ∂y]

        # Construct Divg for a nonuniform grid and general boundaries.
        ∆l⁻¹ = rand.(tuple(N...))  # isfwd = true (false) uses ∆l⁻¹ at dual (primal) locations
        isbloch = [true, false]
        e⁻ⁱᵏᴸ = rand(ComplexF64, 2)

        Divg = create_divg(isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)

        # Test Divg.
        ∂x = (nw = 1; create_∂(nw, isfwd[nw], N, ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        ∂y = (nw = 2; create_∂(nw, isfwd[nw], N, ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        @test Divg == [∂x ∂y]

        # Test Cartesian-component-first ordering.
        Divg_cmpfirst = create_divg(isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=true)
        @test Divg_cmpfirst == Divg[:,r]

        # Test apply_divg!.
        f = F[:]
        mul!(gvec, Divg, f)
        apply_divg!(g, F, Val(:(=)), isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ)  # no keyword argument order_cmpfirst in apply_XXX!()
        @test g[:] ≈ gvec

        mul!(gvec, Divg_cmpfirst, f[r])
        apply_divg!(g, F, Val(:(=)), isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ)  # no keyword argument order_cmpfirst in apply_XXX!()
        @test g[:] ≈ gvec
    end
end  # @testset "create_divg and apply_divg!"

@testset "divergence of curl" begin
    # Construct Curl and Divg for a uniform grid and periodic boundaries.
    N = [3,4,5]
    M = prod(N)
    isfwd = [true, true, true]  # curl(U) and divg(V) are differentiated forward

    Curl = create_curl(isfwd, N, order_cmpfirst=false)
    Divg = create_divg(isfwd, N, order_cmpfirst=false)

    # Construct Divg * Curl.
    A = Divg * Curl

    # Test Divergence of curl.
    @test size(A) == (M,3M)
    @test all(A.==0)
end  # @testset "divergence of curl"

end  # @testset "divergence"


text:
module CrossMappings

using Statistics
using Distributions
using StatsBase

include("convergent_cross_mapping/validate_input.jl")
include("convergent_cross_mapping/crossmapping.jl")
include("convergent_cross_mapping/convergentcrossmapping.jl")

end # module


text:
#functions / variables for global units option. See set_units in ThreeBodyTB.jl

global global_energy_units="eV"
global global_length_units="Å"

const eV=13.6056980659
const Ang=0.529177210903

function convert_energy(energy)
    if global_energy_units=="eV"
        return energy * eV
    else
        return energy
    end
end

function convert_dos(dos)
    if global_energy_units=="eV"
        return dos / eV
    else
        return dos
    end
end

function convert_force(force)
    if global_energy_units=="eV"
        f1 = eV
    else
        f1 = 1.0
    end
    if global_length_units=="Å"
        f2 = 1/Ang
    else
        f2 = 1.0
    end

    return force * f1 * f2

end

function convert_stress(stress)
    if global_energy_units=="eV"
        f1 = eV
    else
        f1 = 1.0
    end
    if global_length_units=="Å"
        f2 = (1/Ang)^3
    else
        f2 = 1.0
    end

    return stress * f1 * f2

end

function convert_length(length)

    if global_length_units=="Å"
        f1 = Ang
    else
        f1 = 1.0
    end

    return length * f1

end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m4_opt_rho0_69(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m4_opt_rho0_69!(A_copy)
end

@inline function exp_ps_m4_opt_rho0_69!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 T2k7
    # Computing Bb2 = x*I+x*A
    coeff1=1.7372701751396924
    coeff2=-0.21601444678956824
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=1.7372701751396924
    coeff2=-0.21601444678956824
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=1.0508101696874834
    coeff2=1.040160435318157
    coeff3=2.474210559924759
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=2.92870970078987
    coeff2=1.395210200929235
    coeff3=0.9208910780611814
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.2866026571168137
    coeff2=-0.06808693499430804
    coeff3=-1.1770525298266241
    coeff4=2.745273367542917
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.6108235957535929
    coeff2=-0.014709720427431594
    coeff3=-0.020859019885585126
    coeff4=0.002429240292310237
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.3880998628583675
    coeff2=2.232643563510906
    coeff3=1.6627704342768925
    coeff4=1.104740854605295
    coeff5=0.5901519971324986
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-1.8244556096069469
    coeff2=1.0277285724068361
    coeff3=1.2223238419486635
    coeff4=0.7913945361070499
    coeff5=0.6037393214424426
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing T2k7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.37176853543253857
    coeff2=5.363140388026199
    coeff3=-1.0137337080611135
    coeff4=-1.6019799776198416
    coeff5=-2.338782896508547
    coeff6=1.6773088417874666
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    return memslots1 # Returning T2k7
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module Utilities

using ForwardDiff, CUDAnative, CuArrays, ..GPUUtils, JuAFEM, IterativeSolvers

export  AbstractPenalty,
        PowerPenalty,
        RationalPenalty,
        TopOptTrace,
        RaggedArray,
        @debug,
        compliance,
        meandiag,
        density,
        find_black_and_white,
        find_varind,
        YoungsModulus,
        PoissonRatio,
        getpenalty,
        getprevpenalty,
        setpenalty!,
        getsolver,
        @params

function getpenalty end
function getprevpenalty end
function setpenalty! end
getsolver(f) = f.solver

# Utilities
include("utils.jl")

# Trace definition
include("traces.jl")

# Penalty definitions
include("penalties.jl")

end


text:
using LinearAlgebra

abstract type AbstractInterface end
export AbstractInterface

abstract type AbstractBoundaryInterface <: AbstractInterface end
export AbstractBoundaryInterface

abstract type AbstractInternalInterface <: AbstractInterface end
export AbstractInternalInterface

struct EmptyInterface <: AbstractInterface end
export EmptyInterface

@with_kw struct IdealGasCatalystInterface{T<:AbstractPhase,N<:AbstractPhase,Q<:AbstractReaction} <: AbstractInternalInterface
    gas::T
    catalyst::N
    reactions::Array{Q,1}
end
export IdealGasCatalystInterface

struct Inlet{Q<:Real,V<:AbstractArray,U<:Real,X<:Real} <: AbstractBoundaryInterface
    F::V
    T::U
    P::X
    dUdt::Q
end

function Inlet(phase::V,conddict::Dict{String,X},F::B) where {V<:AbstractPhase,X<:Real,B<:Real}
    y = makespcsvector(phase,conddict)
    T = conddict["T"]
    P = conddict["P"]
    Fout = F.*y./sum(y)
    dUdt = dot(getEnthalpy.(getfield.(ph.species,:thermo),T),Fout)
    return Inlet(Fout,T,P,dUdt)
end

export Inlet

struct Outlet{V<:Real} <: AbstractBoundaryInterface
    F::V
end
export Outlet


text:
module FlexibleTypeMod

typealias FlexType UInt8

const INTEGER = FlexType(0)
const FLOAT = FlexType(1)
const STRING = FlexType(2)
const VECTOR = FlexType(3)
const LIST = FlexType(4)
const DICT = FlexType(5)
const DATETIME = FlexType(6)
const UNDEFINED = FlexType(7)
const IMAGE = FlexType(8)

using Cxx
import Cxx: CppValue
import Base: show, get, +, -, /, *, >, >=, ==, <, <=
import SFrames.Util: cstring, jlstring

immutable FlexibleType
    data::Int64
    buf::NTuple{4,UInt8}
    typ::FlexType
    buf2::NTuple{3,UInt8}
end

function FlexibleType(data, typ)
    FlexibleType(data, (0x00, 0x00, 0x00, 0x00), typ, (0x00, 0x00, 0x00))
end

function FlexibleType(x::Int)
    FlexibleType(x, INTEGER)
end

function FlexibleType(x::AbstractString)
    icxx"return flexible_type($(cstring(x)));" |> FlexibleType
end

function FlexibleType(c_flex_type::Cxx.CppValue)
    bytes = c_flex_type.data
    # todo make a generated function for this
    data = reinterpret(Int, UInt64(bytes[1]) | UInt64(bytes[2])<<8 | UInt64(bytes[3])<<16 | UInt64(bytes[4])<<24) | UInt64(bytes[5])<<32 | UInt64(bytes[6])<<40 | UInt64(bytes[7])<<48 | UInt64(bytes[8])<<56
    typ = bytes[13]
    FlexibleType(data, typ)
end

function FlexibleType(x::Float64)
    FlexibleType(reinterpret(Int,x), FLOAT)
end

const FLEX_NULL = FlexibleType(0, UNDEFINED)

function Base.convert(::Type{Int}, f::FlexibleType)
    f.data
end

function Base.convert(::Type{Float64}, f::FlexibleType)
    reinterpret(Float64, f.data)
end

function Base.bytestring(f::FlexibleType)
    ref = icxx"((flexible_type*)$(pointer_from_objref(f)))->get<flex_string>();"
    icxx"$ref.size();" |> Int
end

r=FlexibleType("jon")

ptr=icxx"(flexible_type*)($(pointer_from_objref(r)));"
t=icxx"*$ptr;"
icxx"((flex_string)$t).c_str();"

function Base.get(f::FlexibleType)
    if f.typ == INTEGER
        Int(f)
    elseif f.typ == FLOAT
        Float64(f)
    elseif f.typ == UNDEFINED
        nothing
    elseif f.typ == STRING
        "hi"
    end
end

function FlexibleType(x::FlexibleType)
    FlexibleType(x.data, x.typ)
end

function show(io::IO, t::FlexibleType)
    print(io, "FlexibleType(")
    show(io, get(t))
    print(io, ")")
end

#  todo this is just a rough sketch of defining operators on flex types
for op in [:+, :-, :/, :*]
    @eval function $op(f1::FlexibleType, f2::FlexibleType)
        if f1.typ == INTEGER && f2.typ == INTEGER
            FlexibleType($op(Int(f1), Int(f2)))
        elseif f1.typ == UNDEFINED || f2.typ == UNDEFINED
            FLEX_NULL
        else
            FlexibleType($op(Float64(f1), Float64(f2)))
        end
    end
end

for op in [:(==), :>, :>=, :<, :<=]
    @eval function $op(f1::FlexibleType, f2::FlexibleType)
        f1.typ == f2.typ || error("Invalid comparsion")
        $op(f1.data, f2.data)
    end
end

macro dispatch_numeric(funcs...)
    block = Expr(:block)
    for func in funcs
        e = quote
            function $(esc(func))(f::FlexibleType)
                if f.typ == INTEGER
                    $(esc(func))(Int(f)) |> FlexibleType
                elseif f.typ == FLOAT
                    $(esc(func))(Float64(f)) |> FlexibleType
                else
                    error("Invalid type for operation")
                end
            end
        end
        push!(block.args, e)
    end
    block
end
@dispatch_numeric Base.sin Base.cos Base.tan Base.log Base.exp Base.sqrt Base.mod Base.mod1 Base.expm1 (-)

# for (julia_op, c_op) in [
#     (:+, "+"),
#     (:-, "-"),
#     (:/, "/"),
#     (:*, "*"),
#     (:>, ">"),
#     (:>=, ">="),
#     (:(==), "=="),
#     (:<, "<"),
#     (:<=, "<=")
#     ]
#     cstr = "\$(s1.val) $c_op \$(s2.val);"
#     cstr_scalar = "\$(s1.val) $c_op \$s2;"
#     @eval  function $julia_op(s1::FlexibleType, s2::FlexibleType)
#         FlexibleType($(Expr(:macrocall, Symbol("@icxx_str"), cstr)))
#     end
#
#     @eval function $julia_op(s1::FlexibleType, s2::Int)
#         FlexibleType($(Expr(:macrocall, Symbol("@icxx_str"), cstr_scalar)))
#     end
#
#     @eval $julia_op(s1::Int, s2::FlexibleType) = $julia_op(s2, s1)
# end

end


text:
using Documenter, WAV

makedocs(sitename="WAV.jl")


text:
"""
    x_ray_emission(n_cm3::Real, T_K::Real;
                   Emin::Real=5.e4, Emax::Real=1.e10,
                   xH::Real=0.76)

X-Ray emission of a particle with number density `n_cm3` in ``1/cm^3`` and temperature `T` in ``K``.
`Emin` and `Emax` give the minimum and maximum energy of the oberservation.
`xH` gives the hydrogen fraction used in the simulation.

# Returns
X-Ray surface brightness contribution in units of [erg/cm^2/s/Hz].
"""
function x_ray_emission(n_cm3::Real, T_K::Real;
                        Emin::Real=5.e4, Emax::Real=1.e10,
                        xH::Real=0.76)

    prefac = 4.0 * 2.42e-24 / (1 + xH)
    T_eV   = T_K * cgs2eV

    return prefac * sqrt( k_B * T_eV * 1.e-3 ) * n_cm3^2 *
            ( exp( -Emin / (k_B * T_eV) ) - exp( -Emax / (k_B * T_eV) ) )
end


text:
include("../src/HartreeFock.jl")

using .HartreeFock
using BenchmarkTools

# define some unit vector
const n = (
    vec3(  1.,    0., 0.),
    vec3( 1/2,  √3/2, 0.),
    vec3(-1/2,  √3/2, 0.),
    vec3( -1.,    0., 0.),
    vec3(-1/2, -√3/2, 0.),
    vec3( 1/2, -√3/2, 0.)
)

# dCC: bond length between adjacent `C` Atom
# dCH: bond length between `C` and `H`
function Benzene(dCC::Float64, dCH::Float64)::Molecule
    Molecule("C6H6", vcat(
        collect(ntuple(i->Atom("C", dCC * n[i]), 6)),
        collect(ntuple(i->Atom("H", (dCC+dCH) * n[i]), 6))
    ))
end

function Benzene_pes(dCC_range, dCH_range)
    NCC = length(dCC_range)
    NCH = length(dCH_range)
    Emat = zeros(Float64, NCC, NCH)
    Threads.@threads for i in 1:NCC
    # for i in 1:NCC
        for j in 1:NCH
            benzene = Benzene(dCC_range[i] / a0, dCH_range[j] / a0)
            total_energy = hartree_fock(benzene)
            Emat[i, j] = total_energy
        end
    end
    i, j = argmin(Emat).I
    println("minimal energy is $(Emat[i,j]), with dCC = $(dCC_range[i])Å, dCH = $(dCH_range[j])Å")
end

#display(@benchmark Benzene_pes(1.3:0.02:1.5, 1:0.02:1.2))
Benzene_pes(1.3:0.02:1.5, 1:0.02:1.2)


text:
import gmsh

model = gmsh.model
factory = model.occ

gmsh.initialize(ARGS)

# creates the FLTK user interface; this could also be called after the geometry
# is created (or not at all - fltk.run() will do it automatically)
gmsh.fltk.initialize()

# Copied from boolean.py...
model.add("boolean")
gmsh.option.setNumber("Mesh.Algorithm", 6);
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4);
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4);
R = 1.4; Rs = R*.7; Rt = R*1.25
factory.addBox(-R,-R,-R, 2*R,2*R,2*R, 1)
factory.addSphere(0,0,0,Rt, 2)
factory.intersect([(3, 1)], [(3, 2)], 3)
factory.addCylinder(-2*R,0,0, 4*R,0,0, Rs, 4)
factory.addCylinder(0,-2*R,0, 0,4*R,0, Rs, 5)
factory.addCylinder(0,0,-2*R, 0,0,4*R, Rs, 6)
factory.fuse([(3, 4), (3, 5)], [(3, 6)], 7)
factory.cut([(3, 3)], [(3, 7)], 8)
factory.synchronize();
# ...end of copy

# this would be equivalent to gmsh.fltk.run():
#
# gmsh.graphics.draw()
# while true
#     gmsh.fltk.wait()
#     println("just treated an event in the interface")
# end

gmsh.fltk.run()

gmsh.finalize()

max_stars_repo_name: Poofee/fastFEM

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
"""
Implementation of the Kronecker PCA method for covariance estimation

Paper: Tsiligkaridis T, Hero AO.
    Covariance estimation in high dimensions via kronecker product expansions.
     IEEE Transactions on Signal Processing. 2013 Aug 22;61(21):5347-60.

Author: Wayne Wang
Last modified: Apr 1, 2021
"""

function rearrange(A::AbstractMatrix{<:Real},
    m1::Int, n1::Int, m2::Int, n2::Int)
    """
    Rearrange operator of (Van Loan, 1993)

    Inputs:
        A: input matrix comprising m1 × n1 blocks of m2 × n2 matrices
    Output:
        R: rearranged matrix comprising n1 × n2 blocks of m1 × m2 matrices
    """
    R = zeros((m1*n1,m2*n2))
    Aij = zeros((m2,n2))

    nrow = 1
    for j=1:n1
        for i=1:m1
            Aij .= view(A, ((i-1)*m2+1):i*m2, ((j-1)*n2+1):j*n2)
            R[nrow,:] .= Aij[:]
            nrow += 1
        end
    end

    return R
end

function rearrange_adj(R::AbstractMatrix{<:Real},
    m1::Int, n1::Int, m2::Int, n2::Int)
    """
    Adjoint operator of the rearrange operator (Van Loan, 1993)

    Inputs:
        R: rearranged matrix comprising n1 × n2 blocks of m1 × m2 matrices
    Output:
        A: input matrix comprising m1 × n1 blocks of m2 × n2 matrices
    """
    A = zeros((m1*m2,n1*n2))
    Aij = zeros(m2*n2)

    nrow = 1
    for j=1:n1
        for i=1:m1
            Aij .= view(R, nrow, :)
            A[((i-1)*m2+1):i*m2, ((j-1)*n2+1):j*n2] .= reshape(Aij,(m2,n2))
            nrow += 1
        end
    end

    return A
end

function kron_pca(S::AbstractMatrix{<:Real}, k::Int,
    px::Int, py::Int; λ::Real=0)
    """
    PCA on the rearranged sample covariance of data matrix X

    Inputs:
        S: sample covariance matrix d × d, where d = pxpy
        k: top-k PCA/SVD
    Optional:
        λ: threshold value for singular values
    Output:
        A, B: top-k left (A of size px²×k) and right (B of size py²×k)
         singular vectors
    """
    # rearrange the sample covariance
    R = rearrange(S, px, px, py, py)

    # topc-k SVD for rearranged sample covariance
    F = svds(R, nsv=k)[1]
    U = F.U
    Vt = F.Vt
    svdvals = F.S

    # form covariance estiamte
    SigmaH = zeros(size(S))
    for j=1:k
        SigmaH .+= max(svdvals[j]-λ/2,0)*kron(reshape(view(U,:,j),(px,px)),
         reshape(view(Vt,j,:),(py,py)))
    end

    return SigmaH
end

function robust_pca(X::AbstractMatrix{<:Real}, lambdaL::Real, lambdaS::Real,
     method::String; tau::Real=0.5, iter::Int=5, r::Int=5)
    """
    Perform robust PCA as in Nadakuditi et al 2014 (SSP). Proximal
    gradient method.
    """
    # dims
    c = size(X,2)/size(X,1)
    q = minimum(size(X))

    # initialize
    L = zeros(size(X))
    S = spzeros(size(X,1),size(X,2))
    M = deepcopy(X)
    U = zeros((size(X,1),q))
    Vt = zeros((q,size(X,2)))
    s = zeros(q)

    for i=1:iter
        # SVT
        if method == "SVT"
            F = svd(M-S)
            U .= F.U
            Vt .= F.Vt
            s .= F.S
            r = count(i -> i>lambdaL*tau, s)

            if r == 0
                r = 1
            end

            if s[1] == 0
                s[1] = 1e-10
            end
        elseif method == "OptShrink"
            F = svd(M-S)
            U .= F.U
            Vt .= F.Vt
            s .= F.S

            for i = 1:r
                sumV = 1/(q-r)*sum(s[i]./(s[i]^2 .- s[r+1:q].^2))
                D = sumV*(c*sumV+(1-c)/s[i])
                sumVprime = 1/(q-r)*sum((-s[i]^2 .- s[r+1:q].^2)./(s[i]^2 .- s[r+1:q].^2).^2)
                Dprime = (c*sumV+(1-c)/s[i])*sumVprime + (c*sumVprime-(1-c)/s[i]^2)*sumV
                s[i] = -2*D/Dprime
            end
        end

        L = U[:,1:r]*Diagonal(s[1:r])*Vt[1:r,:]

        # Soft
        S .= M .- L
        S .= sign.(S).*max.(abs.(S).-lambdaS*tau,0)

        # Update M
        M .= L .+ S .- tau*(L .+ S .- X)
    end

    return L, S, U[:,1:r], s[1:r], Vt[1:r,:]
end

function robust_kron_pca(S::AbstractMatrix{<:Real}, px::Int, py::Int,
    lambdaL::Real, lambdaS::Real, method::String;
    tau::Real=0.5, iter::Int=5, r::Int=5)
    """
    Robust Kronecker PCA method for covariance estimation
    """
    # rearrange the sample covariance
    R = rearrange(S, px, px, py, py)

    # robust pca for rearranged sample covariance
    L, S, _, _, _ = robust_pca(R, lambdaL, lambdaS, method; tau, iter, r)

    # form covariance estiamte
    SigmaH = rearrange_adj(L+S, px, px, py, py)

    return SigmaH
end

# ## Testing
# using Distributions
# using PDMats
# using BenchmarkTools

# px = py = 25
# N = 100
# X = zeros((px*py,N))

# for i=1:N
#     X[:,i] .= vec(rand(MatrixNormal(zeros((px,py)),ScalMat(px,2.0),ScalMat(py,4.0))))
# end

# S = cov(copy(X'))
# lambdaL = 20*(px^2+py^2+log(max(px,py,N)))/N
# lambdaS = 20*sqrt(log(px*py)/N)

# begin
#     @btime SigmaH = robust_kron_pca(S,px,py,lambdaL,lambdaS,"OptShrink";tau=0.5,r=5)
#     @btime SigmaH_kron_pca = kron_pca(S,5,px,py)
# end

# norm(8.0*I(px*py) - SigmaH)
# norm(8.0*I(px*py) - SigmaH_kron_pca)

# isposdef((SigmaH+SigmaH')/2.0)


text:
using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq

using Symbolics: scalarize

function cart_pole(u, p, t)
    du = similar(u)
    F = -0.2 + 0.5*sin(6*t) # the input
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -(19.62*sin(u[1])+sin(u[1])*cos(u[1])*u[3]^2+F*cos(u[1]))/(2-cos(u[1])^2)
    du[4] = -(sin(u[1])*u[3]^2 + 9.81*sin(u[1])*cos(u[1])+F)/(2-cos(u[1])^2)
    return du
end

u0 = [0.3; 0; 1.0; 0]
tspan = (0.0, 5.0)
dt = 0.1
cart_pole_prob = ODEProblem(cart_pole, u0, tspan)
solution = solve(cart_pole_prob, Tsit5(), saveat = dt)

# Create the differential data
X = solution[:,:]
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = cart_pole(xi, [], solution.t[i])
end
t = solution.t

ddprob = ContinuousDataDrivenProblem(
    X , t, DX = DX[3:4, :], U = (u,p,t) -> [-0.2 + 0.5*sin(6*t)]
)

@variables u[1:4] du[1:2] x[1:1] t
u = scalarize(u)
du = scalarize(du)
x = scalarize(x)

polys = polynomial_basis(u, 2)
push!(polys, sin.(u[1]))
push!(polys, cos.(u[1]))
push!(polys, sin.(u[1])^2)
push!(polys, cos.(u[1])^2)
push!(polys, sin.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*u[3:4].^2...)
push!(polys, sin.(u[1]).*cos.(u[1])...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4].^2...)
implicits = [du;  du[1] .* u; du[2] .* u; du .* cos(u[1]);   du .* cos(u[1])^2; polys]
push!(implicits, x...)
push!(implicits, x[1]*cos(u[1]))
push!(implicits, x[1]*sin(u[1]))

basis= Basis(implicits, [u; du], controls = x,  iv = t)

# Simply use any optimizer you would use for sindy
λ = [1e-4;5e-4;1e-3;2e-3;3e-3;4e-3;5e-3;6e-3;7e-3;8e-3;9e-3;1e-2;2e-2;3e-2;4e-2;5e-2;
6e-2;7e-2;8e-2;9e-2;1e-1;2e-1;3e-1;4e-1;5e-1;6e-1;7e-1;8e-1;9e-1;1;1.5;2;2.5;3;3.5;4;4.5;5;
6;7;8;9;10;20;30;40;50;100;200];

opt = ImplicitOptimizer(λ)
# AICC
ĝ(x) = x[1] <= 1 ? Inf : 2*x[1]-2*log(x[2])
res = solve(ddprob, basis, opt, du, maxiter = 1000, g = ĝ, scale_coefficients = true)
m = metrics(res)

@test m.Sparsity == 10
@test m.AICC < 180.0
@test m.Error < 1.0


text:
abstract type LinearSolver{T} end

factorize!(s::LinearSolver) = error("factorize! not implemented for $(typeof(s))")
LinearAlgebra.ldiv!(s::LinearSolver) = error("ldiv! not implemented for $(typeof(s))")

function LinearSolver(x::AbstractVector{T}; linear_solver = :julia) where {T}
    n = length(x)

    if linear_solver === nothing || linear_solver == :julia
        linear_solver = LUSolver{T}(n)
    elseif linear_solver == :lapack
        linear_solver = LUSolverLAPACK{T}(BlasInt(n))
    else
        @assert typeof(linear_solver) <: LinearSolver{T}
        @assert n == linear_solver.n
    end
    return linear_solver
end


text:
# Autogenerated wrapper script for Ogg_jll for x86_64-apple-darwin
export libogg

JLLWrappers.@generate_wrapper_header("Ogg")
JLLWrappers.@declare_library_product(libogg, "@rpath/libogg.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libogg,
        "lib/libogg.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Tokenize

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

include("token.jl")
include("lexer.jl")

import .Lexers: tokenize
import .Tokens: untokenize

export tokenize, untokenize, Tokens

# disable precompilation when profiling runtime performance, as
# it can lead to wrong traces
include("_precompile.jl")
_precompile_()

end # module