text:
#TODO: use proper data type(s) here
#       to receive message from ROS (and actually from Gazebo)
function callback(msg::Float64MultiArray)
   println(msg)
end
# https://github.com/jdlangs/RobotOS.jl

# TODO: here is an example of pushing data in fixed rate
#       (just like there is a timer)
function loop(pub_obj)
   loop_rate = Rate(1)
   npt = Float64MultiArray()
   npt.data = Array{Float64}(6)
   npt.data[6] = 0.0
   while ! is_shutdown()
     npt.data[1] += 0.2
     npt.data[2] = 0.0
     npt.data[3] = 0.0
     npt.data[4] = 0.0
     npt.data[5] = 0.0
     npt.data[6] += 0.2  # NOTE this might need to be set to 0

     # send information to ROS (then Gazebo)
     publish(pub_obj, npt)
     rossleep(loop_rate)
   end
end

function main()
   #init_node("rosjl_example") # what is this?

   pub = Publisher{Float64MultiArray}("gazebo_client/test", queue_size=10)
   sub = Subscriber{Float64MultiArray}("gazebo_client/sendscan", callback, (), queue_size=10)

   loop(pub)
end

if ! isinteractive()
   main()
end


text:
import .Interpolations: AbstractInterpolation, AbstractInterpolationWrapper, BSplineInterpolation,
                        LanczosInterpolation, ScaledInterpolation, GriddedInterpolation, MonotonicInterpolation,
                        Extrapolation, itptype, coefficients,  itpinfo, value_weights, WeightedAdjIndex,
                        allbetween, interp_getindex, WeightedIndex, weightedindex_parts, maybe_weightedindex,
                        BoundsCheckStyle, _checkbounds
import Base: @propagate_inbounds
export UnSafeInterp, preweight

# a wrapper to turn off the boundary check when we use broadcast
struct UnSafeInterp{T,N,ITPT,IT} <: AbstractInterpolationWrapper{T,N,ITPT,IT}
    itp::ITPT
    UnSafeInterp(itp::AbstractInterpolation{T,N,IT}) where {T,N,IT} = new{T,N,typeof(itp),IT}(itp)
end
Base.parent(usitp::UnSafeInterp) = usitp.itp
@inline (usitp::UnSafeInterp{T,N})(args::Vararg{Any,N}) where {T,N} = @inbounds usitp.itp(args...)

struct WeightedInterp{A}
    coefs::A
end
Base.parent(witp::WeightedInterp) = witp.coefs
@inline (witp::WeightedInterp{<:AbstractArray{T,N}})(args::Vararg{Any,N}) where {T,N} = @inbounds witp.coefs[args...]

# pre_weight
unsaled(r::AbstractUnitRange, x) = @lzb x .- first(r) .+ oneunit(eltype(r))
unsaled(r::AbstractRange, x) = @lzb (x .- first(r)) .* inv(step(r)) .+ oneunit(eltype(r))

@inline getcoefs(x::AbstractInterpolationWrapper) = parent(x) |> getcoefs
@inline getcoefs(x::AbstractInterpolation) = coefficients(x)
@inline getcoefs(x::WeightedInterp) = parent(x)

# preweight
@inline preweight(itp::Extrapolation{T,N}, args::Vararg{Any,N}) where {T,N} =
    throw(ArgumentError("Extrapolation is not supported"))
@inline preweight(sitp::ScaledInterpolation{T,N}, args::Vararg{Any,N}) where {T,N} = begin
    @boundscheck _checkbounds(BoundsCheckStyle(sitp), sitp, args) ||
        Base.throw_boundserror(sitp, args)
    @inbounds preweight(sitp.itp, unsaled.(sitp.ranges, args)...)
end
@inline preweight(itp::BSplineInterpolation{T,N}, args::Vararg{Any,N}) where {T,N} = begin
    @boundscheck _checkbounds(BoundsCheckStyle(itp), itp, args) ||
        Base.throw_boundserror(itp, args)
    function weight_ind(itpflag, knotvec, x)
        makewi(y, ::Any) = begin
            pos, (coefs,) = weightedindex_parts((value_weights,), itpflag, knotvec, y)
            maybe_weightedindex(pos, coefs)
        end
        makewi.(x, getcoefs(itp) |> Ref)
    end
    (WeightedInterp(itp.coefs), weight_ind.(itpinfo(itp)..., args)...)
end

@require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin

    import Adapt: adapt_structure, adapt
    adapt_structure(to, itp::BSplineInterpolation{T,N}) where {T,N} = begin
        coefs, parentaxes, it = itp.coefs, itp.parentaxes, itp.it
        coefs′ = adapt(to, coefs)
        Para = map(typeof, (coefs′,it,parentaxes))
        BSplineInterpolation{T,N,Para...}(coefs′, parentaxes, it)
    end

    adapt_structure(to, itp::LanczosInterpolation{T,N}) where {T,N} = begin
        coefs, parentaxes, it = itp.coefs, itp.parentaxes, itp.it
        coefs′ = adapt(to, coefs)
        parentaxes′ = adapts(to, parentaxes...)
        Para = map(typeof, (it,coefs′,parentaxes′))
        LanczosInterpolation{T,N,Para...}(coefs′, parentaxes′, it)
    end

    adapt_structure(to, itp::ScaledInterpolation{T,N}) where {T,N} = begin
        ranges = itp.ranges
        itp′ = adapt(to, itp.itp)
        IT = itptype(itp)
        ScaledInterpolation{T,N,typeof(itp′),IT,typeof(ranges)}(itp′, ranges)
    end

    adapt_structure(to, itp::Extrapolation{T,N}) where {T,N} = begin
        et = itp.et
        itp′ = adapt(to,  itp.itp)
        IT = itptype(itp)
        Extrapolation{T,N,typeof(itp′),IT,typeof(et)}(itp′, et)
    end

    adapt_structure(to, itp::UnSafeInterp{T,N}) where {T,N} =
        adapt(to, itp.itp) |> UnSafeInterp

    adapt_structure(to, itp::WeightedInterp) =
        adapt(to, itp.coefs) |> WeightedInterp

    # With Adapt v"3.3.1", there's no need to use Ref to force adapt
    # but I think the style hack is still needed.
    broadcasted(itp::Union{WeightedInterp, AbstractInterpolation}, args...) = begin
        args′ = broadcastable.(args)
        style = combine_styles(Ref(getcoefs(itp)), args′...)
        broadcasted(style, itp, args′...)
    end

    #make checkbounds work on gpu
    allbetween(l::Real, xs::AbstractArray, u::Real) = begin
        device(xs) isa GPU && return all(x -> l <= x <= u, xs)
        invoke(allbetween, Tuple{typeof(l),Any,typeof(u)}, l, xs, u)
    end

    # Type interfer and display fix
    import GPUArrays
    import GPUArrays: AbstractGPUArray
    GPUArrays._getindex(A::AbstractGPUArray{T,N}, I::Vararg{Union{Int,WeightedIndex},N}) where {T,N} =
        interp_getindex(A, I, ntuple(Returns(0), Val(N))...)

end

ext: jl


text:
using CovarianceMatrices
using Random
using GLM
using DataFrames
using BenchmarkTools

Random.seed!(1)

df = DataFrame(y = randn(1000))
for j in Symbol.("x".*string.(collect(1:5)))
    df[j] = randn(1000)
end

frm = @formula(y ~ x1 + x2 + x3 + x4 + x5)
lm1 = glm(frm, df, Normal(), IdentityLink())
k_fix = TruncatedKernel(2)
k_opt = TruncatedKernel()

suite = BenchmarkGroup()
suite["HAC"] = BenchmarkGroup(["Optimal Uncached", "Fixed Uncached"])
suite["HAC"]["Optimal Uncached"] = @benchmarkable vcov(lm1, $k_opt)
suite["HAC"]["Fixed Uncached"] = @benchmarkable vcov(lm1, $k_fix)


text:
import MsgPack
import UUIDs: UUID
import HTTP
import Sockets

"Will hold all 'response handlers': functions that respond to a WebSocket request from the client. These are defined in `src/webserver/Dynamic.jl`."
const funkies = Dict{Symbol,Function}()

macro expose(funcdef::Expr)
    quote
        funky = $(esc(funcdef))
        funkies[nameof(funky)] = funky
    end
end

"Attempts to find the MIME pair corresponding to the extension of a filename. Defaults to `text/plain`."
function mime_fromfilename(filename)
    # This bad boy is from: https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
    mimepairs = Dict(".aac" => "audio/aac", ".bin" => "application/octet-stream", ".bmp" => "image/bmp", ".css" => "text/css", ".csv" => "text/csv", ".eot" => "application/vnd.ms-fontobject", ".gz" => "application/gzip", ".gif" => "image/gif", ".htm" => "text/html", ".html" => "text/html", ".ico" => "image/vnd.microsoft.icon", ".jpeg" => "image/jpeg", ".jpg" => "image/jpeg", ".js" => "text/javascript", ".json" => "application/json", ".jsonld" => "application/ld+json", ".mjs" => "text/javascript", ".mp3" => "audio/mpeg", ".mpeg" => "video/mpeg", ".oga" => "audio/ogg", ".ogv" => "video/ogg", ".ogx" => "application/ogg", ".opus" => "audio/opus", ".otf" => "font/otf", ".png" => "image/png", ".pdf" => "application/pdf", ".rtf" => "application/rtf", ".sh" => "application/x-sh", ".svg" => "image/svg+xml", ".tar" => "application/x-tar", ".tif" => "image/tiff", ".tiff" => "image/tiff", ".ttf" => "font/ttf", ".txt" => "text/plain", ".wav" => "audio/wav", ".weba" => "audio/webm", ".webm" => "video/webm", ".webp" => "image/webp", ".woff" => "font/woff", ".woff2" => "font/woff2", ".xhtml" => "application/xhtml+xml", ".xml" => "application/xml", ".xul" => "application/vnd.mozilla.xul+xml", ".zip" => "application/zip")
    file_extension = getkey(mimepairs, '.' * split(filename, '.')[end], ".txt")
    MIME(mimepairs[file_extension])
end

function assetresponse(path)
    try
        @assert isfile(path)
        response = HTTP.Response(200, read(path, String))
        push!(response.headers, "Content-Type" => string(mime_fromfilename(path)))
        push!(response.headers, "Access-Control-Allow-Origin" => "*")

        response
    catch e
        HTTP.Response(404, "Not found!: $(e)")
    end
end

function serve_onefile(path)
    return request::HTTP.Request -> assetresponse(normpath(path))
end

function serve_asset(req::HTTP.Request)
    reqURI = req.target |> HTTP.URIs.unescapeuri |> HTTP.URI

    filepath = joinpath(pwd(), relpath(reqURI.path, "/"))
    assetresponse(filepath)
end

# https://github.com/JuliaWeb/HTTP.jl/issues/382
const streamtokens = WeakKeyDict{IO,Channel}()

function create_flushtoken()
    flushtoken = Channel{Nothing}(1)
    put!(flushtoken, nothing)
    return flushtoken
end

const MSG_DELIM = "IUUQ.km jt ejggjdvmu vhi" |> codeunits |> collect # riddle me this, Julius

function write_serialized(io::IO, x::Any)
    write(io, MsgPack.pack(x), MSG_DELIM)

end

function Base.endswith(vec::Vector{T}, suffix::Vector{T}) where T
    local liv = lastindex(vec)
    local lis = lastindex(suffix)
    liv >= lis && (view(vec, (liv - lis + 1):liv) == suffix)
end

function Base.readuntil(stream::HTTP.WebSockets.WebSocket, delim::Vector{UInt8})
    data = UInt8[]
    while !endswith(data, MSG_DELIM)
        if (eof(stream))
            if isempty(data)
                @warn "What is this"
                return data
            end
            @warn "Unexpected eof after" data
            return data
        end
        push!(data, readavailable(stream)...)
    end
    return data[1:end - length(delim)]
end

"""Start a MARGO server _synchronously_ (i.e. blocking call) on `http://localhost:[port]/`.

This will start a WebSocket server. MARGO Notebooks will be started dynamically (by user input)."""
function run(host, port::Integer)
    hostIP = parse(Sockets.IPAddr, host)
    serversocket = Sockets.listen(hostIP, UInt16(port))
    servertask = @async HTTP.serve(hostIP, UInt16(port), stream=true, server=serversocket) do http::HTTP.Stream
        # messy messy code so that we can use the websocket on the same port as the HTTP server

        if HTTP.WebSockets.is_upgrade(http.message)
            try
                HTTP.WebSockets.upgrade(http) do clientstream::HTTP.WebSockets.WebSocket
                    if !isopen(clientstream)
                        return
                    end
                    while !eof(clientstream)
                        # This stream contains data received over the WebSocket.
                        # It is formatted and encoded by connection.mjs
                        try
                            parentbody = let
                                data = readuntil(clientstream, MSG_DELIM)
                                MsgPack.unpack(data)
                            end
                            # Pass the message to one of the user-defined response functions
                            process_ws_message(parentbody, clientstream)
                        catch ex
                            if ex isa InterruptException
                                rethrow(ex)
                            elseif ex isa HTTP.WebSockets.WebSocketError
                                # that's fine!
                            else
                                bt = stacktrace(catch_backtrace())
                                @warn "Reading WebSocket client stream failed for unknown reason:" exception = (ex, bt)
                            end
                        end
                    end
                end
            catch ex
                if ex isa InterruptException
                    rethrow(ex)
                elseif ex isa Base.IOError
                    # that's fine!
                elseif ex isa ArgumentError && occursin("stream is closed", ex.msg)
                    # that's fine!
                else
                    bt = stacktrace(catch_backtrace())
                    @warn "HTTP upgrade failed for unknown reason" exception = (ex, bt)
                end
            end
        else
            # not a WS connection, so we serve a local file:
            request::HTTP.Request = http.message
            HTTP.closeread(http)
            response_body = serve_asset(http.message)

            request.response::HTTP.Response = response_body
            request.response.request = request

            HTTP.startwrite(http)
            write(http, request.response.body)
            HTTP.closewrite(http)
        end
    end

    # print address to console:
    root_url = get(ENV, "PLUTO_ROOT_URL", "/")
    address = if root_url == "/"
        hostPretty = (hostStr = string(hostIP)) == "127.0.0.1" ? "localhost" : hostStr
        portPretty = Int(port)
        "http://$(hostPretty):$(portPretty)/"
    else
        root_url
    end
    println("Go to $address to start writing ~ have fun!")
    println()
    println("Press Ctrl+C in this terminal to stop MARGO")
    println()

    # create blocking call:
    try
        wait(servertask)
    catch e
        if isa(e, InterruptException)
            println("\n\nClosing MARGO... Restart Julia for a fresh session. \n\nHave a nice day! 🎈")
            close(serversocket)
        else
            rethrow(e)
        end
    end
end

run(port::Integer=1234; kwargs...) = run("127.0.0.1", port; kwargs...)

# MsgPack returns numbers in the most efficient type (Int8, ..., UInt64),
# we always want Int64 to keep things simple
function withmorebits(d::Dict)
    Dict((p.first => withmorebits(p.second)
        for p in d))
end
function withmorebits(x::T) where T <: Union{Signed,Unsigned}
    Int64(x)
end
function withmorebits(x::T) where T <: AbstractFloat
    Float64(x)
end
function withmorebits(x::Vector)
    withmorebits.(x)
end
function withmorebits(x::Any)
    x
end

"Respond to a message from JS."
function process_ws_message(parentbody::Dict{Any,Any}, clientstream::HTTP.WebSockets.WebSocket)
    client_id = Symbol(parentbody["client_id"])

    # every stream has a `token` which only one async task can have at one time (like a semaphore)
    # this ensures that no two tasks write to the same client stream at the same time
    flushtoken = let
        t = get(streamtokens, clientstream, nothing)
        if t === nothing
            streamtokens[clientstream] = create_flushtoken()
        else
            t
        end
    end

    messagetype = Symbol(parentbody["type"])
    request_id = Symbol(parentbody["request_id"])

    body = parentbody["body"]
    body_parsed = (Symbol(p.first) => withmorebits(p.second) for p in body)

    # Find the requested function and call it with the given arguments:
    result = funkies[messagetype](;body_parsed...)

    # We use a blocking channel (`flushtoken`) to ensure that at most one async process is writing to the HTTP stream at once. HTTP.jl doesn't do this automatically :(
    # This can be improved by using one token per stream, instead of one shared token. (I've tested it and it works.)
    token = take!(flushtoken)
    try
        if clientstream !== nothing
            if isopen(clientstream)
                clientstream.frame_type = HTTP.WebSockets.WS_BINARY
                # Send back the original request_id so that JS can return it to the original (async) request
                write_serialized(clientstream, Dict(
                    :initiator_id => client_id,
                    :request_id => request_id,
                    :body => result,
                ))
            else
                @info "Client $(client_id) stream closed."
                put!(flushtoken, token)
            end
        end
    catch ex
        bt = stacktrace(catch_backtrace())
        if ex isa Base.IOError
            # client socket closed, so we return false (about 5 lines after this one)
        else
            @warn "Failed to write to WebSocket of $(client_id) " exception = (ex, bt)
        end
    end
    put!(flushtoken, token)
end


text:
using MolecularGraph: UndirectedGraph
import Base.convert

function convert(::Type{T}, arr::Array{T,1}) where {T<:Number}
    if size(arr,1) > 1
        return T.(arr)
    end
    return T(arr[1])
end
function convert(::Type{Array{T}}, t::T) where {T<:Number}
    return [t]
end
function convert(::Type{T}, arr::Array{T,1}) where {T<:AbstractString}
    if size(arr,1) > 1
        return T.(arr)
    end
    return T(arr[1])
end
convert(::Type{String}, i::Int) = "$i"
function convert(::Type{String}, f::T) where T<:Union{Float16,Float32,Float64}
	"$f"
end
function varcall(name::String, body::Any)
    name=Symbol(name)
    @eval (($name) = ($body))
	return Symbol(name)
end
function collectall(args...; maxdepth = 5, currentdepth = 1)
	ff = []
	f1 = collect([args...])
	f2 = collect([Base.Iterators.flatten(f1)...])

    for x in f2
		if typeof(x) <: Union{AbstractArray,AbstractRange} && currentdepth < maxdepth
			xx = try
                collect([Base.Iterators.flatten(x)...])
            catch
                collect(x)
            end
			push!(ff,collectall(xx; maxdepth, currentdepth = currentdepth+1)...)
		else
			push!(ff,x)
		end
	end

	typ = typeof(ff[1])
    try
        ff = Vector{typ}(ff)
    catch

    end

	return ff
end
function collectkeys(args)
    return keys(args) |> collect
end
function collectvals(args)
    return values(args) |> collect
end
function reversekv(dict::AbstractDict{K,V}) where {K,V}
	vkdict = [x[2].=>x[1] for x in dict]
    if typeof(dict) <: OrderedDict
        return OrderedDict{V,K}(vkdict)
    end
	return Dict{V,K}(vkdict)
end
pdbS(x; kwargs...) = try
        retrievepdb(x; kwargs...)
    catch
        println("$(x) doesn't work, sorry")
    end
pdbM(x; model = "1", group = "ATOM", kwargs...) = try
        read("$pdbdir\\$(x).pdb", MIToS.PDB.PDBFile, model = "1", group = "ATOM", kwargs...)
    catch
        pdbS(x; kwargs...)
    end
anynan(x) = any(isnan.(x))
indexshift(idxs, shift=1.0) = try
		float.(idxs).+=shift .|> Int
	catch
		float.(idxs).+=shift
end
function tryint(number)
    return (try
        Int64(number)
    catch
        number
    end)
end
function tryfloat(number)
    return (try
        Float64(number)
    catch
        number
    end)
end
function tryfloat32(number)
    return (try
        Float32(number)
    catch
        number
    end)
end
function unitrange(arr::AbstractVector{T}) where {T<:Int}
    start_value = arr[1]
    end_value = arr[end]
    end_value == start_value && error("the start and end points are the same value")
	end_value < start_value && error("the last value is lower than the first")
    for i in 2:size(arr,1)
        if arr[i]-arr[i-1] != 1
            error("inconsistent step for unit range")
        end
    end
    return UnitRange(start_value,end_value)
end
function steprange(arr::AbstractVector{T}) where {T<:Real}
    start_value = arr[1]
    end_value = arr[end]
	start_value == end_value && error("the start and end points are the same value")
    step = 1
    for i in 2:size(arr,1)
        if i == 2
            step = arr[i]-arr[i-1]
        end
        if arr[i]-arr[i-1] != step
            error("inconsistent step for step range")
        end
    end
    return StepRangeLen(start_value,step,length(arr))
end
splatrange(range) = [(range...)]
function splatranges(ranges...)
    splattedrange = []
    for range in ranges
        splattedrange = vcat(splattedrange, splatrange(range))
    end
    return eval([Int64.(splattedrange)...])
end
function internaldistances(atms::AbstractVector{AbstractAtom})
    internaldists = zeros(Float64, (size(atms,1),size(atms,1)))
    for (i,x) in enumerate(atms)
        for (j,y) in enumerate(atms)
            internaldists[i,j] = Distances.euclidean(BioStructures.coords(x),BioStructures.coords(y))
        end
    end
    return internaldists
end
function internaldistances(vals::AbstractArray{Number})
    internaldists = zeros(Float64, (size(vals,1),size(vals,1)))
    for (i,x) in enumerate(vals)
        for (j,y) in enumerate(vals)
            internaldists[i,j] = Distances.euclidean(x,y)
        end
    end
    return internaldists
end
function internaldistances(vals::AbstractArray)
    internaldists = zeros(Float64, (size(vals,1),size(vals,1)))
    for i in 1:size(vals,1)
        for j in 1:size(vals,1)
            internaldists[i,j] = Distances.euclidean(vals[i,:],vals[j,:])
        end
    end
    return internaldists
end
function internaldistorders(internaldists::AbstractArray)
    orders = zeros(Float64, (size(internaldists,1),size(internaldists,2)))
    for i = 1:size(internaldists,1)
        orders[i,:] = sortperm(sortperm(internaldists[i,:]))
    end
    return orders
end
function resmass(res::BioStructures.Residue)
    total = 0.0
    for atm in res
        total+=atomicmasses["$(element(atm, strip=true))"]
    end
    return total
end
function resvdw(res::BioStructures.Residue; scale = 1.0)
    total = 0.0
    for atm in res
        total+=vdwrad["$(element(atm, strip=true))"]
    end
    return total*scale
end
function makeclrgrad(vec::AbstractArray{T}, colrmap::AbstractArray) where T<:Real
    softmaxvec = Flux.softmax(vec)
    scalefactor = size(colrmap,1) / maximum(softmaxvec)
    colorindices = round.(Int64, softmaxvec .* scalefactor)
    indexedcolors = colrmap[colorindices]
    return indexedcolors
end
∑(x) = sum(x)
function centerofpoints(points::AbstractArray{T}) where T <: Number
    xs = points[:,1]
    ys = points[:,2]
    zs = points[:,3]
    return centerofpoints = [ ∑(xs)/size(xs,1), ∑(ys)/size(ys,1), ∑(zs)/size(zs,1) ] |> _g
end
function centerofmass(atms::AbstractArray{T}) where {T<:AbstractAtom}
    masses = []
    positions = coordarray(atms)
    xs = positions[1,:]
    ys = positions[2,:]
    zs = positions[3,:]
    totalmass = 0.0
    for i = 1:length(atms)
        if element(atms[i]) == "C"
            push!(masses, 12.0107)
            totalmass += 12.0107
        elseif element(atms[i]) == "N"
            push!(masses, 14.0067)
            totalmass += 14.0067
        elseif element(atms[i]) == "H"
            push!(masses,1.0079)
            totalmass += 1.0079
        elseif element(atms[i]) == "O"
            push!(masses, 15.9994)
            totalmass += 15.9994
        elseif element(atms[i]) == "S"
            push!(masses, 32.065)
            totalmass += 32.065
        else
            push!(masses, 0.0)
            totalmass += 0.0
        end
    end
    return centerofmass = [ ∑(masses.*xs)/totalmass, ∑(masses.*ys)/totalmass, ∑(masses.*zs)/totalmass ]
end
function surfacearea(coordinates, connectivity)
    totalarea = 0.0
    for i = 1:size(connectivity,1)
        totalarea += area(GeometryBasics.Point3f0.(coordinates[connectivity[i,1],:],
						coordinates[connectivity[i,2],:], coordinates[connectivity[i,3],:]))
    end
    return totalarea
end
function linesegs(arr::AbstractArray{T,3}) where {T<:AbstractFloat}
    new_arr::AbstractArray{Point3f0} = []
    for i in 1:size(arr,1)
        push!(new_arr, Makie.Point3f0(arr[i,1,:]))
        push!(new_arr, Makie.Point3f0(arr[i,2,:]))
    end
    return new_arr |> combinedims |> transpose |> collect
end
transposed(arr::AbstractArray) = transpose(arr) |> collect
_t(arr::AbstractArray) = transposed(arr)
_v(arr::AbstractArray) = reverse(arr; dims = 1)
_h(arr::AbstractArray) = reverse(arr; dims = 2)
function _stripkeys(dict::AbstractDict)
    ks = string.(strip.(keys(dict)))
    return ks
end
_stripallkeys(dicts::AbstractArray) =  _stripkeys.(dicts)
elecolors = Dict( "C" => :gray,
                  "N" => :blue,
                  "H" => :white,
                  "O" => :red,
                  "S" => :yellow,
				  "X" => :gray,
				  "ZN" => :gray,
				  "CL" => :gray
)
aquacolors = Dict("C" => RGBf0(0.5,0.5,0.5),
                  "N" => RGBf0(0.472,0.211,0.499),
                  "H" => RGBf0(0.65,0.96,0.70),
                  "O" => RGBf0(0.111,0.37,0.999),
                  "S" => RGBf0(0.992,0.753,0.525),
				  "X" => RGBf0(0.5,0.5,0.5),
				  "ZN" => RGBf0(0.5,0.5,0.5),
				  "CL" => RGBf0(0.5,0.5,0.5)
)
SMILESaa = OrderedDict("A" => "N[C@@]([H])(C)C(=O)O",
        "R" => "N[C@@]([H])(CCCNC(=N)N)C(=O)O",
        "N" => "N[C@@]([H])(CC(=O)N)C(=O)O",
        "D" => "N[C@@]([H])(CC(=O)O)C(=O)O",
        "C" => "N[C@@]([H])(CS)C(=O)O",
        "Q" => "N[C@@]([H])(CCC(=O)N)C(=O)O",
        "E" => "N[C@@]([H])(CCC(=O)O)C(=O)O",
        "G" => "NCC(=O)O",
        "H" => "N[C@@]([H])(CC1=CN=C-N1)C(=O)O",
        "I" => "N[C@@]([H])([C@]([H])(CC)C)C(=O)O",
        "L" => "N[C@@]([H])(CC(C)C)C(=O)O",
        "K" => "N[C@@]([H])(CCCCN)C(=O)O",
        "M" => "N[C@@]([H])(CCSC)C(=O)O",
        "F" => "N[C@@]([H])(Cc1ccccc1)C(=O)O",
        "P" => "N1[C@@]([H])(CCC1)C(=O)O",
        "S" => "N[C@@]([H])(CO)C(=O)O",
        "T" => "N[C@@]([H])([C@]([H])(O)C)C(=O)O",
        "W" => "N[C@@]([H])(CC(=CN2)C1=C2C=CC=C1)C(=O)O",
        "Y" => "N[C@@]([H])(Cc1ccc(O)cc1)C(=O)O",
        "V" => "N[C@@]([H])(C(C)C)C(=O)O")
#
function protsmiles(aas::T) where {T}
    sm = ""
    if T<:String
        aas = string.([aas...])
    elseif T<:AbstractResidue
        aas = [resletterdict[aas[i].name] for i in 1:size(aas,1)]
    end

    for i in 1:(size(aas,1)-1)
        sm = sm*(SMILESaa[aas[i]][1:end-1])
    end
    sm = sm*(SMILESaa[aas[end]])

    return sm
end
function download_file(url::AbstractString, filename::AbstractString;
                        headers::Dict{String,String}=Dict{String,String}(),
                        kwargs...)
    HTTP.open("GET", url, headers; kwargs...) do stream
        open(filename, "w") do fh
            write(fh, stream)
        end
    end
    return filename
end
function download_file(url::AbstractString;
                        headers::Dict{String,String}=Dict{String,String}(),
                        kwargs...)
    download_file(url, tempname(); headers=headers, kwargs...)
end
function downloadpfam(pfamcode::String; filename::String="$pfamcode.stockholm.gz", kwargs...)
    @assert endswith(filename,".gz") "filename must end with the .gz extension"
    if occursin(r"^PF\d{5}$"i, pfamcode)
        number = pfamcode[3:end]
        download_file("http://pfam.xfam.org/family/PF$(number)/alignment/full/gzipped",
        filename; kwargs...)
    else
        throw(ErrorException("$pfamcode is not a correct Pfam code"))
    end
end
#---
# Functions for Kidera factors, from:
# Kenta Nakai, Akinori Kidera, Minoru Kanehisa, Cluster analysis of amino acid indices for prediction of protein structure and function,
# Protein Engineering, Design and Selection, Volume 2, Issue 2, July 1988, Pages 93–100, https://doi.org/10.1093/protein/2.2.93
kideradict = Dict(
    "A" => [-1.56,-1.67,-0.97,-0.27,-0.93,-0.78,-0.2,-0.08,0.21,-0.48],
    "R" => [0.22,1.27,1.37,1.87,-1.7,0.46,0.92,-0.39,0.23,0.93],
    "N" => [1.14,-0.07,-0.12,0.81,0.18,0.37,-0.09,1.23,1.1,-1.73],
    "D" => [0.58,-0.22,-1.58,0.81,-0.92,0.15,-1.52,0.47,0.76,0.7],
    "C" => [0.12,-0.89,0.45,-1.05,-0.71,2.41,1.52,-0.69,1.13,1.1],
    "Q" => [-0.47,0.24,0.07,1.1,1.1,0.59,0.84,-0.71,-0.03,-2.33],
    "E" => [-1.45,0.19,-1.61,1.17,-1.31,0.4,0.04,0.38,-0.35,-0.12],
    "G" => [1.46,-1.96,-0.23,-0.16,0.1,-0.11,1.32,2.36,-1.66,0.46],
    "H" => [-0.41,0.52,-0.28,0.28,1.61,1.01,-1.85,0.47,1.13,1.63],
    "I" => [-0.73,-0.16,1.79,-0.77,-0.54,0.03,-0.83,0.51,0.66,-1.78],
    "L" => [-1.04,0,-0.24,-1.1,-0.55,-2.05,0.96,-0.76,0.45,0.93],
    "K" => [-0.34,0.82,-0.23,1.7,1.54,-1.62,1.15,-0.08,-0.48,0.6],
    "M" => [-1.4,0.18,-0.42,-0.73,2,1.52,0.26,0.11,-1.27,0.27],
    "F" => [-0.21,0.98,-0.36,-1.43,0.22,-0.81,0.67,1.1,1.71,-0.44],
    "P" => [2.06,-0.33,-1.15,-0.75,0.88,-0.45,0.3,-2.3,0.74,-0.28],
    "S" => [0.81,-1.08,0.16,0.42,-0.21,-0.43,-1.89,-1.15,-0.97,-0.23],
    "T" => [0.26,-0.7,1.21,0.63,-0.1,0.21,0.24,-1.15,-0.56,0.19],
    "W" => [0.3,2.1,-0.72,-1.57,-1.16,0.57,-0.48,-0.4,-2.3,-0.6],
    "Y" => [1.38,1.48,0.8,-0.56,0,-0.68,-0.31,1.03,-0.05,0.53],
    "V" => [-0.74,-0.71,2.04,-0.4,0.5,-0.81,-1.07,0.06,-0.46,0.65],
    "B" => [0.86,-0.145,-0.85,0.81,-0.37,0.26,-0.805,0.85,0.93,-0.515],
    "Z" => [-0.96,0.215,-0.77,1.135,-0.105,0.495,0.44,-0.165,-0.19,-1.225],
    "X" => [0,0,0,0,0,0,0,0,0,0],
    "J" => [-0.885,-0.08,0.775,-0.935,-0.545,-1.01,0.065,-0.125,0.555,-0.425],
    "-" => [0,0,0,0,0,0,0,0,0,0],
    "*" => [0,0,0,0,0,0,0,0,0,0],
    "." => [0,0,0,0,0,0,0,0,0,0],
    "?" => [0,0,0,0,0,0,0,0,0,0],
	"ALA" => [-1.56,-1.67,-0.97,-0.27,-0.93,-0.78,-0.2,-0.08,0.21,-0.48],
    "ARG" => [0.22,1.27,1.37,1.87,-1.7,0.46,0.92,-0.39,0.23,0.93],
    "ASN" => [1.14,-0.07,-0.12,0.81,0.18,0.37,-0.09,1.23,1.1,-1.73],
    "ASP" => [0.58,-0.22,-1.58,0.81,-0.92,0.15,-1.52,0.47,0.76,0.7],
    "CYS" => [0.12,-0.89,0.45,-1.05,-0.71,2.41,1.52,-0.69,1.13,1.1],
    "GLN" => [-0.47,0.24,0.07,1.1,1.1,0.59,0.84,-0.71,-0.03,-2.33],
    "GLU" => [-1.45,0.19,-1.61,1.17,-1.31,0.4,0.04,0.38,-0.35,-0.12],
    "GLY" => [1.46,-1.96,-0.23,-0.16,0.1,-0.11,1.32,2.36,-1.66,0.46],
    "HIS" => [-0.41,0.52,-0.28,0.28,1.61,1.01,-1.85,0.47,1.13,1.63],
    "ILE" => [-0.73,-0.16,1.79,-0.77,-0.54,0.03,-0.83,0.51,0.66,-1.78],
    "LEU" => [-1.04,0,-0.24,-1.1,-0.55,-2.05,0.96,-0.76,0.45,0.93],
    "LYS" => [-0.34,0.82,-0.23,1.7,1.54,-1.62,1.15,-0.08,-0.48,0.6],
    "MET" => [-1.4,0.18,-0.42,-0.73,2,1.52,0.26,0.11,-1.27,0.27],
    "PHE" => [-0.21,0.98,-0.36,-1.43,0.22,-0.81,0.67,1.1,1.71,-0.44],
    "PRO" => [2.06,-0.33,-1.15,-0.75,0.88,-0.45,0.3,-2.3,0.74,-0.28],
    "SER" => [0.81,-1.08,0.16,0.42,-0.21,-0.43,-1.89,-1.15,-0.97,-0.23],
    "THR" => [0.26,-0.7,1.21,0.63,-0.1,0.21,0.24,-1.15,-0.56,0.19],
    "TRP" => [0.3,2.1,-0.72,-1.57,-1.16,0.57,-0.48,-0.4,-2.3,-0.6],
    "TYR" => [1.38,1.48,0.8,-0.56,0,-0.68,-0.31,1.03,-0.05,0.53],
    "VAL" => [-0.74,-0.71,2.04,-0.4,0.5,-0.81,-1.07,0.06,-0.46,0.65],
    "ASX" => [0.86,-0.145,-0.85,0.81,-0.37,0.26,-0.805,0.85,0.93,-0.515],
    "GLX" => [-0.96,0.215,-0.77,1.135,-0.105,0.495,0.44,-0.165,-0.19,-1.225],
    "XAA" => [0,0,0,0,0,0,0,0,0,0],
    "XLE" => [-0.885,-0.08,0.775,-0.935,-0.545,-1.01,0.065,-0.125,0.555,-0.425] )
kdict2 = Dict(
    "A" => [-1.67],
    "R" => [1.27],
    "N" => [-0.07],
    "D" => [-0.22],
    "C" => [-0.89],
    "Q" => [0.24],
    "E" => [0.19],
    "G" => [-1.96],
    "H" => [0.52],
    "I" => [-0.16],
    "L" => [0],
    "K" => [0.82],
    "M" => [0.18],
    "F" => [0.98],
    "P" => [-0.33],
    "S" => [-1.08],
    "T" => [-0.7],
    "W" => [2.1],
    "Y" => [1.48],
    "V" => [-0.71],
    "B" => [-0.145],
    "Z" => [0.215],
    "X" => [0],
    "J" => [-0.08],
    "-" => [0],
    "*" => [0],
    "." => [0],
    "?" => [0],
	"ALA" => [-1.67],
    "ARG" => [1.27],
    "ASN" => [-0.07],
    "ASP" => [-0.22],
    "CYS" => [-0.89],
    "GLN" => [0.24],
    "GLU" => [0.19],
    "GLY" => [-1.96],
    "HIS" => [0.52],
    "ILE" => [-0.16],
    "LEU" => [0],
    "LYS" => [0.82],
    "MET" => [0.18],
    "PHE" => [0.98],
    "PRO" => [-0.33],
    "SER" => [-1.08],
    "THR" => [-0.7],
    "TRP" => [2.1],
    "TYR" => [1.48],
    "VAL" => [-0.71],
    "ASX" => [-0.145],
    "GLX" => [0.215],
    "XAA" => [0],
    "XLE" => [-0.08] )
#
kfweights = [3,3,3,3,2,1,1,1,1,1]
weightedkideradict = OrderedDict(collectkeys(kideradict) .=> [kfweights.*(collectvals(kideradict)[i]) for i in 1:size(collectvals(kideradict),1)])
function _kdict(str::String)
    if length(str) == 3
        kideradict3["$str"]
    elseif length(str) == 1
        kideradict["$str"]
    else
        throw(ErrorException("can't get kdict for $str"))
    end
end
_kdict(c::Char) = _kdict(string(c))
function _wkdict(str::AbstractString)
    if length(str) == 3
        weightedkideradict["$(str |> string)"]
    elseif length(str) == 1
        weightedkideradict["$(str |> string)"]
    elseif length(str) > 3
        return [weightedkideradict[(string(str[i]))] for i in 1:length(str)]
    else
        throw(ErrorException("can't get kdict for $str"))
    end
end
_wkdict(c::Char) = _wkdict(string(c))


text:
# Creating a patch with an optional parameter
@testset "patch with optional parameter" begin
    hourvalue(h::Hour=Hour(0)) = Dates.value(h)

    patch = @patch hourvalue(h::Hour=Hour(21)) = 2 * Dates.value(h)
    apply(patch) do
        @test (@mock hourvalue()) == 42
        @test (@mock hourvalue(Hour(4))) == 8
    end
end

# Creating a patch with a keyword parameter
@testset "patch with keyword parameter" begin
    hourvalue(; hour::Hour=Hour(0)) = Dates.value(hour)

    patch = @patch hourvalue(; hour::Hour=Hour(21)) = 2 * Dates.value(hour)
    apply(patch) do
        @test (@mock hourvalue()) == 42

        # Test @mock calls with keyword arguments
        @test (@mock hourvalue(hour=Hour(4))) == 8      #:kw
        @test (@mock hourvalue(; hour=Hour(4))) == 8    #:parameters
    end
end


text:
using Statistics
using NCDatasets
using Plots
using Flux
using OceanParameterizations
using Oceananigans.Grids
using BSON
using OrdinaryDiffEq, DiffEqSensitivity
include("lesbrary_data.jl")
include("data_containers.jl")
include("animate_prediction.jl")

# Training Datasets
train_files = ["strong_wind", "strong_wind_weak_heating"]
output_gif_directory = "Output"

PATH = pwd()

𝒟train = data(train_files,
                    scale_type=ZeroMeanUnitVarianceScaling,
                    animate=false,
                    animate_dir="$(output_gif_directory)/Training")

uw_NN_model = BSON.load(joinpath(PATH, "NDEs", "vw_NDE_2sims_100.bson"))[:neural_network]
vw_NN_model = BSON.load(joinpath(PATH, "NDEs", "vw_NDE_2sims_100.bson"))[:neural_network]
wT_NN_model = BSON.load(joinpath(PATH, "NDEs", "vw_NDE_2sims_100.bson"))[:neural_network]

# Manually fixing the boundary conditions in the ocean profile
function predict_NDE(NN, x, top, bottom)
    interior = NN(x)
    return [top; interior; bottom]
end

f = 1f-4
H = Float32(abs(𝒟train.uw.z[end] - 𝒟train.uw.z[1]))
τ = Float32(abs(𝒟train.t[:,1][end] - 𝒟train.t[:,1][1]))
Nz = 32
u_scaling = 𝒟train.scalings["u"]
v_scaling = 𝒟train.scalings["v"]
T_scaling = 𝒟train.scalings["T"]
uw_scaling = 𝒟train.scalings["uw"]
vw_scaling = 𝒟train.scalings["vw"]
wT_scaling = 𝒟train.scalings["wT"]
μ_u = Float32(u_scaling.μ)
μ_v = Float32(v_scaling.μ)
σ_u = Float32(u_scaling.σ)
σ_v = Float32(v_scaling.σ)
σ_T = Float32(T_scaling.σ)
σ_uw = Float32(uw_scaling.σ)
σ_vw = Float32(vw_scaling.σ)
σ_wT = Float32(wT_scaling.σ)
uw_weights, re_uw = Flux.destructure(uw_NN_model)
vw_weights, re_vw = Flux.destructure(vw_NN_model)
wT_weights, re_wT = Flux.destructure(wT_NN_model)
uw_top = Float32(𝒟train.uw.scaled[1,1])
uw_bottom = Float32(uw_scaling(-1e-3))
vw_top = Float32(𝒟train.vw.scaled[1,1])
vw_bottom = Float32(𝒟train.vw.scaled[end,1])
wT_top = Float32(𝒟train.wT.scaled[1,1])
wT_bottom₁ = Float32(𝒟train.wT.scaled[end,1])
wT_bottom₂ = Float32(wT_scaling(-4e-8))
size_uw_NN = length(uw_weights)
size_vw_NN = length(vw_weights)
size_wT_NN = length(wT_weights)

p₁ = [f; τ; H; μ_u; μ_v; σ_u; σ_v; σ_T; σ_uw; σ_vw; σ_wT; uw_top; uw_bottom; vw_top; vw_bottom; wT_top; wT_bottom₁; uw_weights; vw_weights; wT_weights]

# Differentiation Operator to take the vertical gradients of fluxes
D_cell = Float32.(Dᶜ(Nz, 1 / Nz))

# Coupled Neural Differential Equations
function NDE_nondimensional_flux(x, p, t)
    f, τ, H, μ_u, μ_v, σ_u, σ_v, σ_T, σ_uw, σ_vw, σ_wT, uw_top, uw_bottom, vw_top, vw_bottom, wT_top, wT_bottom = p[1:17]
    Nz = 32
    uw_weights = p[18:18 + size_uw_NN - 1]
    vw_weights = p[18 + size_uw_NN:18 + size_uw_NN + size_vw_NN - 1]
    wT_weights = p[18 + size_uw_NN + size_vw_NN:18 + size_uw_NN + size_vw_NN + size_wT_NN - 1]
    uw_NN = re_uw(uw_weights)
    vw_NN = re_vw(vw_weights)
    wT_NN = re_wT(wT_weights)
    A = - τ / H
    B = f * τ
    u = x[1:Nz]
    v = x[Nz + 1:2 * Nz]
    T = x[2 * Nz + 1:96]
    dx₁ = A .* σ_uw ./ σ_u .* D_cell * predict_NDE(uw_NN, x, uw_top, uw_bottom) .+ B ./ σ_u .* (σ_v .* v .+ μ_v)
    dx₂ = A .* σ_vw ./ σ_v .* D_cell * predict_NDE(vw_NN, x, vw_top, vw_bottom) .- B ./ σ_v .* (σ_u .* u .+ μ_u)
    dx₃ = A .* σ_wT ./ σ_T .* D_cell * predict_NDE(wT_NN, x, wT_top, wT_bottom)
    return [dx₁; dx₂; dx₃]
end

# Data Subsampling
function time_window(t, uvT, trange)
    return (Float32.(t[trange]), Float32.(uvT[:,trange]))
end

start_index = 1
end_index = 100

timesteps = start_index:5:end_index
uvT₁ = Float32.(𝒟train.uvT_scaled[:,start_index])
uvT₂ = Float32.(𝒟train.uvT_scaled[:,289 + start_index])

t_train, uvT_train₁ = time_window(𝒟train.t, 𝒟train.uvT_scaled, timesteps)
_, uvT_train₂ = time_window(𝒟train.t, 𝒟train.uvT_scaled[:, 290:end], timesteps)
t_train = Float32.(t_train ./ τ)
tspan_train = (t_train[1], t_train[end])

opt_NDE = Tsit5()

prob₁ = ODEProblem(NDE_nondimensional_flux, uvT₁, tspan_train, p₁, saveat=t_train)
prob₂ = ODEProblem(NDE_nondimensional_flux, uvT₂, tspan_train, p₁, saveat=t_train)
sol₁ = solve(prob₁, opt_NDE)
sol₂ = solve(prob₂, opt_NDE)

function loss_NDE_NN()
    p₁ = [f; τ; H; μ_u; μ_v; σ_u; σ_v; σ_T; σ_uw; σ_vw; σ_wT; uw_top; uw_bottom; vw_top; vw_bottom; wT_top; wT_bottom₁; uw_weights; vw_weights; wT_weights]
    p₂ = [f; τ; H; μ_u; μ_v; σ_u; σ_v; σ_T; σ_uw; σ_vw; σ_wT; uw_top; uw_bottom; vw_top; vw_bottom; wT_top; wT_bottom₂; uw_weights; vw_weights; wT_weights]

    _sol₁ = Array(solve(prob₁, opt_NDE, p=p₁, reltol=1f-5, saveat=t_train, sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))
    _sol₂ = Array(solve(prob₂, opt_NDE, p=p₂, reltol=1f-5, saveat=t_train, sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))

    loss = mean(Flux.mse(_sol₁, uvT_train₁) + Flux.mse(_sol₂, uvT_train₂))
    return loss
end

function cb_NDE()
    p₁ = [f; τ; H; μ_u; μ_v; σ_u; σ_v; σ_T; σ_uw; σ_vw; σ_wT; uw_top; uw_bottom; vw_top; vw_bottom; wT_top; wT_bottom₁; uw_weights; vw_weights; wT_weights]
    p₂ = [f; τ; H; μ_u; μ_v; σ_u; σ_v; σ_T; σ_uw; σ_vw; σ_wT; uw_top; uw_bottom; vw_top; vw_bottom; wT_top; wT_bottom₂; uw_weights; vw_weights; wT_weights]

    _sol₁ = Array(solve(prob₁, opt_NDE, p=p₁, reltol=1f-5, saveat=t_train, sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))
    _sol₂ = Array(solve(prob₂, opt_NDE, p=p₂, reltol=1f-5, saveat=t_train, sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))

    loss = mean(Flux.mse(_sol₁, uvT_train₁) + Flux.mse(_sol₂, uvT_train₂))
    @info loss
    return (_sol₁, _sol₂)
end

function save_NDE_weights()
    uw_NN_params = Dict(:weights => uw_weights)
    bson(joinpath(PATH, "Output", "uw_NDE_weights_2DaySuite_2Sims.bson"), uw_NN_params)

    vw_NN_params = Dict(:weights => vw_weights)
    bson(joinpath(PATH, "Output", "vw_NDE_weights_2DaySuite_2Sims.bson"), vw_NN_params)

    wT_NN_params = Dict(:weights => wT_weights)
    bson(joinpath(PATH, "Output", "wT_NDE_weights_2DaySuite_2Sims.bson"), wT_NN_params)
end

function train_NDE(epochs)
    for i in 1:epochs
        @info "epoch $i/$epochs"
        Flux.train!(loss_NDE_NN, Flux.params(uw_weights, vw_weights, wT_weights), Iterators.repeated((), 2), ADAM(0.01), cb=Flux.throttle(cb_NDE, 5))
        if i % 5 == 0
            save_NDE_weights()
        end
    end
    save_NDE_weights()
end

# Train the Neural Differential Equations
train_NDE(2000)


text:
using CorpusLoaders
using Test
using Base.Iterators
using MultiResolutionIterators

@testset "Basic use" begin
    dataset = load(WikiGold())

    @test length(dataset) > 0
    @test minimum(length.(dataset)) > 0
end

@testset "Using flatten_levels" begin
    dataset = load(WikiGold())
    docs = dataset[1:5]

    words = full_consolidate(flatten_levels(docs, (!lvls)(WikiGold, :word)))
    @test length(words) > length(docs)
    @test length(words) > sum(length.(docs))
    @test typeof(words) == Vector{CorpusLoaders.NerOnlyTaggedWord}

    plain_words = word.(words)
    @test typeof(plain_words) <: Vector{String}

    ner_tags = CorpusLoaders.named_entity.(words)
    @test typeof(plain_words) <: Vector{String}
end


text:
# assign index (1..n) to each elements of the set
# returns the tuple:
#   - the vector of elements ordered by their index
#   - the dictionary from element to its index
function _index_elements(all_elms::Set{T}) where T
    ix2elm = collect(all_elms)
    sort!(ix2elm)
    # map from element to its index
    elm2ix = Dict{T,Int}()
    sizehint!(elm2ix, length(ix2elm))
    for i in 1:length(ix2elm)
        elm2ix[ix2elm[i]] = i
    end
    return ix2elm, elm2ix
end

# prepares disjoint tiles for the sets collection
# returns the tuple of :
#   - transposed version of the next sparse mask (set×element)
#   - sparse mask representation of sets collection (element×set)
#   - sparse mask representation of tiles collecttion (element×tile)
#   - sparse mask representation of sets composition using tiles (tile×set)
function _prepare_tiles(sets, elm2ix::Dict{<:Any, Int})
    if isempty(sets)
        # no sets -- no tiles
        return SparseMaskMatrix(0, length(elm2ix)), SparseMaskMatrix(length(elm2ix), 0),
               SparseMaskMatrix(length(elm2ix), 0), SparseMaskMatrix(0, 0)
    end

    # build set-X-element membership matrix
    elmXset = fill(false, length(elm2ix), length(sets))
    for (i, s) in enumerate(sets)
        set_elms = view(elmXset, :, i)
        for e in s
            ix = get(elm2ix, e, 0)
            if ix > 0
                @inbounds set_elms[ix] = true
            end
        end
    end
    setXelm = elmXset'
    # squash element-X-set matrix (duplicate rows) into membership-mask-to-tile dict
    sets2elms = Dict{Vector{Int}, Vector{Int}}()
    sizehint!(sets2elms, length(elm2ix))
    for i in 1:size(setXelm, 2)
        @inbounds set_ixs = findall(view(setXelm, :, i))
        push!(get!(() -> Vector{Int}(), sets2elms, set_ixs), i)
    end

    # build tile-X-set and element-X-tile membership matrices
    set2tile_ixs = [Vector{Int}() for _ in 1:length(sets)]
    tile_elm_ranges = Vector{Int}()
    elm_ixs = Vector{Int}()
    for (set_ixs, tile_elm_ixs) in sets2elms
        push!(tile_elm_ranges, length(elm_ixs)+1)
        tile_ix = length(tile_elm_ranges)
        append!(elm_ixs, tile_elm_ixs)

        for set_ix in set_ixs
            push!(set2tile_ixs[set_ix], tile_ix)
        end
    end
    push!(tile_elm_ranges, length(elm_ixs)+1)

    set_tile_ranges = Vector{Int}()
    tile_ixs = Vector{Int}()
    for (set_ix, set_tile_ixs) in enumerate(set2tile_ixs)
        push!(set_tile_ranges, length(tile_ixs)+1)
        append!(tile_ixs, set_tile_ixs)
    end
    push!(set_tile_ranges, length(tile_ixs)+1)

    return SparseMaskMatrix(setXelm), SparseMaskMatrix(elmXset),
        SparseMaskMatrix(length(elm2ix), length(tile_elm_ranges)-1,
                         tile_elm_ranges, elm_ixs), # elmXtile
        SparseMaskMatrix(length(tile_elm_ranges)-1, length(sets),
                         set_tile_ranges, tile_ixs) # tileXset
end

# calculate the sizes of sets given the sizes of tiles
function _set_sizes(tileXset::SparseMaskMatrix,
                    tile_sizes::AbstractVector{Int})
    @assert length(tile_sizes) == size(tileXset, 1)
    set_sizes = zeros(Int, size(tileXset, 2))
    for six in eachindex(set_sizes)
        tile_ixs = view(tileXset, :, six)
        isempty(tile_ixs) && continue
        @inbounds set_sizes[six] = sum(tix -> tile_sizes[tix], tile_ixs)
    end
    return set_sizes
end

"""
    SetMosaic{T,S}

Represents a collection of (potentially overlapping) sets as
a "mosaic" of non-overlapping "tiles".

# Type parameters
* `T`: type of set elements
* `S`: type of set keys
"""
struct SetMosaic{T,S}
    ix2elm::Vector{T}           # element index to element
    elm2ix::Dict{T, Int}        # element to its index

    ix2set::Vector{S}           # set index to the ID of the most relevant set (if there are duplicates)
    set2ix::Dict{S, Int}        # set ID to index

    set_sizes::Vector{Int}
    set_relevances::Vector{Float64}  # 0..1 set relevance score

    setXelm::SparseMaskMatrix   # set×element membership
    elmXset::SparseMaskMatrix   # element×set membership (transpose of setXelm)
    elmXtile::SparseMaskMatrix  # element-to-tile mask
    tileXset::SparseMaskMatrix  # rows = tiles, cols = sets from collection, true if tile is a subset of a set

    setXset_scores::Matrix{Float64}

    function SetMosaic(setids::AbstractVector{S},
                       sets::AbstractVector{<:Union{AbstractVector{T},Set{T}}},
                       all_elms::Set{T},
                       set_relevances::Union{AbstractVector{<:Real}, Nothing} = nothing;
                       setXset_nextra_elms::Integer=0) where {T, S}
        length(setids) == length(sets) ||
            throw(ArgumentError("Number of set IDs does not match the number of sets"))
        (set_relevances === nothing) || (length(set_relevances) == length(sets)) ||
            throw(ArgumentError("Number of set relevance scores does not match the number of sets"))
        ix2elm, elm2ix = _index_elements(all_elms)
        setXelm, elmXset, elmXtile, tileXset = _prepare_tiles(sets, elm2ix)

        # elementwise sets comparison
        function setels_compare(iset, jset)
            lendiff = length(iset) - length(jset) # prioritize length since we need sorting just to find duplicated
            (lendiff != 0) && return lendiff > 0 ? 1 : -1
            (length(iset) == 0) && return 0 # both empty
            for i in eachindex(iset)
                @inbounds eldiff = iset[i] - jset[i]
                (eldiff != 0) && return eldiff > 0 ? 1 : -1
            end
            return 0
        end
        # sort sets by their tiles and relevance, so that duplicate sets are adjacent
        function sets_isless(i, j)
            els_diff = setels_compare(view(tileXset, :, i), view(tileXset, :, j))
            (els_diff != 0) && return els_diff < 0
            rel_diff = set_relevances === nothing ? 0.0 : set_relevances[i] - set_relevances[j]
            (rel_diff != 0.0) && return rel_diff > 0.0 # descending relevance
            return setids[i] < setids[j]
        end
        @inbounds set_sort = sortperm(1:size(tileXset, 2), lt=sets_isless)
        ix2new = fill(0, length(setids))
        new2ix = Vector{Int}()
        if !isempty(set_sort) # remove duplicate sets keeping the ones with higher relevance and id
            lastset_tiles = view(tileXset, :, 1)
            for setix in set_sort
                curset_tiles = view(tileXset, :, setix)
                if isempty(new2ix) || curset_tiles != lastset_tiles
                    lastset_tiles = curset_tiles
                    push!(new2ix, setix)
                end
                ix2new[setix] = length(new2ix)
            end
            neword = sortperm(new2ix)
            new2ix = new2ix[neword]
            ix2new = invperm(neword)[ix2new]
        end

        tile_sizes = Int[length(view(elmXtile, :, tile_ix)) for tile_ix in 1:size(elmXtile, 2)]
        tileXset = tileXset[:, new2ix]
        set_sizes = _set_sizes(tileXset, tile_sizes)
        new{T, S}(ix2elm, elm2ix, setids[new2ix],
                  Dict(id => ix2new[ix] for (ix, id) in enumerate(setids)),
                  set_sizes, set_relevances === nothing ? fill(1.0, length(new2ix)) : set_relevances[new2ix],
                  setXelm[new2ix, :], elmXset[:, new2ix], elmXtile, tileXset,
                  _setXset_scores(tileXset, length(all_elms) + setXset_nextra_elms, set_sizes, tile_sizes))
    end
end

Base.eltype(::Type{SetMosaic{T, S}}) where {T, S} = T
Base.eltype(mosaic::SetMosaic) = eltype(typeof(mosaic))

setidtype(::Type{SetMosaic{T, S}}) where {T, S} = S
setidtype(mosaic::SetMosaic) = setidtype(typeof(mosaic))

"""
    SetMosaic(sets, [all_elms::Set{T}],
              [set_relevances::AbstractVector{Float64}];
              [setXset_nextra_elms=0]) where {T, S} -> SetMosaic{T, S}

Construct `SetMosaic` for a given sets collection.

# Arguments
  * `sets`: a mapping from the set ids (of type `S`) to its elements (`Set{T}`).
     Could be either a dictionary or a vector (then *S* is *Int*).
  * `setXset_nextra_elms`: added to the number of total elements to adjust
    redundancy scores (so that overlaps of very big terms are still penalized)
"""
SetMosaic(sets::AbstractVector{Set{T}}, all_elms::AbstractSet{T} = foldl(union!, sets, init=Set{T}()),
          set_relevances::Union{AbstractVector{Float64}, Nothing} = nothing;
          kwargs...) where {T} =
    SetMosaic(collect(1:length(sets)), sets, all_elms, set_relevances; kwargs...)

SetMosaic(sets::Dict{S, Set{T}}, all_elms::Set{T} = foldl(union!, values(sets), init=Set{T}()),
          set_relevances::Union{Dict{S, Float64}, Nothing} = nothing;
          kwargs...) where {S, T} =
    SetMosaic(collect(keys(sets)), collect(values(sets)), all_elms,
              set_relevances === nothing ? nothing :
              get.(Ref(set_relevances), keys(sets), 1.0); kwargs...)

"""
    nelements(mosaic::SetMosaic) -> Int

The number of distinct elements (i.e. genes) in the sets collection.
"""
nelements(mosaic::SetMosaic) = length(mosaic.ix2elm)

"""
    nsets(mosaic::SetMosaic) -> Int

The number of tiles (pairwise disjoint sets) in the mosaic representation of the set collection.
"""
ntiles(mosaic::SetMosaic) = size(mosaic.tileXset, 1)

"""
    nsets(mosaic::SetMosaic) -> Int

The number of sets in the collection.
"""
nsets(mosaic::SetMosaic) = size(mosaic.tileXset, 2)

"""
    tile(mosaic::SetMosaic, i::Integer) -> AbstractVector{Int}

Get the `i`-th mosaic tile as the vector of indices of its elements.
"""
tile(mosaic::SetMosaic, tile_ix::Integer) = view(mosaic.elmXtile, :, tile_ix)

"""
    set(mosaic::SetMosaic, i::Integer) -> AbstractVector{Int}

Get the `i`-th set as the vector of indices of its elements.
"""
set(mosaic::SetMosaic, set_ix::Integer) = view(mosaic.elmXset, :, set_ix)

"""
    set(mosaic::SetMosaic, i::Integer) -> Int

Get the size of the `i`-th set.
"""
setsize(mosaic::SetMosaic, set_ix::Integer) = mosaic.set_sizes[set_ix]


text:
function checkSVDDone(S::Vector, thresh::Float64)
  N = length(S)
  (N <= 1 || thresh < 0.0) && return (true, 1)
  S1t = S[1] * thresh
  start = 2
  while start <= N
    (S[start] < S1t) && break
    start += 1
  end
  if start >= N
    return (true, N)
  end
  return (false, start)
end

function svd_recursive(M::AbstractMatrix; thresh::Float64=1E-3, north_pass::Int=2)
  Mr, Mc = size(M)
  if Mr > Mc
    V, S, U = svd_recursive(transpose(M))
    conj!(U)
    conj!(V)
    return U, S, V
  end

  #rho = BLAS.gemm('N','T',-1.0,M,M) #negative to sort eigenvalues greatest to smallest
  rho = -M * M' #negative to sort eigenvalues in decreasing order
  D, U = eigen(Hermitian(rho), 1:size(rho, 1))

  Nd = length(D)

  V = M' * U

  V, R = qr_positive(V)
  for n in 1:Nd
    D[n] = R[n, n]
  end

  (done, start) = checkSVDDone(D, thresh)

  done && return U, D, V

  u = view(U, :, start:Nd)
  v = view(V, :, start:Nd)

  b = u' * (M * v)
  bu, bd, bv = svd_recursive(b; thresh=thresh, north_pass=north_pass)

  u .= u * bu
  v .= v * bv
  view(D, start:Nd) .= bd

  return U, D, V
end

# TODO: maybe move to another location?
# Include options for other svd algorithms
function polar(M::AbstractMatrix)
  U, S, V = svd(M) # calls LinearAlgebra.svd(_)
  return U * V', V * Diagonal(S) * V'
end


text:
module NKModel

#%% cell
#=
# Solving a New Keynesian model with Julia

This file is part of a computational appendix that accompanies the paper.

> MATLAB, Python, Julia: What to Choose in Economics?
>
> Coleman, Lyon, Maliar, and Maliar (2017)

For some details regarding the model solved in this file, please see the
[companion notebook that describes the model](http://bookshelf.quantecon.org/submission/59fa1b45145fc3772b0cef82).

In order to run the codes in this file you will need to install and configure a
few Julia packages. We recommend downloading
[JuliaPro](https://juliacomputing.com/products/juliapro.html) and/or following
the instructions on
[quantecon.org](https://lectures.quantecon.org/py/getting_started.html).

Once your Julia installation is up and running, there are a few additional packages
you will need in order to run the code here. We use the
[InstantiateFromURL](https://github.com/QuantEcon/InstantiateFromURL.jl) pacakge to make
it easy to get all required packages set up in your Julia environment.  If you have not
previously installed the `InstantiateFromURL` package you will need to uncomment the line
`pkg"add InstantiateFromURL"` and run the cell below.
=#

#%% cell
using Pkg
# pkg"add InstantiateFromURL"
using InstantiateFromURL: activate_github_path
activate_github_path("sglyon/CLMMJuliaPythonMatlab", path="NewKeynesian/julia", activate=true, force=true)

#%% cell
#=
## Julia Code

The Julia version of our algorithm is implemented as a few functions defined on
a core type named `Model`. This type is itself composed of three different types
that hold the model parameters, steady state, and grids needed to describe the
numerical model. Before we get to the types, we need to bring in some
dependencies:
=#

#%% cell
# for plotting the output
using PlotlyJS

# for constructing Sobol sequences
using Sobol
# for basis matrix of complete monomials and monomial quadrature rules
using BasisMatrices, QuantEcon
using MAT
using Parameters

# standard library components
using Printf  # printing messages
using Statistics  # computing `mean`
using Random: seed!
using LinearAlgebra: diagm, cholesky, dot
using InteractiveUtils: versioninfo
using DelimitedFiles: writedlm

seed!(42);  # set seed for reproducibility

include("model.jl")
include("solution.jl")

#%% cell
#=
### Running the code

Now that we've done all the hard work to define the model, its solution and
simulation, and accuracy checks, let's put things together and run the code!
=#

#%% cell
function ensurefile(url, localpath)
    if !isfile(localpath)
        println("Downloading $url to $localpath")
        download(url, localpath)
    end
end

function main(m=Model(); io::IO=stdout)
    ensurefile("https://github.com/sglyon/CLMMJuliaPythonMatlab/raw/master/NewKeynesian/Sobol_grids.mat", "Sobol_grids.mat")
    ensurefile("https://github.com/sglyon/CLMMJuliaPythonMatlab/raw/master/NewKeynesian/epsi_test_NK.mat", "epsi_test_NK.mat")
    ensurefile("https://github.com/sglyon/CLMMJuliaPythonMatlab/raw/master/NewKeynesian/random_grids.mat", "random_grids.mat")
    coefs, solve_time = solve(m);

    # simulate the model
    ts = time(); sim = Simulation(m, coefs); simulation_time = time() - ts

    # check accuracy
    tr = time(); resids = Residuals(m, coefs, sim); resids_time = time() - tr

    err_by_eq = max_E(resids)
    l1 = mean(resids)
    l∞ = max(resids)
    tot_time = solve_time + simulation_time + resids_time
    round3(x) = round(x, digits=3)
    round2(x) = round(x, digits=2)

    println(io, "Solver time (in seconds): $(solve_time)")
    println(io, "Simulation time (in seconds): $(simulation_time)")
    println(io, "Residuals time (in seconds): $(resids_time)")
    println(io, "total time (in seconds): $(tot_time)")
    println(io, "\nAPPROXIMATION ERRORS (log10):");
    println(io, "\ta) mean error in the model equations: $(round3(l1))");
    println(io, "\tb) max error in the model equations: $(round3(l∞))");
    println(io, "\tc) max error by equation:$(round3.(err_by_eq))");

    solve_time, simulation_time, resids_time, coefs, sim, resids

    t = 1:100
    p1 = plot([scatter(x=t, y=sim.S[t], name="S"), scatter(x=t, y=sim.F[t], name="F")],
              Layout(title="Figure 1a. S and F"))

    p2 = plot([scatter(x=t, y=sim.Y[t], name="Y"), scatter(x=t, y=sim.Yn[t], name="Yn")],
              Layout(title="Figure 1b. Output and natural output"))

    p3 = plot([scatter(x=t, y=sim.C[t], name="C"), scatter(x=t, y=sim.L[t], name="L")],
              Layout(title="Figure 1c. Consumption and labor"))

    p4 = plot([scatter(x=t, y=sim.δ[t], name="δ"), scatter(x=t, y=sim.R[t], name="R"),
               scatter(x=t, y=sim.π[t], name="π")],
              Layout(title="Figure 1d. Distortion, interest rate and inflation"))

    p = [p1 p2; p3 p4]
    p.plot.layout["width"] = 1000
    p.plot.layout["height"] = 600

    p, solve_time, simulation_time, resids_time, coefs, sim, resids, l1, l∞
end

#%% cell
function build_paper_table()
    # call main once to precompile all routines
    main()

    to_csv = []
    push!(to_csv, ["pi_star", "zlb", "degree", "solve_time", "l_1", "l_inf"])

    this_dir = @__DIR__

    open(joinpath(this_dir, "output.log"), "w") do f
        for params in [Dict(:πstar=>1.0, :zlb=>false),
                       Dict(:πstar=>1.0, :zlb=>true),
                       Dict(:πstar=>1 + 0.0598/4, :zlb=>false)]
            for deg in 1:5
                println("Starting $params and deg=$deg")
                m = Model(;grid_kind=:sobol, deg=deg, params...)

                println(f, "Starting $params and deg=$deg")
                solve_time, l1, linf = main(m, io=f)[[2, 8, 9]]
                push!(to_csv, [params[:πstar], params[:zlb], deg, solve_time, l1, linf])
                @show to_csv
                println(f, "\n"^5)
                flush(f)
            end
        end
        versioninfo(f, verbose=true)
    end

    writedlm(
        joinpath(this_dir, "output_osx.csv"),
        permutedims(hcat(to_csv...), (2, 1)), ','
    )
end

#%% cell
# if "table" in ARGS
#     build_paper_table()
# else
#     results = main()
#     results[1]  # show plot
# end

end  # module

max_stars_repo_name: sglyon/CLMMJuliaPythonMatlab

max_issues_repo_name: sglyon/CLMMJuliaPythonMatlab

max_forks_repo_name: sglyon/CLMMJuliaPythonMatlab


text:
module test_julia_endian

using Test

@test Base.ENDIAN_BOM == 0x04030201 # Little-endian machine

end # module test_julia_endian


text:
using Microclimate, Unitful, Test
using Microclimate: range_interpolator, increment_interpolator, interp_layer, layer_sizes,
                    layer_props, layer_bounds, max_height, lin_interp
using Unitful: W, m, °C, K, s, g, L, kPa

const layerincrements = (0.0, 0.025, 0.05, 0.1, 0.2, 0.3, 0.5, 1.0) .* m
const layerrange = (0.01, 1.20) .* m
const nextlayer = 2.0m

rad = [1000.0, 900.0]W*m^-2
snow = [1.0, 0.9]m
airt = [25.0 24.0
        30.0 29.0]°C .|> K
ws = [1.0 2.0
      3.0 4.0]m*s^-1
rh = [0.7 0.5
      0.8 0.6]
soilt = [20.0 19.0 18.0 17.0 16.0 15.0 14.0 13.0
         12.0 11.0 10.0 9.0  8.0  7.0  6.0  5.0]°C .|> K
soilwp = [-100.0 -90.0 -120.0 -80.0 -110.0 -70.0 -130.0 -60.0
          -100.0 -90.0 -120.0 -80.0 -110.0 -70.0 -130.0 -60.0]kPa
soilwc = [0.2 0.3 0.2 0.3 0.2 0.3 0.2 0.3
          0.4 0.5 0.4 0.5 0.4 0.5 0.4 0.5]

@testset "MicroclimPoint" begin
    env = MicroclimPoint{100,layerincrements,nextlayer,layerrange}(rad, snow, airt, rh, ws, soilt, soilwp, soilwc)
    half = (1.2m - 0.01m) / 2 + 0.01m

    @test all(layer_sizes(env) .≈ (0.0125, 0.025, 0.0375, 0.075, 0.1, 0.15, 0.35, 0.75) .* m)
    @test all(layer_bounds(env) .≈ (0.0125, 0.0375, 0.075, 0.15, 0.25, 0.4, 0.75, 1.5) .* m)
    @test max_height(env) == 1.5m
    @test layer_props(env) == layer_sizes(env) ./ 1.5m

    @test radiation(env, 1, range_interpolator(env, 1.2m)) == 1000.0W*m^-2
    @test snowdepth(env, 1, range_interpolator(env, 1.2m)) == 1.0m

    @testset "ranges interpolate" begin
        @test airtemperature(env, 1, range_interpolator(env, 0.01m)) |> °C == 25.0°C
        @test airtemperature(env, 1, range_interpolator(env, half)) |> °C == 24.5°C
        @test airtemperature(env, 1, range_interpolator(env, 1.2m)) |> °C == 24.0°C
        @test windspeed(env, 2, range_interpolator(env, 1.2m)) == 4.0m*s^-1
        @test relhumidity(env, 2, range_interpolator(env, 1.2m)) == 0.6
    end

    @testset "increments interpolate" begin
        @test soiltemperature(env, 1, increment_interpolator(env, 0.0m)) |> °C == 20.0°C
        @test soiltemperature(env, 1, increment_interpolator(env, 0.075m)) |> °C == 17.5°C
        @test soiltemperature(env, 2, increment_interpolator(env, 0.1m)) |> °C == 9.0°C
        @test soiltemperature(env, 1, increment_interpolator(env, 0.40m)) |> °C == 14.50°C
        @test soiltemperature(env, 1, increment_interpolator(env, 0.5m)) |> °C == 14.0°C
        @test soiltemperature(env, 2, increment_interpolator(env, 1.0m)) |> °C == 5.0°C
        @test soilwaterpotential(env, 1, increment_interpolator(env, 0.0m)) == -100kPa
        @test soilwatercontent(env, 2, increment_interpolator(env, 0.0m)) == 0.4
    end

    @testset "float indices are linearly interpolated in the time dimension" begin
        @test windspeed(env, 1.5, range_interpolator(env, 1.2m)) == 3.0m*s^-1
        @test relhumidity(env, 1.75, range_interpolator(env, 1.2m)) == 0.575
        @test relhumidity(env, 1.5, range_interpolator(env, half)) == 0.65
    end

    @testset "handles heights larger than available" begin
        @test soiltemperature(env, 1, increment_interpolator(env, 2.0m)) == 13.0°C
        @test airtemperature(env, 1, range_interpolator(env, 2.0m)) |> °C == 24.0°C
    end

    @testset "weighted mean" begin
        @test weightedmean(env, soilwatercontent(env), 1, 0.0125m) ≈ 0.2
        mid = sum([0.2, 0.3, 0.2, 0.3] .* (layer_sizes(env)[1:4] ./ 0.15m))
        @test weightedmean(env, soilwatercontent(env), 1, 0.15m) ≈ mid
        max = sum([0.2, 0.3, 0.2, 0.3, 0.2, 0.3, 0.2, 0.3] .* layer_props(env))
        @test weightedmean(env, soilwatercontent(env), 1, 1.5m) ≈ max
        @test weightedmean(env, soilwatercontent(env), 1, 2.0m) ≈ max
        @test mean_soilwatercontent(env, 1, 1.5m) ≈ max
        mean_soilwaterpotential(env, 1, 1.5m)
        mean_soiltemperature(env, 1, 1.5m)
    end

    @testset "max" begin
        @test layermax(soiltemperature(env), 1, increment_interpolator(env, 2.0m)) |> °C == 20.0°C
    end

end

@testset "MicroclimInstant" begin
    point = MicroclimPoint{100,layerincrements,nextlayer,layerrange}(
        rad, snow, airt, rh, ws, soilt, soilwp, soilwc
    )
    instant = MicroclimInstant(point, 1, 1.2m)
    @test radiation(instant) == 1000W*m^-2
    @test snowdepth(instant) == 1.0m

    @test airtemperature(instant) == 24°C
    @test relhumidity(instant) == 0.5
    @test windspeed(instant) == 2m/s

    instant = MicroclimInstant(point, 2, 1.0m)
    @test soiltemperature(instant) |> °C == 5°C
    @test soilwaterpotential(instant) == -60.0kPa
    @test soilwatercontent(instant) == 0.5
end


text:
import Mosek
try
    Pkg.installed("JuMP")
    using JuMP
    using FactCheck
    using Base.Test

    solver = Mosek.MosekSolver(QUIET=false)
    lp_solvers         = Any[solver]
    ip_solvers         = Any[solver]
    quad_solvers       = Any[solver]
    quad_mip_solvers   = Any[solver]
    convex_nlp_solvers = Any[solver]
    sdp_solvers        = Any[solver]
    semi_solvers       = Any[]
    sos_solvers        = Any[]
    lazy_solvers, cut_solvers, heur_solvers, info_solvers = Any[], Any[], Any[], Any[]
    minlp_solvers      = Any[]

    # Solver-dependent tests
    include(joinpath(Pkg.dir("JuMP"),"test","model.jl"))
    include(joinpath(Pkg.dir("JuMP"),"test","probmod.jl"))
    include(joinpath(Pkg.dir("JuMP"),"test","qcqpmodel.jl"))
    include(joinpath(Pkg.dir("JuMP"),"test","nonlinear.jl"))
    include(joinpath(Pkg.dir("JuMP"),"test","sdp.jl"))

    FactCheck.exitstatus()
catch
    # If convex not installed, skip this test
end


text:
function test_ad(test_function, Δoutput, inputs...; atol=1e-6, rtol=1e-6)

    # Verify that the forwards-pass produces the correct answer.
    output, pb = Zygote.pullback(test_function, inputs...)
    @test output ≈ test_function(inputs...)

    # Compute the adjoints using AD and FiniteDifferences.
    dW_ad = pb(Δoutput)
    dW_fd = FiniteDifferences.j′vp(central_fdm(5, 1), test_function, Δoutput, inputs...)

    # Compare AD and FiniteDifferences results.
    @testset "$(typeof(test_function)) argument $n" for n in eachindex(inputs)
        @test dW_ad[n] ≈ dW_fd[n] atol=atol rtol=rtol
    end
end


text:
using QSM: GaussSeidel
using QSM: BackwardSweep, ForwardSweep, SymmetricSweep
using QSM: Poisson2, __norm_residual, residual!

@testset "Gauss-Seidel" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(16, 17, 18), (19, 20, 21)]
        maxit = 2*prod(sz)
        R = CartesianIndices(map(n -> 2:n-1, sz))

        m = zeros(Bool, sz)
        b = zeros(T, sz)
        x0 = zeros(T, sz)

        m[R] .= true
        b[R] .= randn(size(R))
        x0[R] .= randn(size(R))

        A = Poisson2(T, m, Tuple(rand(T, 3)))
        x1 = mul!(copy(x0), A, x0)

        for S in (SymmetricSweep, ForwardSweep, BackwardSweep)
            for RB in (false, true)
                solver = GaussSeidel(S(); iter=maxit, rb=RB)

                y = solver(copy(x0), A, b)
                @test __norm_residual(b, A, y) / norm(b) ≤ sqrt(eps(T))
                @test norm(residual!(copy(y), b, A, y)) / norm(b) ≤ sqrt(eps(T))

                # Fixed-point test.
                # Relaxation should not alter the exact solution.
                # Use exact solution as initial guess.
                solver = GaussSeidel{S, 1, RB}()
                @test solver(copy(x0), A, x1) ≈ x0

                solver = GaussSeidel{S, 25, RB}()
                @test solver(copy(x0), A, x1) ≈ x0
            end
        end
    end
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbdegelev(nurbs::NURBS{I,F},ntimes::Union{I, Vector{I}}
              )::NURBS{I,F} where {I<:Integer,F<:AbstractFloat}

Degree elevates the NURBS curve or surface. This function uses the B-Spline
function [`bspdegelev`](@ref).

# Arguments:
- `nurbs`: NURBS structure (curve or surface see [`nrbmak`](@ref))
- `ntimes`: number or pair of numbers indicating how many times the structure
  shall be degree elevated along U direction or V direction

# Output:
- `enurbs`: new NURBS structure with degrees elevated

# Examples:
```julia
julia> ecrv = nrbdegelev(crv,utimes);
julia> esrf = nrbdegelev(srf,[utimes;vtimes]);
```

Increase the NURBS surface degree twice along the V direction.
```julia
julia> esrf = nrbdegelev(srf,[0;2]);
```

# Reference
- see also: [`bspdegelev`](@ref)
"""
function nrbdegelev(nurbs::NURBS{I,F},ntimes::Union{I, Vector{I}}
                   )::NURBS{I,F} where {I<:Integer,F<:AbstractFloat}
# TODO: this ported version does not do all the tests of the matlab version,
#       instead the type features of julia are used -> check if the matlab
#       errors get caught

degree = nurbs.order .- 1;

if typeof(nurbs)<:NURBS2D
  if typeof(ntimes)<:I || length(ntimes) != 2
    throw(ArgumentError("ntimes must be a integer vector of lenght 2 !"));
  end # if
  # NURBS represents a surface
  dim, num1, num2 = size(nurbs.coefs);
  knots = Vector{Vector{F}}(undef,2);

  # Degree elevate along the v direction
  if ntimes[2] == 0
    coefs = nurbs.coefs;
    knots[2] = nurbs.knots[2];
  else
    coefs = reshape(nurbs.coefs,(4*num1,num2));
    coefs, knots[2] = bspdegelev(degree[2],coefs,nurbs.knots[2],ntimes[2]);
    num2 = size(coefs,2);
    coefs = reshape(coefs,(4, num1, num2));
  end # if

  # Degree elevate along the u direction
  if ntimes[1] == 0
    knots[1] = nurbs.knots[1];
  else
    coefs = permutedims(coefs,[1,3,2]);
    coefs = reshape(coefs,(4*num2,num1));
    coefs, knots[1] = bspdegelev(degree[1],coefs,nurbs.knots[1],ntimes[1]);
    coefs = reshape(coefs,(4,num2,size(coefs,2)));
    coefs = permutedims(coefs,[1,3,2]);
  end # if

  # build the new nurbs surface
  return nrbmak(coefs,knots);
else
  if !(typeof(ntimes)<:I)
    # TODO; well it should be a unsigned integer
    throw(ArgumentError("Valid input of ntimes for curves is a integer!"));
  end # if

  # NURBS represents a curve
  if ntimes == 0
    coefs = nurbs.coefs;
    knots = nurbs.knots;
  else
    coefs, knots = bspdegelev(degree,nurbs.coefs,nurbs.knots,ntimes);
  end # if

  # build the new nurbs curve
  return nrbmak(coefs,knots);
end # if
end # nrbdegelev


text:
if VERSION < v"1.1"
    isnothing(x) = x === nothing
end

struct ADTestFunction
    name::String
    f::Function
    x::Vector
end
struct DistSpec{Tθ<:Tuple, Tx}
    name::Union{Symbol, Expr}
    θ::Tθ
    x::Tx
end

vectorize(v::Number) = [v]
vectorize(v::Diagonal) = v.diag
vectorize(v::Vector{<:Matrix}) = mapreduce(vec, vcat, v)
vectorize(v) = vec(v)
pack(vals...) = reduce(vcat, vectorize.(vals))

@generated function unpack(x, vals...)
    unpacked = []
    ind = :(1)
    for (i, T) in enumerate(vals)
        if T <: Number
            push!(unpacked, :(x[$ind]))
            ind = :($ind + 1)
        elseif T <: Vector{<:Matrix}
            push!(unpacked, :(unpack_vec_of_mats(x[$ind:$ind+sum(length, vals[$i])-1], vals[$i])))
            ind = :($ind + sum(length, vals[$i]))
        elseif T <: Vector
            push!(unpacked, :(x[$ind:$ind+length(vals[$i])-1]))
            ind = :($ind + length(vals[$i]))
        elseif T <: Diagonal
            push!(unpacked, :(Diagonal(x[$ind:$ind+size(vals[$i],1)-1])))
            ind = :($ind + size(vals[$i],1))
        elseif T <: Matrix
            push!(unpacked, :(reshape(x[$ind:($ind+length(vals[$i])-1)], size(vals[$i]))))
            ind = :($ind + length(vals[$i]))
        else
            throw("Unsupported argument")
        end
    end
    return quote
        @assert $ind == length(x) + 1
        return ($(unpacked...),)
    end
end
function unpack_vec_of_mats(x, val)
    ind = 1
    return map(1:length(val)) do i
        out = reshape(x[ind : ind + length(val[i]) - 1], size(val[i]))
        ind += length(val[i])
        out
    end
end

function get_function(dist::DistSpec, inds, val)
    syms = []
    args = []
    for (i, a) in enumerate(dist.θ)
        if i in inds
            sym = gensym()
            push!(syms, sym)
            push!(args, sym)
        else
            push!(args, a)
        end
    end
    if val
        sym = gensym()
        push!(syms, sym)
        expr = quote
            ($(syms...),) -> begin
                temp_args = ($(args...),)
                temp_dist = $(dist.name)(temp_args...)
                temp_x = $(sym)
                if temp_dist isa UnivariateDistribution && temp_x isa AbstractArray
                    temp = logpdf.(temp_dist, temp_x)
                else
                    temp = logpdf(temp_dist, temp_x)
                end
                if temp isa AbstractVector
                    return sum(temp)
                else
                    return temp
                end
            end
        end
        if length(inds) == 0
            f = eval(:(x -> $expr(unpack(x, $(dist.x))...)))
            return ADTestFunction(string(expr), f, pack(dist.x))
        else
            f = eval(:(x -> $expr(unpack(x, $(dist.θ[inds]...), $(dist.x))...)))
            return ADTestFunction(string(expr), f, pack(dist.θ[inds]..., dist.x))
        end
    else
        @assert length(inds) > 0
        expr = quote
            ($(syms...),) -> begin
                temp_args = ($(args...),)
                temp_dist = $(dist.name)(temp_args...)
                temp_x = $(dist.x)
                if temp_dist isa UnivariateDistribution && temp_x isa AbstractArray
                    temp = logpdf.(temp_dist, temp_x)
                else
                    temp = logpdf(temp_dist, temp_x)
                end
                if temp isa AbstractVector
                    return sum(temp)
                else
                    return temp
                end
            end
        end
        expr = :(x -> $expr(unpack(x, $(dist.θ[inds]...))...))
        f = eval(expr)
        return ADTestFunction(string(expr), f, pack(dist.θ[inds]...))
    end
end
function get_all_functions(dist::DistSpec, continuous=false)
    fs = []
    if length(dist.θ) == 0
        push!(fs, get_function(dist, (), true))
    else
        for inds in combinations(1:length(dist.θ))
            push!(fs, get_function(dist, inds, false))
            if continuous
                push!(fs, get_function(dist, inds, true))
            end
        end
    end
    return fs
end

# Taken from Turing.jl
function get_stage()
    if "STAGE" in keys(ENV)
        return ENV["STAGE"]
    else
        return "all"
    end
end

const zygote_counter = Ref(0)

function test_ad(f, at = 0.5; rtol = 1e-8, atol = 1e-8)
    stg = get_stage()
    if stg == "all"
        isarr = isa(at, AbstractArray)
        reverse_tracker = Tracker.data(Tracker.gradient(f, at)[1])
        reverse_zygote = Zygote.gradient(f, at)[1]
        reverse_diff = ReverseDiff.gradient(f, at)
        if isarr
            forward = ForwardDiff.gradient(f, at)
            @test isapprox(reverse_tracker, forward, rtol=rtol, atol=atol)
            @test isapprox(reverse_zygote, forward, rtol=rtol, atol=atol)
            @test isapprox(reverse_diff, forward, rtol=rtol, atol=atol)
        else
            forward = ForwardDiff.derivative(f, at)
            finite_diff = central_fdm(5,1)(f, at)
            @test isapprox(reverse_tracker, forward, rtol=rtol, atol=atol)
            @test isapprox(reverse_tracker, finite_diff, rtol=rtol, atol=atol)
            @test isapprox(reverse_zygote, finite_diff, rtol=rtol, atol=atol)
            @test isapprox(reverse_diff, forward, rtol=rtol, atol=atol)
        end
    elseif stg == "ForwardDiff_Tracker"
        isarr = isa(at, AbstractArray)
        reverse_tracker = Tracker.data(Tracker.gradient(f, at)[1])
        if isarr
            forward = ForwardDiff.gradient(f, at)
            @test isapprox(reverse_tracker, forward, rtol=rtol, atol=atol)
        else
            forward = ForwardDiff.derivative(f, at)
            finite_diff = central_fdm(5,1)(f, at)
            @test isapprox(reverse_tracker, forward, rtol=rtol, atol=atol)
            @test isapprox(reverse_tracker, finite_diff, rtol=rtol, atol=atol)
        end
    elseif stg == "Zygote"
        zygote_counter[] += 1
        isarr = isa(at, AbstractArray)
        if mod(zygote_counter[], 50) == 0
            Zygote.refresh()
        end
        reverse_zygote = Zygote.gradient(f, at)[1]
        forward = ForwardDiff.gradient(f, at)
        @test isapprox(reverse_zygote, forward, rtol=rtol, atol=atol)
    elseif stg == "ReverseDiff"
        reverse_diff = ReverseDiff.gradient(f, at)
        forward = ForwardDiff.gradient(f, at)
        @test isapprox(reverse_diff, forward, rtol=rtol, atol=atol)
    else
        throw("Unsupported test stage.")
    end
end

_to_cov(B) = B * B' + Matrix(I, size(B)...)


text:
using PyCall: PyObject

using PyFortran90Namelists: f90nml, @pyinterface

export Parser, reads

mutable struct Parser
    o::PyObject
end
Parser() = Parser(f90nml.Parser())

Base.read(
    p::Parser,
    nml_fname::AbstractString,
    nml_patch_in = nothing,
    patch_fname = nothing,
) = p.read(nml_fname, nml_patch_in, patch_fname)
reads(p::Parser, nml_string::AbstractString) = p.reads(nml_string)

@pyinterface Parser


text:
using ComparativeAutograder

#setupcode = @quoteandexecute begin
#    struct Fab
#        f::Function
#        a::Float64
#        b::Float64
#    end
#end

tests = Array{TestCase, 1}()

for _ in 1:10
    push!(
        tests,
        TestCase(
            (sin, -randn(), randn())
        )
    )
end

testsuite = TestSuite(
    "calculate100",
    tests,
)

runtestsuite(testsuite)


text:
using LinearAlgebra
using Plots

export plotnormal2d

"""
    plotnormal2d(p, μ::AbstractVector, Σ::AbstractMatrix; ncontours=2, args...)

Plot the contours of a 2d Normal density. `ncontours` is the number of
contour line to plots.
"""
function plotnormal2d(p, μ::AbstractVector, Σ::AbstractMatrix; ncontours=2, color="blue", label="", args...)
    λ, U = eigen(Σ)
    for i in 1:ncontours
        B = U * diagm(i * sqrt.(λ))
        θ = range(0, stop=2 * pi, length=1000)
        circle = hcat(sin.(θ), cos.(θ))'
        contour = B * circle .+ μ
        plot!(p, contour[1, :], contour[2, :]; color=color, label=i > 1 ? "" : label, args...)
    end
    p
end

max_stars_repo_name: fnlandini/BayesianModels

max_issues_repo_name: fnlandini/BayesianModels

max_forks_repo_name: fnlandini/BayesianModels


text:
using FactorGraph
using FactorGraph: deleteat, insert
using Test

@testset "eliminate" begin
    g = demo_fg()
    for i=9:-1:1
        g = eliminate(g, i)
    end
    @test nv(g) == 1
    @test ne(g) == 0
    @test size(g.tensors[]) == ()

    g2 = demo_fg2()
    for i=9:-1:1
        g2 = eliminate(g2, i)
    end
    @test nv(g2) == 1
    @test ne(g2) == 0
    @test size(g2.tensors[]) == (2,2)
end

@testset "deleteat" begin
    a = reshape(1:16, 4,4)
    b = [1  13; 2  14; 3  15; 4  16]
    @test deleteat(a, 2, (2,3)) == b

    b = [1  5   9  13; 3  7  11  15; 4  8  12  16]
    @test deleteat(a, 1, 2) == b
end

@testset "basic simple tn" begin
    g = demo_fg()
    @test issimple_graph(g)
    @test ne(g) == 9
    @test nv(g) == 7
    @test vertices(g) == 1:7
    @test vertices(g, 3) == [2,4]
    @test edges(g, 3) == [2,4,5]
    @test neighbors(g, 3) == [2,5,6]
    @test neighbors(g, (2,3)) |> sort == [1,4,5,6]
    @test occupied_legs(g, 3) |> sort == [1,2,3]
    @test dangling_legs(g, 3) == []
    @test dangling_legs(g) == repeat([Int[]], 7)
    @test occupied_legs(g) .|> sort == [[1], [1,2,3], [1,2,3], [1,2,3], [1,2], [1,2,3,4],[1,2]]
  end

@testset "basic multi-graph tn" begin
    g = demo_fg2()
    @test !issimple_graph(g)
    @test ne(g) == 9
    @test nv(g) == 7
    @test vertices(g) == 1:7
    @test vertices(g, 9) == [7]
    @test edges(g, 3) |> sort == [2,4,5,6]
    @test neighbors(g, 3) == [2,5,6]
    @test occupied_legs(g, 6) |> sort == [1,3,4]
    @test dangling_legs(g, 6) == [2]
    @test dangling_legs(g) == [[2],[],[],[],[],[2],[]]
    @test occupied_legs(g) .|> sort == [[1], [1,2,3,4], [1,2,3,4], [1,2], [1,2], [1,3,4],[1,2]]
    @test filter(ie -> isloop(g, ie), g |> edges) == [9]
    # eliminate a selfloop
    @test dangling_legs(eliminate(g, 9), 7) == []
    g2 = add_vertex(g, randn(2, 2), (9,nothing))
    @test dangling_legs(g2, 8) == [2]
    @test !any(isloop.(Ref(g2), edges(g2)))
end


text:
export  AbstractMeasure,
        AbstractCanonicalMeasure,
        AbstractQuad,
        AbstractOrthoPoly,
        AbstractCanonicalOrthoPoly,
        AbstractTensor

abstract type AbstractQuad{T<:Real} end

abstract type AbstractMeasure end
abstract type AbstractCanonicalMeasure <: AbstractMeasure end

abstract type AbstractOrthoPoly{M<:AbstractMeasure, Q<:AbstractQuad} end
abstract type AbstractCanonicalOrthoPoly{V<:AbstractVector{<:Real}, M, Q} <: AbstractOrthoPoly{M, Q} end

abstract type AbstractTensor{T<:AbstractOrthoPoly} end


text:
# Adapted from the Python astroquery.sdss module, which is released under
# a 3-clause BSD style license.

using HTTPClient: get
import JSON

const TIMEOUT = 60.0

"Base URL for catalog-related queries like SQL and Cross-ID."
const SKYSERVER_BASEURL = "http://skyserver.sdss.org"

"Base URL for downloading data products like spectra and images."
const SAS_BASEURL = "http://data.sdss3.org/sas"

function get_query_url(data_release)
    suffix = (data_release < 11)? "sql.asp": "x_sql.aspx"
    return string(SKYSERVER_BASEURL, "/dr$(data_release)/en/tools/search/",
                  suffix)
end

imaging_url_suffix(base, dr, rerun, run, camcol, band, field) =
    @sprintf("%s/dr%d/boss/photoObj/frames/%d/%d/%d/frame-%s-%06d-%d-%04d.fits.bz2",
             base, dr, rerun, run, camcol, band, run, camcol, field)

"""Remove comments and newlines from SQL statement."""
function sanitize_query(stmt::AbstractString)
    lines = split(stmt, '\n')
    for i in 1:length(lines)
        r = search(lines[i], "--")
        if r != 0:-1
            lines[i] = lines[i][1:(first(r)-1)]
        end
    end
    return join(lines, ' ')
end

"""
query_sdss_sql(sql_query; data_release=12, timeout=60.0)

Query the SDSS database. Returns an IOBuffer with the contents.
Use `bytestring()` to convert this to a string if desired.
"""
function query_sdss_sql(sql_query::AbstractString;
                        data_release=12, timeout=TIMEOUT)

    url = get_query_url(data_release)
    payload = [("cmd", sanitize_query(sql_query)), ("format", "csv")]
    res = get(url; query_params=payload, request_timeout=timeout)
    if res.http_code != 200
        error("request for $url returned status $(res.http_code)")
    end

    return res.body
end


text:
using SphereSurfaceHistogram, GeometryBasics
using Test

@testset "SphereSurfaceHistogram.jl" begin

    # These should pass without issues
    # If they don't the binning fails for these edge cases.
    @testset "Edge Cases" begin
        B = SSHBinner(1_000)
        push!(B, [0., 0., 1.])
        @test true
        push!(B, [0., 0., -1.])
        @test true
        push!(B, [0., 1., 0.])
        @test true
        push!(B, [0., -1., 0.])
        @test true
        push!(B, [1., 0., 0.])
        @test true
        push!(B, [-1., 0., 0.])
        @test true
    end

    @testset "Flat Histogram" begin
        for method in [
            SphereSurfaceHistogram.partition_sphere1,
            SphereSurfaceHistogram.partition_sphere2
        ]
            for _ in 1:10
                N = floor(Int64, 100_000rand()) + 1_000
                B = SSHBinner(N, method=method)

                # bin_positions (should) return center positions of bins in order
                # check if each point maps to the correct bin
                bin_centers = bin_positions(B)
                center_matched_bin = Bool[]
                for (i, v) in enumerate(bin_centers)
                    push!(B, v)
                    push!(center_matched_bin, B.bins[i] == 1)
                end
                @test all(center_matched_bin)
            end
        end
    end

    @testset "Indexing" begin
        B = SSHBinner(10_000)
        append!(B, [[cos(pi), sin(pi), 0.0] for _ in 1:1000])

        # true position is 4969
        @test B.bins[4969] == 1000

        # indexing
        @test B[pi/2, pi] == 1000
        @test sum(B[pi/2, :]) == 1000
        @test sum(B[:, pi]) == 1000
    end
end


text:
## @lc code=start
using LeetCode

function merge_nodes!(head::ListNode)::ListNode
    node, cur = head, head.next
    while !isnothing(cur)
        if cur.val == 0 && !isnothing(cur.next)
            node.next = cur
            node = node.next
        else
            node.val += cur.val
        end
        cur = cur.next
    end
    node.next = nothing
    return head
end
## @lc code=end

## @lc test=start
@testset "2181.merge-nodes-in-between-zeros.jl" begin
    @test merge_nodes!(ListNode{Int}([0,3,1,0,4,5,2,0])) == ListNode{Int}([4, 11])
    @test merge_nodes!(ListNode{Int}([0,1,0,3,0,2,2,0])) == ListNode{Int}([1,3,4])
end
## @lc test=end

## @lc info=start
# link: [solution to 2181](https://leetcode-cn.com/problems/merge-nodes-in-between-zeros/solution/chuang-jian-xin-jie-dian-yuan-di-xiu-gai-dplq/)
## @lc info=end


text:
using Healpix
using HTTP

const halfpi=1.570796326794896619231321691639751442099
const HTMfunc_Pr=3.1415926535897932385E0/180.0

nside = 16
res = Healpix.Resolution(nside)
println(res)

function search_gaia_db(ra, dec, radius)
    θ = halfpi - dec * HTMfunc_Pr
    ϕ = ra * HTMfunc_Pr
    hpx = Healpix.ang2pixRing(res, θ, ϕ)
    println("The pixel index is $hpx")
    return "GAIA DR2 WebQL::<ra:$ra, dec:$dec, radius:$radius> --> hpx:$(hpx-1)"
end

server = "127.0.0.1"
port = 8081
htdocs = "htdocs/"

println("server: ", server, " port: ", port, " htdocs: ", htdocs)

HTTP.serve(server, port) do request::HTTP.Request
    #@show request
    #@show request.method
    #@show HTTP.header(request, "Content-Type")
    #@show HTTP.payload(request)
    println("incoming request for ", request.target)
    try
        if request.target == "/"
            return HTTP.Response(200, read(htdocs * "index.html"))
        else
            file = htdocs * HTTP.unescapeuri(request.target[2:end])
            if isfile(file)
                HTTP.Response(200, read(file))
            else
                if occursin("GAIAWebQL.html?", request.target)
                    pos = findlast(isequal('?'), request.target)
                    params = split(request.target[pos+1:end], "&")
                    println("params:", params)
                    ra = nothing
                    dec = nothing
                    radius = nothing
                    for pair in params
                        if occursin("=", pair)
                            tmp = split(pair, "=")
                            key = tmp[1]
                            value = tryparse(Float64, tmp[2])
                            if key == "ra"
                                ra = value
                            end
                            if key == "dec"
                                dec = value
                            end
                            if key == "radius"
                                radius = value
                            end
                        end
                    end
                    if ra != nothing && dec != nothing && radius != nothing
                        resp = search_gaia_db(ra, dec, radius)
                        return HTTP.Response(resp)
                    else
                        return HTTP.Response(404)
                    end
                else
                    HTTP.Response(404)
                end
            end
        end
    catch e
        return HTTP.Response(404, "Error: $e")
    end
 end

exit()


text:
using MimiFAIRv1_6_2 # imports the model and its components, as referenced by MimiFAIRv1_6_2.component_name
using MimiFUND
using Mimi
using DataFrames
using VegaLite

include(joinpath(@__DIR__, "helper.jl")) # path relative to this file

function get_fundfair162(;ar6_scenario::String = "ssp245",
                            FAIR_first = 1750,
                            FAIR_last = 2300,
                            FUND_first = 1950,
                            FUND_last = 2300)

    FAIR_len = length(FAIR_first:FAIR_last)
    FUND_len = length(FUND_first:FUND_last)

    # start with default MimiFUND model
    m = MimiFUND.get_model()

    # add new dimensions relevant to FAIR (1) minor greenhouse gases and
    # (2) ozone-depleting substances
    other_ghg_names = ["CF4", "C2F6", "C6F14", "HFC23", "HFC32", "HFC43_10", "HFC125", "HFC134a", "HFC143a", "HFC227ea", "HFC245fa", "SF6"]
    ods_names       = ["CFC_11", "CFC_12", "CFC_113", "CFC_114", "CFC_115", "CARB_TET", "MCF", "HCFC_22", "HCFC_141B", "HCFC_142B", "HALON1211", "HALON1202", "HALON1301", "HALON2402", "CH3BR", "CH3CL"]

    set_dimension!(m, :other_ghg, other_ghg_names) # FAIR dims - Set index for Kyoto and ozone ozone-depleting gases.
    set_dimension!(m, :ozone_depleting_substances, ods_names) # FAIR dims - Set index for Kyoto and ozone ozone-depleting gases.

    # set the time dimension of the model to match the full extent of FAIR
    set_dimension!(m, :time, collect(FAIR_first:FAIR_last))

    # add the MimiFAIR components to the MimiFUND model after the emissions
    # component, noting they will take the first and last of the model which is now
    # FAIR_first and FAIR_last

    add_comp!(m, MimiFAIRv1_6_2.ch4_cycle; after = :emissions);
    add_comp!(m, MimiFAIRv1_6_2.n2o_cycle; after = :ch4_cycle);
    add_comp!(m, MimiFAIRv1_6_2.co2_cycle; after = :n2o_cycle);
    add_comp!(m, MimiFAIRv1_6_2.other_ghg_cycles; after = :co2_cycle);
    add_comp!(m, MimiFAIRv1_6_2.o3_depleting_substance_cycles; after = :other_ghg_cycles);
    add_comp!(m, MimiFAIRv1_6_2.co2_forcing; after = :o3_depleting_substance_cycles);
    add_comp!(m, MimiFAIRv1_6_2.ch4_forcing; after = :co2_forcing);
    add_comp!(m, MimiFAIRv1_6_2.n2o_forcing; after = :ch4_forcing);
    add_comp!(m, MimiFAIRv1_6_2.o3_forcing; after = :n2o_forcing);
    add_comp!(m, MimiFAIRv1_6_2.aerosol_direct_forcing; after = :o3_forcing);
    add_comp!(m, MimiFAIRv1_6_2.aerosol_indirect_forcing; after = :aerosol_direct_forcing);
    add_comp!(m, MimiFAIRv1_6_2.other_ghg_forcing; after = :aerosol_indirect_forcing);
    add_comp!(m, MimiFAIRv1_6_2.o3_depleting_substance_forcing; after = :other_ghg_forcing);
    add_comp!(m, MimiFAIRv1_6_2.contrails_forcing; after = :o3_depleting_substance_forcing);
    add_comp!(m, MimiFAIRv1_6_2.bc_snow_forcing; after = :contrails_forcing);
    add_comp!(m, MimiFAIRv1_6_2.landuse_forcing; after = :bc_snow_forcing);
    add_comp!(m, MimiFAIRv1_6_2.total_forcing; after = :landuse_forcing);
    add_comp!(m, MimiFAIRv1_6_2.temperature; after = :total_forcing);

    # set all the FAIR component parameters and make their internal connections
    update_MimiFAIR162_params!(m; start_year = FAIR_first, end_year = FAIR_last, ar6_scenario = ar6_scenario)

    # connect FUND and FAIR

    #
    # FUND Emissions --> Convert Mtons to Gtons --> FAIR CO₂ Cycle
    #

    # destination: FAIR component :co2_cycle parameter :E_CO₂
    # previous source: exogenous parameter
    # new source: FUND :emissions variable :mco2 (which first runs through a multiplier component)

    add_comp!(m, Mimi.multiplier; after = :emissions, first = FUND_first, last = FUND_last);
    set_param!(m, :multiplier, :multiply, fill(1/1000, FAIR_len)) # convert Mtons C coming out of FUND to Gtons C going into FAIR
    connect_param!(m, :multiplier, :input, :emissions, :mco2)

    # Subset AR6 emissions to proper years.
    ar6_emissions_raw = DataFrame(load(joinpath(@__DIR__, "..", "data", "model_data", "AR6_emissions_"*ar6_scenario*"_1750_2300.csv")))
    # ar6_emissions_raw = DataFrame(load(joinpath(@__DIR__, "data", "model_data", "AR6_emissions_"*ar6_scenario*"_1750_2300.csv")))
    emission_indices = indexin(collect(FAIR_first:FAIR_last), ar6_emissions_raw.Year)
    ar6_emissions = ar6_emissions_raw[emission_indices, :]

    FAIR_CO₂_backup = (ar6_emissions.FossilCO2 .+ ar6_emissions.OtherCO2)
    connect_param!(m, :co2_cycle, :E_co2, :multiplier, :output, FAIR_CO₂_backup, backup_offset = 1)

    # A note here is that FAIR v1.6.2 splits land use CO2 emissions and fossil CO2 emission, and feeds
    # their sum into the `co2_cycle`, but just the land use emissions into the `landuse_forcing`
    # component in MimiFAIRv1_6_2. Given FUND does not explicilty provide the land use emissions
    # (AFOLU), we do not modify the latter connection, and thus FAIR will run with its default
    # land use emissions for the `landuse_forcing` component.

    #
    # FAIR Temperature --> FUND CO₂ Cycle, Biodiversity, Ocean and Regional Climate Components
    #

    # destination: FUND components :climateco2cycle, :biodiversity, :ocean --> parameter :temp
    # previous source: FUND component :climatedynamics --> variable :temp
    # new source: FAIR :temperature --> variable :T
    connect_param!(m, :climateco2cycle, :temp, :temperature, :T)
    connect_param!(m, :biodiversity, :temp, :temperature, :T)
    connect_param!(m, :ocean, :temp, :temperature, :T)

    # destination: FUND component :climateregional, parameter :inputtemp
    # previous source: FUND component :climatedynamics --> variable :temp
    # new source: FAIR component :temperature --> variable :T
    connect_param!(m, :climateregional, :inputtemp, :temperature, :T)

    return m

end


text:
import Base: cmp, push!

struct Op
    val     :: Symbol
    unary   :: Bool
end
Op(x::Symbol) = Op(x, false)
Op(x::Vector) = Op(Symbol(x[1]))
Op(x, c::Ctx) = Op(x)

# A symbol can have a different precedence if it acts as a unary operator
Un(x::Symbol) = Op(x, true)
Un(x::Vector) = Un(Symbol(x[1]))
Un(x, c::Ctx) = [nil, Un(x)] # Dummy left operand

l1(x) = length(x) == 1

U(x, c) = Call(x[1], x[2:end]) # Unary
I(x, c) = Call(:at, x)         # Indexing
C(x, c) = Call(nil, x)         # Temporary marker for call

G(x, c) = (!l1(x)
              || first(x) isa Symbol) ?
              Call(x[1], x[2:end]) : x         # Group
B(x, c) = l1(x) ? x : Call(:do, x)             # Block
E(x, c) = l1(x) ? x : infixtoprefix(x, c)      # Infix operators

Na(x, c) = Symbol(x[1])                        # Name
Nu(x, c) = parse(Float64, x[1])                # Number
St(x, c) = unescape_string(x[1][2:end-1])      # String
Ln(x, c) = G(x, c)                             # Line

function combine(x, c)
    isempty(x)     && return nil
    length(x) == 1 && return x[1]
    lhs, rhs = x
    head(rhs) ≡ nil && return Call(lhs, args(rhs))
    insert!(rhs, 2, lhs)
    return rhs
end

struct Ops <: OpPolicy
    prec    :: Order
end
Ops() = Ops(Order(Op))

isop(op,     ops::Ops) = false
isop(op::Op, ops::Ops) = true
node(op, lhs, rhs, ops::Ops) = Call(op, Any[lhs, rhs])

cmp(a::Op, b::Op, c::Ops) = cmp(a, b, c.prec)
push!(ops::Ops, lt::Pair{Op}) = push!(ops.prec, lt)

cleanup(x, c) = (@assert length(x) == 1; cleanup(first(x)))
cleanup(x) = x
cleanup(x::Op) = x.val
cleanup(x::Ex) = Call(cleanup(head(x)), cleanup.(args(x)[args(x) .≢ nil]))


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "OneOrthogonal" begin

##
@info("--------- Testing OneOrthogonal ----------")
using LinearAlgebra, Test
using ACE
using ACE.OrthPolys: OrthPolyBasis, OneOrthogonal
using JuLIP: evaluate, evaluate_d, evaluate!, evaluate_d!
using JuLIP: evaluate
##
# Assemble basis
N = 10  # degree
Nquad = 1000
dt =  2 / Nquad
tdf = range(-1.0+dt/2, 1.0-dt/2, length=Nquad)
ww = ones(Nquad) / Nquad
P = OrthPolyBasis(N,  0, 1.0, 2, -1.0, tdf, ww)
Q = OneOrthogonal(P)

# Evaluate basis
p = Matrix{Float64}(undef,Nquad,N)
q = Matrix{Float64}(undef,Nquad,N)
for i = 1:Nquad
    evaluate!(@view(p[i,:]), nothing, P, tdf[i])
    evaluate!(@view(q[i,:]), nothing, Q, tdf[i])
end

@info("Test spanning property")
@test rank(p'*q) == N

@info("Test orthogonality properties")
@test q'*Diagonal(ww)*q ≈ I
@test norm(q[:,1:end-1]'*Diagonal(ww)*one.(tdf)) < sqrt(eps())

@info("Test derivative")
ε = sqrt(eps())
q1 = evaluate(Q, -ε/2)
q2 = evaluate(Q, ε/2)
dq = evaluate_d(Q,0.0)
@test (q2.-q1)./ε ≈ dq

##

@info("Generate many-body 1-orth basis")

deg = 10  # degree
r0, rcut = 0.5, 2.5
trans = PolyTransform(2, r0)
t0, tcut = trans(r0), trans(rcut)
tdf = range(t0, tcut, length=1000)
ww = ones(Nquad) / Nquad
J = OrthPolyBasis(deg,  0, trans(r0), 2, trans(rcut), tdf, ww)
J1 = OneOrthogonal(J)
P = ACE.TransformedPolys(J1, trans, r0, rcut)

# SHIPBasis(spec::BasisSpec, J; T=Float64, kwargs...)
spec = SparseSHIP(2, deg; wL = 1.5,
                  filterfcn = ν -> ACE.OrthPolys.filter_oneorth(ν, J1))
aceB = SHIPBasis(spec, P)
pureB = ACE.PureBasis(aceB)

I2 = findall(aceB.aalists[1].len .== 1)
B2 = ACE.Utils.SubBasis(aceB, I2)

##

pureB = ACE.PureBasis(aceB, 2)

##

Rs = ACE.rand_vec(aceB.J, 3)
Zs = zeros(Int16, 5)
b = evaluate(aceB, Rs, Zs, 0)
bp = pureB(Rs)
# @assert length(b) == length(aceB) == length(bp)

G, A_dot_1 = let nsamples = 10_000
   G = zeros(length(pureB), length(pureB))
   A_dot_1 = zeros(length(I2))
   for _ = 1:nsamples
      Rs = ACE.rand_vec(aceB.J, 2)
      b = pureB(Rs)
      G += b * b' / nsamples
      A_dot_1 += b[I2] / nsamples
   end
   G, A_dot_1
end

##
D = Diagonal( diag(G).^(-0.5) )
Gsc = D * G * D
@show cond(Gsc[1:9,1:9])  # this shows the A-basis is orthogonal
@show cond(Gsc)           # but not the ∏A basis
@show norm(A_dot_1, Inf)

end


text:
# Practical Machine learning
# Bayesian learning - Naive Bayes example
# Chapter 9

using StatsBase
using NaiveBayes

function test_multinomial()
    print("testing MultinomialNB... ")
    m = MultinomialNB([:a, :b, :c], 5)
    X = [1 2 5 2;
         5 3 -2 1;
         0 2 1 11;
         6 -1 3 3;
         5 7 7 1]
    y = [:a, :b, :a, :c]

    fit(m, X, y)
    @assert predict(m, X) == y
    println("OK")
end

function test_gaussian()
    print("testing GaussianNB... ")
    n_obs = 100
    m = GaussianNB([:a, :b, :c], 5)
    X = randn(5, n_obs)
    y = sample([:a, :b, :c], n_obs)

    fit(m, X, y)
    accuracy = sum(predict(m, X) .== y) / n_obs
    println(accuracy)
    println("OK")
end


text:
import DynamicalBilliards.specular!
@inline function specular!(p::AbstractParticle{T}, r::FiniteWall{T})::Void where {T}
    n = normalvec(r, p.pos)
    φ = atan2(n[2], n[1]) + 0.95(π*rand() - π/2) #this cannot be exactly π/2
    p.vel = SVector{2,T}(cos(φ), sin(φ))
    return nothing
end

import DynamicalBilliards.collisiontime
function collisiontime(p::Particle{T}, w::FiniteWall{T})::T where {T}
    n = normalvec(w, p.pos)
    denom = dot(p.vel, n)
    # case of velocity pointing away of wall:
    denom ≥ 0 && return Inf
    posdot = dot(w.sp-p.pos, n)
    # Case of particle starting behind finite wall:
    posdot ≥ 0 && return Inf
    colt = posdot/denom
    intersection = p.pos .+ colt .* p.vel
    dfc = norm(intersection - w.center)
    if dfc > w.width/2
        return Inf
    else
        return colt
    end
end

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics