text:
function JSServe.jsrender(session::Session, scene::Scene)
    three, canvas = WGLMakie.three_display(session, scene)
    return canvas
end

const WEB_MIMES = (MIME"text/html", MIME"application/vnd.webio.application+html",
                   MIME"application/prs.juno.plotpane+html")
for M in WEB_MIMES
    @eval begin
        function AbstractPlotting.backend_show(::WGLBackend, io::IO, m::$M, scene::Scene)
            three = nothing
            inline_display = JSServe.with_session() do session, request
                three, canvas = three_display(session, scene)
                return canvas
            end
            Base.show(io, m, inline_display)
            return three
        end
    end
end

function scene2image(scene::Scene)
    three = nothing
    session = nothing
    inline_display = JSServe.with_session() do s, request
        session = s
        three, canvas = three_display(s, scene)
        return canvas
    end
    electron_display = display(inline_display)
    task = @async wait(session.js_fully_loaded)
    tstart = time()
    # Jeez... Base.Event was a nice idea for waiting on
    # js to be ready, but if anything fails, it becomes unkillable -.-
    while !istaskdone(task)
        sleep(0.01)
        (time() - tstart > 30) && error("JS Session not ready after 30s waiting")
    end
    # HMMMMPFH... This is annoying - we really need to find a way to have
    # devicePixelRatio work correctly
    img_device_scale = AbstractPlotting.colorbuffer(three)
    return ImageTransformations.imresize(img_device_scale, reverse(size(scene)))
end

function AbstractPlotting.backend_show(::WGLBackend, io::IO, m::MIME"image/png",
                                       scene::Scene)
    img = scene2image(scene)
    return FileIO.save(FileIO.Stream(FileIO.format"PNG", io), img)
end

function AbstractPlotting.backend_show(::WGLBackend, io::IO, m::MIME"image/jpeg",
                                       scene::Scene)
    img = scene2image(scene)
    return FileIO.save(FileIO.Stream(FileIO.format"JPEG", io), img)
end

function AbstractPlotting.backend_showable(::WGLBackend, ::T, scene::Scene) where {T<:MIME}
    return T in WEB_MIMES
end

function session2image(sessionlike)
    s = JSServe.session(sessionlike)
    to_data = js"document.querySelector('canvas').toDataURL()"
    picture_base64 = JSServe.evaljs_value(s, to_data)
    picture_base64 = replace(picture_base64, "data:image/png;base64," => "")
    bytes = JSServe.Base64.base64decode(picture_base64)
    return ImageMagick.load_(bytes)
end

function AbstractPlotting.colorbuffer(screen::ThreeDisplay)
    return session2image(screen)
end

function three_display(session::Session, scene::Scene)
    update!(scene)

    serialized = serialize_scene(scene)
    smaller_serialized = replace_dublicates(serialized)
    # smaller_serialized = [serialized, []]
    JSServe.register_resource!(session, smaller_serialized[1])
    width, height = size(scene)
    canvas = DOM.um("canvas", width=width, height=height)
    comm = Observable(Dict{String,Any}())
    scene_data = Observable(smaller_serialized)
    context = JSObject(session, :context)

    setup = js"""
    function setup([scenes, duplicates]){
        const canvas = $(canvas)
        const renderer = $(WGL).threejs_module(canvas, $comm, $width, $height)
        $(WGL).set_duplicate_references(duplicates)
        const three_scenes = scenes.map($(WGL).deserialize_scene)
        const cam = new $(WGLMakie.THREE).PerspectiveCamera(45, 1, 0, 100)
        console.log(three_scenes)
        $(WGL).start_renderloop(renderer, three_scenes, cam)

        function get_plot(plot_uuid) {
            for (const idx in three_scenes) {
                const plot = three_scenes[idx].getObjectByName(plot_uuid)
                if (plot) {
                    return plot
                }
            }
            return undefined;
        }

        function add_plot(scene, plot) {
            const mesh = $(WGL).deserialize_plot(plot);
        }
        const context = {
            three_scenes,
            add_plot,
            get_plot,
            renderer
        }
        put_on_heap($(uuidstr(context)), context);
    }
    """

    JSServe.onjs(session, scene_data, setup)
    WGLMakie.connect_scene_events!(session, scene, comm)
    WGLMakie.mousedrag(scene, nothing)
    scene_data[] = scene_data[]

    canvas_width = lift(x -> [round.(Int, widths(x))...], pixelarea(scene))
    onjs(session, canvas_width, js"""function update_size(canvas_width){
        const context = $(context);
        const w_h = deserialize_js(canvas_width);
        context.renderer.setSize(w_h[0], w_h[1]);
        var canvas = $(canvas)
        canvas.style.width = w_h[0];
        canvas.style.height = w_h[1];
    }""")
    connect_scene_events!(session, scene, comm)
    mousedrag(scene, nothing)
    three = ThreeDisplay(context)
    push!(scene.current_screens, three)
    return three, canvas
end

function AbstractPlotting.backend_display(::WGLBackend, scene::Scene)
    three = nothing
    inline_display = JSServe.with_session() do s, request
        session = s
        three, canvas = three_display(s, scene)
        return canvas
    end
    electron_display = display(inline_display)
    return three
end


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

@doc raw"""A scope selector represents the AND of the selectors represented by the scoped-resource selector requirements.

    IoK8sApiCoreV1ScopeSelector(;
        matchExpressions=nothing,
    )

    - matchExpressions::Vector{IoK8sApiCoreV1ScopedResourceSelectorRequirement} : A list of scope selector requirements by scope of the resources.
"""
mutable struct IoK8sApiCoreV1ScopeSelector <: SwaggerModel
    matchExpressions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ScopedResourceSelectorRequirement} } # spec name: matchExpressions

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

const _property_map_IoK8sApiCoreV1ScopeSelector = Dict{Symbol,Symbol}(Symbol("matchExpressions")=>Symbol("matchExpressions"))
const _property_types_IoK8sApiCoreV1ScopeSelector = Dict{Symbol,String}(Symbol("matchExpressions")=>"Vector{IoK8sApiCoreV1ScopedResourceSelectorRequirement}")
Base.propertynames(::Type{ IoK8sApiCoreV1ScopeSelector }) = collect(keys(_property_map_IoK8sApiCoreV1ScopeSelector))
Swagger.property_type(::Type{ IoK8sApiCoreV1ScopeSelector }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ScopeSelector[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ScopeSelector }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ScopeSelector[property_name]

function check_required(o::IoK8sApiCoreV1ScopeSelector)
    true
end

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

ext: jl


text:
using Gurobi, Test

@testset "GRB String Attributes" begin
    # enviroment and problem
    env = Gurobi.Env()
    m = Gurobi.Model(env, "test")

    # read mps problem in Gurobi
    read_model(m, joinpath(@__DIR__, "test_get_strarray.mps"))

    # test functions
    @test Gurobi.get_strattrarray(m, "VarName", 1, Gurobi.num_vars(m)) == ["XONE","YTWO","ZTHREE"]
    @test Gurobi.get_strattrarray(m, "ConstrName", 1, Gurobi.num_constrs(m)) == ["LIM1","LIM2","MYEQN"]
end


text:
########################################
# Zulip structures
########################################
mutable struct ZulipMessage
    id::Int
    type::String
    stream_id::Int
    subject::String
    display_recipient::Any
    sender_full_name::String
    sender_id::Int
end
ZulipMessage() = ZulipMessage(-1, "private", -1, "", "", "", -1)
StructTypes.StructType(::Type{ZulipMessage}) = StructTypes.Mutable()

mutable struct ZulipRequest
    data::String
    token::String
    message::ZulipMessage
end
ZulipRequest() = ZulipRequest("", "", ZulipMessage())
StructTypes.StructType(::Type{ZulipRequest}) = StructTypes.Mutable()

########################################
# Zulip Client
########################################
struct ZulipClient
    baseep::String
    ep::String
    headers::Vector{Pair{String, String}}
end

const ZulipOpts = Ref(ZulipClient("", "", []))

function ZulipClient(; email = "", apikey = "", ep = "https://julialang.zulipchat.com", api_version = "v1", use_globally = true)
    if isempty(email) || isempty(apikey)
        throw(ArgumentError("Arguments email and apikey should not be empty."))
    end

    key = base64encode(email * ":" * apikey)
    headers = ["Authorization" => "Basic " * key, "Content-Type" => "application/x-www-form-urlencoded"]
    endpoint = ep * "/api/" * api_version * "/"

    client = ZulipClient(ep, endpoint, headers)
    if use_globally
        ZulipOpts[] = client
    end

    return client
end

function query(client::ZulipClient, apimethod, params; method = "POST")
    params = HTTP.URIs.escapeuri(params)
    url = client.ep * apimethod
    JSON3.read(HTTP.request(method, url, client.headers, params).body)
end

function sendMessage(; client = ZulipOpts[], params...)
    query(client, "messages", Dict(params))
end

function updateMessage(msg_id; client = ZulipOpts[], params...)
    query(client, "messages/" * string(msg_id), Dict(params), method = "PATCH")
end


text:
using StatGeochem
using HDF5

# Load dataset as dict
print("Loading SRTM\n")
@time srtm = get_srtm15plus()

# # Calculate average slope
# print("Calculating slope\n")
# @time slope = ave_slope_earth(srtm["elevation"], srtm["x_lon_cntr"], srtm["y_lat_cntr"], srtm["cellsize"], minmatval=-12000, maxmatval=9000)
#
# # Save results
# fid = h5open("srtm15plus_aveslope.h5","w");
# g = g_create(fid, "vars");
# print("Saving to HDF5:\n")
# g["y_lat_cntr"] = srtm["y_lat_cntr"]
# g["x_lon_cntr"] = srtm["x_lon_cntr"]
# g["cellsize"] = srtm["cellsize"]
# g["scalefactor"] = srtm["scalefactor"]
# @time g["slope","compress",3] = slope;
# close(fid)

# Calculate maximum slope
print("Calculating slope\n")
@time slope = max_slope_earth(srtm["elevation"], srtm["x_lon_cntr"], srtm["y_lat_cntr"], srtm["cellsize"], minmatval=-12000)

# Save results
fid = h5open("srtm15plus_maxslope.h5","w");
g = g_create(fid, "vars");
print("Saving to HDF5:\n")
g["y_lat_cntr"] = srtm["y_lat_cntr"]
g["x_lon_cntr"] = srtm["x_lon_cntr"]
g["cellsize"] = srtm["cellsize"]
g["scalefactor"] = srtm["scalefactor"]
@time g["slope","compress",3] = slope;
close(fid)


text:
using CineFiles
using Test

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


text:
"Where to download cifar from"
cifarurl = "http://www.cs.toronto.edu/~kriz"

"Where to download cifar to"
cifardir = joinpath(@__DIR__,"cifar")

"cifar10() => (xtrn,ytrn,xtst,ytst,labels)"
function cifar10(;
                 tgz="cifar-10-binary.tar.gz",
                 dir="cifar-10-batches-bin",
                 trn=["data_batch_$i.bin" for i=1:5],
                 tst=["test_batch.bin"],
                 lbl="batches.meta.txt",
                 )
    global _cifar10_xtrn, _cifar10_ytrn, _cifar10_xtst, _cifar10_ytst, _cifar10_lbls
    if !(@isdefined _cifar10_xtrn)
        _cifar10_xtrn, _cifar10_ytrn, _cifar10_xtst, _cifar10_ytst, _cifar10_lbls = _cifar_read_tgz(tgz,dir,trn,tst,lbl)
    end
    return _cifar10_xtrn, _cifar10_ytrn, _cifar10_xtst, _cifar10_ytst, _cifar10_lbls
end

"cifar100() => (xtrn,ytrn,xtst,ytst,labels)"
function cifar100(;
                  tgz="cifar-100-binary.tar.gz",
                  dir="cifar-100-binary",
                  trn=["train.bin"],
                  tst=["test.bin"],
                  lbl="fine_label_names.txt",
                  )
    global _cifar100_xtrn, _cifar100_ytrn, _cifar100_xtst, _cifar100_ytst, _cifar100_lbls
    if !(@isdefined _cifar100_xtrn)
        _cifar100_xtrn, _cifar100_ytrn, _cifar100_xtst, _cifar100_ytst, _cifar100_lbls = _cifar_read_tgz(tgz,dir,trn,tst,lbl)
    end
    return _cifar100_xtrn, _cifar100_ytrn, _cifar100_xtst, _cifar100_ytst, _cifar100_lbls
end

"Utility to view a cifar image, requires the Images package"
cifarview(x,i)=colorview(RGB,permutedims(x[:,:,:,i],(3,2,1)))

function _cifar_read_tgz(tgz,dir,trn,tst,labels)
    @info("Reading $tgz...")
    if !isdir(cifardir)
        mkpath(cifardir)
    end
    dirpath = joinpath(cifardir,dir)
    if !isdir(dirpath)
        tgzpath = joinpath(cifardir,tgz)
        if !isfile(tgzpath)
            url = "$cifarurl/$tgz"
            download(url, tgzpath)
        end
        run(`tar xzf $tgzpath -C $cifardir`)
    end
    xtrn,ytrn = _cifar_read_files(dirpath,trn)
    xtst,ytst = _cifar_read_files(dirpath,tst)
    lbls = readlines(joinpath(dirpath,labels))
    return xtrn,ytrn,xtst,ytst,lbls
end

function _cifar_read_files(dir,files)
    xs,ys = [],[]
    for file in files
        x,y = _cifar_read_file(dir,file)
        push!(xs,x); push!(ys,y)
    end
    return cat(xs..., dims=4), vcat(ys...)
end

function _cifar_read_file(dir,file)
    a = read(joinpath(dir,file))
    d = occursin("cifar-100",dir) ? 1 : 0
    a = reshape(a, (3073+d, div(length(a),3073+d)))
    y = a[1+d,:] .+ 0x1 # first row (second for cifar100) is Int8 index representation of correct answers
    x = a[2+d:end,:] # rows 2:end (3:end for cifar100) give 32,32,3,N images
    # y = full(sparse(y,1:length(y),1f0,10,length(y))) # one-hot vector representation
    # maybe convert y to int?
    x = reshape(x ./ 255f0, (32,32,3,div(length(x),3072)))
    return (x,y)
end

nothing


text:
# To ensure a build, it isn't sufficient to modify highs_common.jl.
# You also need to update a line in this file:
#     Last updated: 2022-04-22

include("../highs_common.jl")

script = build_script(shared_libs = "ON")

products = [
    LibraryProduct("libhighs", :libhighs),
    ExecutableProduct("highs", :highs),
]

platforms = expand_cxxstring_abis(platforms)

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version = v"6",
    julia_compat = "1.6",
)

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using Dates

# Change this so e.g. `julia_path("1.1")` gives the right binary for Julia
julia_path(v) = "julia" * replace(v, "." => "p")

println_flush(io, x) = (println(io, x); flush(io))
println_flush(x) = io->println_flush(io, x)

function print_code(cmd)
    return io->begin
        println_flush(io, "\n" * raw"```julia")
        r = run(pipeline(cmd, stdout = io, stderr = io); wait = false)
        wait(r)
        flush(io)
        println_flush(io, raw"```" * "\n")
    end
end

mutable struct LazyBuffer
    buffer::IOBuffer
    fill!
    resolved::Bool
    LazyBuffer(f) = new(IOBuffer(), f, false)
end

LazyBuffer(s::AbstractString) = LazyBuffer(io->println(io, s))

function resolve!(L::LazyBuffer)
    L.resolved && return
    L.fill!(L.buffer)
    L.resolved = true
    return
end

let
    versions = ["1.1", "1.2", "1.3"]
    envs = [("ZygoteMaster", "Zygote master"), ("ZygoteRelease", "Zygote release")]
    methods = readdir("Methods")
    buffers = LazyBuffer[]

    setup_dict = Dict()
    push!(buffers, LazyBuffer("# Attempts to differentiate `eigmin` in Julia\n"))
    push!(buffers, LazyBuffer("Generated by `run.jl`; uses OS-dependent paths that likely need to be updated."))

    # populate buffers
    @info "Lazily defining buffers"  Dates.now()
    for method in methods
        method_path = joinpath("Methods", method)
        push!(buffers, LazyBuffer("## Method $method\n"))
        push!(buffers, LazyBuffer("Code:"))
        push!(buffers, LazyBuffer(print_code(`cat $method_path`)))
        push!(buffers, LazyBuffer("Results:"))

        for (env, envname) in envs
            for version in versions
                binary = julia_path(version)
                push!(buffers, LazyBuffer("\n### Julia $version, with $envname"))
                B = LazyBuffer(print_code(`$binary --project=$env $method_path`))
                push!(buffers, B)
                list = get!(setup_dict, (env, version), [])
                push!(list, B)
            end
        end
    end
    push!(buffers, LazyBuffer("\n## Version information\n"))

    for (env, envname) in envs
        for version in versions
            binary = julia_path(version)
            push!(buffers, LazyBuffer("\n### Julia $version, with $envname\n"))
            push!(buffers, LazyBuffer("Version information:"))
            B = LazyBuffer(print_code(`$binary --project=$env verinfo.jl`))
            push!(buffers, B)
            list = get!(setup_dict, (env, version), [])
            push!(list, B)
        end
    end

    @info "Done defining buffers. Starting to resolve buffers"  Dates.now()

    # Resolve in parallel without overlapping setups
    for (env, envname) in envs
        for version in versions
            binary = julia_path(version)
            @info "Starting to resolve environment" env version Dates.now()
            run(`$binary --project=$env precompile.jl`)
            @info "Done precompile" env version Dates.now()
            list = setup_dict[(env, version)]
            @sync begin
                for B in list
                    @async resolve!(B)
                end
            end
            @info "Finished resolving environment" env version  Dates.now()
        end
    end

    @info "Resolving the remaining buffers" Dates.now()
    resolve!.(buffers) # resolve the rest

    @info "Writing to `README.md`"  Dates.now()
    # write buffers to file in order
    open("README.md", "w") do file
        for L in buffers
            print(file, String(take!(L.buffer)))
        end
    end
    @info "Done!" Dates.now()
end

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
module AtBashM

export atbash

using Main.AffineM: affine
using Main.UtilsM: removeExtras

"""
    atbash(text::String, alphabet::String)

This is secretly an affine cipher with `a` and `b` set to the length of the alphabet - 1, so that's how this function treats it.
"""
function atbash(text::String, alphabet::String)

    return affine(
        removeExtras(text, alphabet),
        alphabet,
        length(alphabet) - 1,
        length(alphabet) - 1,
    )

end # function

end # module

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


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

@doc raw"""
    ObjsBotProfile(;
        app_id=nothing,
        deleted=nothing,
        icons=nothing,
        id=nothing,
        name=nothing,
        team_id=nothing,
        updated=nothing,
    )

    - app_id::DefsAppId
    - deleted::Bool
    - icons::ObjsBotProfileIcons
    - id::DefsBotId
    - name::String
    - team_id::DefsTeam
    - updated::Int32
"""
mutable struct ObjsBotProfile <: SwaggerModel
    app_id::Any # spec type: Union{ Nothing, DefsAppId } # spec name: app_id
    deleted::Any # spec type: Union{ Nothing, Bool } # spec name: deleted
    icons::Any # spec type: Union{ Nothing, ObjsBotProfileIcons } # spec name: icons
    id::Any # spec type: Union{ Nothing, DefsBotId } # spec name: id
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    team_id::Any # spec type: Union{ Nothing, DefsTeam } # spec name: team_id
    updated::Any # spec type: Union{ Nothing, Int32 } # spec name: updated

    function ObjsBotProfile(;app_id=nothing, deleted=nothing, icons=nothing, id=nothing, name=nothing, team_id=nothing, updated=nothing)
        o = new()
        validate_property(ObjsBotProfile, Symbol("app_id"), app_id)
        setfield!(o, Symbol("app_id"), app_id)
        validate_property(ObjsBotProfile, Symbol("deleted"), deleted)
        setfield!(o, Symbol("deleted"), deleted)
        validate_property(ObjsBotProfile, Symbol("icons"), icons)
        setfield!(o, Symbol("icons"), icons)
        validate_property(ObjsBotProfile, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsBotProfile, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ObjsBotProfile, Symbol("team_id"), team_id)
        setfield!(o, Symbol("team_id"), team_id)
        validate_property(ObjsBotProfile, Symbol("updated"), updated)
        setfield!(o, Symbol("updated"), updated)
        o
    end
end # type ObjsBotProfile

const _property_map_ObjsBotProfile = Dict{Symbol,Symbol}(Symbol("app_id")=>Symbol("app_id"), Symbol("deleted")=>Symbol("deleted"), Symbol("icons")=>Symbol("icons"), Symbol("id")=>Symbol("id"), Symbol("name")=>Symbol("name"), Symbol("team_id")=>Symbol("team_id"), Symbol("updated")=>Symbol("updated"))
const _property_types_ObjsBotProfile = Dict{Symbol,String}(Symbol("app_id")=>"DefsAppId", Symbol("deleted")=>"Bool", Symbol("icons")=>"ObjsBotProfileIcons", Symbol("id")=>"DefsBotId", Symbol("name")=>"String", Symbol("team_id")=>"DefsTeam", Symbol("updated")=>"Int32")
Base.propertynames(::Type{ ObjsBotProfile }) = collect(keys(_property_map_ObjsBotProfile))
Swagger.property_type(::Type{ ObjsBotProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsBotProfile[name]))}
Swagger.field_name(::Type{ ObjsBotProfile }, property_name::Symbol) =  _property_map_ObjsBotProfile[property_name]

function check_required(o::ObjsBotProfile)
    (getproperty(o, Symbol("app_id")) === nothing) && (return false)
    (getproperty(o, Symbol("deleted")) === nothing) && (return false)
    (getproperty(o, Symbol("icons")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("team_id")) === nothing) && (return false)
    (getproperty(o, Symbol("updated")) === nothing) && (return false)
    true
end

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


text:
# IJulia hooks for displaying plots with RCall

# TODO: create a special graphics device. This would allow us to not accidently close devices opened by users, and display plots immediately as they appear.

ijulia_mime = nothing

"""
Set options for R plotting with IJulia.

The first argument should be a MIME object: currently supported are
* `MIME("image/png")` [default]
* `MIME("image/svg+xml")`

The remaining arguments (keyword only) are passed to the appropriate R graphics
device: see the relevant R help for details.
"""
function ijulia_setdevice(m::MIME;kwargs...)
    global ijulia_mime
    rcall_p(:options, rcalljl_device=rdevicename(m))
    rcall_p(:options, rcalljl_options=Dict(kwargs))
    ijulia_mime = m
    nothing
end
ijulia_setdevice(m::MIME"image/png") = ijulia_setdevice(m, width=6*72, height=5*72)
ijulia_setdevice(m::MIME"image/svg+xml") = ijulia_setdevice(m, width=6, height=5)

rdevicename(m::MIME"image/png") = :png
rdevicename(m::MIME"image/svg+xml") = :svg

function ijulia_displayfile(m::MIME"image/png", f)
    open(f) do f
        d = read(f)
        display(m,d)
    end
end

function _fix_svg_ids(svg_str::String)
    # R svg images use named defs, which cause problem when used inline, see
    # https://github.com/jupyter/notebook/issues/333
    # we get around this by renaming the elements.
    rand_id = randstring()
    svg_str = replace(svg_str, "id=\"glyph" => "id=\"glyph"*rand_id)
    replace(svg_str, "href=\"#glyph" => "href=\"#glyph"*rand_id)
end

function ijulia_displayfile(m::MIME"image/svg+xml", f)
    open(f) do f
        svg_str = read(f, String)
        display(m, _fix_svg_ids(svg_str))
    end
end

"""
Called after cell evaluation.
Closes graphics device and displays files in notebook.
"""
function ijulia_displayplots()
    if rcopy(Int,rcall_p(Symbol("dev.cur"))) != 1
        rcall_p(Symbol("dev.off"))
        for fn in sort(readdir(ijulia_file_dir))
            ffn = joinpath(ijulia_file_dir,fn)
            ijulia_displayfile(ijulia_mime,ffn)
            rm(ffn)
        end
    end
end

# cleanup after error
function ijulia_cleanup()
    if rcopy(Int,rcall_p(Symbol("dev.cur"))) != 1
        rcall_p(Symbol("dev.off"))
    end
    for fn in readdir(ijulia_file_dir)
        ffn = joinpath(ijulia_file_dir,fn)
        rm(ffn)
    end
end

ijulia_file_dir = ""

function ijulia_init()
    global ijulia_file_dir
    ijulia_file_dir = mktempdir()
    ijulia_file_fmt = joinpath(ijulia_file_dir,"rij_%03d")
    rcall_p(:options,rcalljl_filename=ijulia_file_fmt)

    reval_p(rparse_p("""
        options(device = function(filename=getOption('rcalljl_filename'), ...) {
            args <- c(filename = filename, getOption('rcalljl_options'))
            do.call(getOption('rcalljl_device'), modifyList(args, list(...)))
        })
        """))

    Main.IJulia.push_postexecute_hook(ijulia_displayplots)
    Main.IJulia.push_posterror_hook(ijulia_cleanup)
    ijulia_setdevice(MIME"image/png"())
end


text:
const T_fw_frz_Kelvin = 273.15

const T_fw_frz =  0.0       # Freeze point of freshwater in Celcius
const T_sw_frz = -1.8       # Freeze point of seawater in Celcius

const c_p_sw = 3996.0     # J / kg / C   copied from models/csm_share/shr/shr_const_mod.F90
const ρ_sw   = 1026.0     # kg / m^3     copied from models/csm_share/shr/shr_const_mod.F90
const ρ_si   = 0.917e3    # kg / m^3     copied from models/csm_share/shr/shr_const_mod.F90
const ρ_fw   = 1000.0     # kg / m^3     copied from models/csm_share/shr/shr_const_mod.F90
const g      = 9.80616    # m / s^2      copied from models/csm_share/shr/shr_const_mod.F90
const Re     = 6.37122e6  # m            copied from models/csm_share/shr/shr_const_mod.F90
const Ω      = 2π / (86400 / (1 + 1/365)) # rad / s
const Hf_sw = 3.337e5     # J / kg  latent heat of fusion  copied from models/csm_share/shr/shr_const_mod.F90

const ρcp_sw  = ρ_sw * c_p_sw

const missing_value = 1e20

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
ENV["QSG_RENDER_LOOP"] = "basic"
using CxxWrap # for safe_cfunction
using QML
using GR

qmlfile = joinpath(dirname(Base.source_path()), "qml_ex.qml")

type Parameters
  nbins::Float64
end

parameters = Parameters(30)
w, h = (600, 450)

# Called from QQuickPaintedItem::paint with the QPainterRef as an argument
function paint(p::QML.QPainterRef, item::QML.JuliaPaintedItemRef)
  global w, h

  ENV["GKSwstype"] = 381
  ENV["GKSconid"] = split(repr(p.cpp_object), "@")[2]

  dev = device(p)
  r = effectiveDevicePixelRatio(window(item))
  w, h = width(dev) / r, height(dev) / r

  plt = gcf()
  plt[:size] = (w, h)

  nbins = Int64(round(parameters.nbins))
  hexbin(randn(1000000), randn(1000000),
         nbins=nbins, xlim=(-5,5), ylim=(-5,5), title="nbins: $nbins")

  return
end

function mousePosition(eventx, eventy)
  if w > h
    xn = eventx / w
    yn = (h - eventy) / w
  else
    xn = eventx / h
    yn = (h - eventy) / h
  end
  x, y = ndctowc(xn, yn)
  "($(round(x,4)), $(round(y,4)))"
end

# Convert to cfunction, passing the painter as void*
paint_cfunction = safe_cfunction(paint, Void, (QML.QPainterRef, QML.JuliaPaintedItemRef))

# paint_cfunction becomes a context property
@qmlapp qmlfile paint_cfunction parameters

@qmlfunction mousePosition

exec()


text:
"""
    Spectral{F,M,V1,V2,V3}

Create a spectral layer, whose forward pass is given by

    σ.(B*(λ_in.*x) .- λ_out.*(B*x) .+ b)

where `B` is a static matrix, and `λ_in`, `λ_out` and `b` are trainable vectors of parameters.

These layers can be used as drop-in replacement to traditional dense layers, and are usually as good as dense layers for similar amounts of parameters (but spectral layers are wider for the same number of parameters).

# Examples
```jldoctest
julia> sp = Spectral(5, 2)
Spectral(5 => 2)
```
"""
struct Spectral{F, M, V1, V2, V3}
    B::M
    λ_in::V1
    λ_out::V2
    b::V3
    σ::F
end

function Spectral(in::Integer, out::Integer, σ=identity;
    init=Flux.glorot_uniform,
    bias=true
    )
    B = init(out, in)
    λ_in = init(in)
    λ_out = init(out)
    b = bias ? init(out) : bias
    Spectral(B, λ_in, λ_out, b, σ)
end

Flux.trainable(s::Spectral) = (s.λ_in, s.λ_out, s.b,)

Flux.@functor Spectral

function (s::Spectral)(x::AbstractVecOrMat)
    B, λ_in, λ_out, b, σ = s.B, s.λ_in, s.λ_out, s.b, s.σ
    σ.(B*(λ_in.*x) .- λ_out.*(B*x) .+ b)
end

function Base.show(io::IO, sp::Spectral)
    print(io, "Spectral(", size(sp.B, 2), " => ", size(sp.B, 1))
    sp.σ == identity || print(io, ", ", sp.σ)
    sp.b == false && print(io, "; bias=false")
    print(io, ")")
end

"""
    spectralize(d::Dense)

Create a spectral layer with the same dimensions as `d`.
"""
function spectralize(d::Dense)
    in, out = size(d.weight)
    Spectral(in, out, d.σ)
end


text:
include("types/types.jl")
include("utilities/utl.jl")

function potential_energy(A::AbstractSystem, snap::SNAPParams)
    # Produce NeighborList
    nnlist = neighborlist(A, snap)
    # Get number of coefficients
    num_coeff = get_num_coeffs(snap.twojmax, length(snap.elements), snap.chem_flag)

    # Initialize SNAP Bispectrum, dBispectrum, and Stress arrays
    B = zeros(num_coeff)

    for  (i, ai) in enumerate(A)
        i_element = findall(x->x==Symbol(atomic_symbol(ai)), snap.elements)[1]

        runtime_arrays = initialize_runtime_arrays(i, nnlist, A, snap)
        # These must be done with each new configuration
        compute_ui(i_element, snap, runtime_arrays)
        compute_zi(snap, runtime_arrays)
        compute_bi(snap, runtime_arrays)

        B += runtime_arrays.blist
    end
    return B
end

function force(A::AbstractSystem, snap::SNAPParams)
    number_of_particles = length(A.particles)
    # Produce NeighborList
    nnlist = neighborlist(A, snap)
    println(nnlist)
    # Get number of coefficients
    num_coeff = get_num_coeffs(snap.twojmax, length(snap.elements), snap.chem_flag)

    # Initialize SNAP Bispectrum, dBispectrum, and Stress arrays
    dB = [zeros(num_coeff*length(snap.elements), 3) for i = 1:number_of_particles]
    W = [zeros(num_coeff*length(snap.elements), 6) for i = 1:number_of_particles]

    for  (i, ai) in enumerate(A.particles)
        i_element = findall(x->x==Symbol(ai.element.symbol), snap.elements)[1]

        runtime_arrays = initialize_runtime_arrays(i, nnlist, A, snap)
        # These must be done with each new configuration
        compute_ui(i_element, snap, runtime_arrays)
        compute_zi(snap, runtime_arrays)
        compute_bi(snap, runtime_arrays)

        i_offset = num_coeff * (i_element - 1)
        # Compute Forces and Stresses
        for ind = 1:length(runtime_arrays.indij)
            ij = runtime_arrays.indij[ind]
            ii = ij[1]
            jj = ij[2]
            rij = runtime_arrays.rij[ind]
            wj = runtime_arrays.wj[ind]
            rcut = runtime_arrays.rcutij[ind]
            j_element = findall(x->x==Symbol(A.particles[jj].element.symbol), snap.elements)[1]
            ## Need to zero-out dulist and dblist each time

            # Now compute dulist for (i,j)
            compute_duidrj(rij, wj, rcut, jj, snap, runtime_arrays)
            compute_dbidrj(j_element, snap, runtime_arrays)

            dB[ii][(i_offset+1):(num_coeff+i_offset), :] += runtime_arrays.dblist
            dB[jj][(i_offset+1):(num_coeff+i_offset), :] -= runtime_arrays.dblist

        end
    end
    return dB
end

function virial_stress(A::AbstractSystem, snap::SNAPParams)
    number_of_particles = length(A.particles)
    # Produce NeighborList
    nnlist = neighborlist(A, snap)
    println(nnlist)
    # Get number of coefficients
    num_coeff = get_num_coeffs(snap.twojmax, length(snap.elements), snap.chem_flag)

    # Initialize SNAP Bispectrum, dBispectrum, and Stress arrays
    W = [zeros(num_coeff*length(snap.elements), 6) for i = 1:number_of_particles]

    for  (i, ai) in enumerate(A.particles)
        i_element = findall(x->x==Symbol(ai.element.symbol), snap.elements)[1]

        runtime_arrays = initialize_runtime_arrays(i, nnlist, A, snap)
        # These must be done with each new configuration
        compute_ui(i_element, snap, runtime_arrays)
        compute_zi(snap, runtime_arrays)
        compute_bi(snap, runtime_arrays)

        B[i] = runtime_arrays.blist

        i_offset = num_coeff * (i_element - 1)
        # Compute Forces and Stresses
        for ind = 1:length(runtime_arrays.indij)
            ij = runtime_arrays.indij[ind]
            ii = ij[1]
            jj = ij[2]
            rij = runtime_arrays.rij[ind]
            wj = runtime_arrays.wj[ind]
            rcut = runtime_arrays.rcutij[ind]
            j_element = findall(x->x==Symbol(A.particles[jj].element.symbol), snap.elements)[1]
            ## Need to zero-out dulist and dblist each time

            # Now compute dulist for (i,j)
            compute_duidrj(rij, wj, rcut, jj, snap, runtime_arrays)
            compute_dbidrj(j_element, snap, runtime_arrays)

            dB[ii][(i_offset+1):(num_coeff+i_offset), :] += runtime_arrays.dblist
            dB[jj][(i_offset+1):(num_coeff+i_offset), :] -= runtime_arrays.dblist

            xi, yi, zi = ustrip.(A.particles[ii].position)
            xj, yj, zj = ustrip.(A.particles[jj].position)

            W[ii][(i_offset+1):(num_coeff+i_offset), :] += reshape([runtime_arrays.dblist[:, 1]*xi;
                         runtime_arrays.dblist[:, 2]*yi;
                         runtime_arrays.dblist[:, 3]*zi;
                         runtime_arrays.dblist[:, 2]*zi;
                         runtime_arrays.dblist[:, 1]*zi;
                         runtime_arrays.dblist[:, 1]*yi], num_coeff, 6)
            W[jj][(i_offset+1):(num_coeff+i_offset), :] -= reshape([runtime_arrays.dblist[:, 1]*xj;
                         runtime_arrays.dblist[:, 2]*yj;
                         runtime_arrays.dblist[:, 3]*zj;
                         runtime_arrays.dblist[:, 2]*zj;
                         runtime_arrays.dblist[:, 1]*zj;
                         runtime_arrays.dblist[:, 1]*yj], num_coeff, 6)
        end
    end
    return sum(W)
end

function virial(A::AbstractSystem, snap::SNAPParams)
    v = virial(A, snap)
    return sum(v, dims = 2)
end

function compute_snap(A::AbstractSystem, snap::SNAPParams)
    number_of_particles = length(A.particles)
    # Produce NeighborList
    nnlist = neighborlist(A, snap)
    # Get number of coefficients
    num_coeff = get_num_coeffs(snap.twojmax, length(snap.elements), snap.chem_flag)

    # Initialize SNAP Bispectrum, dBispectrum, and Stress arrays
    B = [zeros(num_coeff) for i = 1:number_of_particles]
    dB = [zeros(num_coeff*length(snap.elements), 3) for i = 1:number_of_particles]
    W = [zeros(num_coeff*length(snap.elements), 6) for i = 1:number_of_particles]

    for  (i, ai) in enumerate(A)
        i_element = findall(x->x==Symbol(ai.atomic_symbol), snap.elements)[1]

        runtime_arrays = initialize_runtime_arrays(i, nnlist, A, snap)
        # These must be done with each new configuration
        compute_ui(i_element, snap, runtime_arrays)
        compute_zi(snap, runtime_arrays)
        compute_bi(snap, runtime_arrays)

        B[i] = runtime_arrays.blist

        i_offset = num_coeff * (i_element - 1)
        # Compute Forces and Stresses
        for ind = 1:length(runtime_arrays.indij)
            ij = runtime_arrays.indij[ind]
            ii = ij[1]
            jj = ij[2]
            rij = runtime_arrays.rij[ind]
            wj = runtime_arrays.wj[ind]
            rcut = runtime_arrays.rcutij[ind]
            j_element = findall(x->x==Symbol(A[jj].atomic_symbol), snap.elements)[1]
            ## Need to zero-out dulist and dblist each time

            # Now compute dulist for (i,j)
            compute_duidrj(rij, wj, rcut, jj, snap, runtime_arrays)
            compute_dbidrj(j_element, snap, runtime_arrays)

            dB[ii][(i_offset+1):(num_coeff+i_offset), :] += runtime_arrays.dblist
            dB[jj][(i_offset+1):(num_coeff+i_offset), :] -= runtime_arrays.dblist

            xi, yi, zi = ustrip.(position(A[ii]))
            xj, yj, zj = ustrip.(position(A[jj]))

            W[ii][(i_offset+1):(num_coeff+i_offset), :] += reshape([runtime_arrays.dblist[:, 1]*xi;
                         runtime_arrays.dblist[:, 2]*yi;
                         runtime_arrays.dblist[:, 3]*zi;
                         runtime_arrays.dblist[:, 2]*zi;
                         runtime_arrays.dblist[:, 1]*zi;
                         runtime_arrays.dblist[:, 1]*yi], num_coeff, 6)
            W[jj][(i_offset+1):(num_coeff+i_offset), :] -= reshape([runtime_arrays.dblist[:, 1]*xj;
                         runtime_arrays.dblist[:, 2]*yj;
                         runtime_arrays.dblist[:, 3]*zj;
                         runtime_arrays.dblist[:, 2]*zj;
                         runtime_arrays.dblist[:, 1]*zj;
                         runtime_arrays.dblist[:, 1]*yj], num_coeff, 6)
        end
    end
    return B, dB, W
end


text:
module TestShow
    using DataTables
    using Compat
    using Base.Test
    import Compat.String
    dt = DataTable(A = 1:3, B = ["x", "y", "z"])

    io = IOBuffer()
    show(io, dt)
    show(io, dt, true)
    showall(io, dt)
    showall(io, dt, true)

    subdt = view(dt, [2, 3]) # dt[dt[:A] .> 1.0, :]
    show(io, subdt)
    show(io, subdt, true)
    showall(io, subdt)
    showall(io, subdt, true)

    dtvec = DataTable[dt for _=1:3]
    show(io, dtvec)
    showall(io, dtvec)

    gd = groupby(dt, :A)
    show(io, gd)
    showall(io, gd)

    dtr = DataTableRow(dt, 1)
    show(io, dtr)

    dt = DataTable(A = Vector{String}(3))

    A = DataTables.StackedVector(Any[[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    show(io, A)
    A = DataTables.RepeatedVector([1, 2, 3], 5, 1)
    show(io, A)
    A = DataTables.RepeatedVector([1, 2, 3], 1, 5)
    show(io, A)

    #Test show output for REPL and similar
    dt = DataTable(Fish = ["Suzy", "Amir"], Mass = [1.5, Nullable()],
                   E = NullableCategoricalArray(["a", Nullable()]))
    io = IOBuffer()
    show(io, dt)
    str = String(take!(io))
    @test str == """
    2×3 DataTables.DataTable
    │ Row │ Fish │ Mass  │ E     │
    ├─────┼──────┼───────┼───────┤
    │ 1   │ Suzy │ 1.5   │ a     │
    │ 2   │ Amir │ #NULL │ #NULL │"""

    # Test computing width for Array{String} columns
    dt = DataTable(Any[["a"]], [:x])
    io = IOBuffer()
    show(io, dt)
    str = String(take!(io))
    @test str == """
    1×1 DataTables.DataTable
    │ Row │ x │
    ├─────┼───┤
    │ 1   │ a │"""
end


text:
#!/usr/bin/env julia

using ArgParse
using Base.Iterators
using Printf

function main()
    args = parse_arguments()
    answers = load_answers()
    prob_strs = select_problems(answers, args["prob_nums"], args["skip"])
    pass_count = 0
    for prob_str in prob_strs
        if test_problem(answers, prob_str)
            pass_count += 1
        end
    end
    println("-"^40)
    println("Solutions passed $(pass_count)/$(length(prob_strs)) tests")
end

"Parses script arguments from the command line and returns them as a `Dict`."
function parse_arguments()
    settings = ArgParseSettings()
    @add_arg_table settings begin
        "-s", "--skip"
            help = "test all *except* the following problems"
            action = :store_true
        "prob_nums"
            help = "number of a specific problem to be tested"
            arg_type = Int
            nargs = '*'
            metavar = "prob_num"
    end
    return parse_args(settings)
end

"Reads problem numbers and answers from `stdin` and returns them as a `Dict`."
function load_answers()
    answers = Dict{String, Int}()
    for line in eachline()
        prob_num, answer = split(line)
        answers[prob_num] = parse(Int, answer)
    end
    return answers
end

"Returns a list of strings representing the problem numbers to be tested."
function select_problems(
    answers::Dict{String, Int},
    prob_nums::Array{Int, 1},
    should_skip::Bool
)
    if isempty(prob_nums)
        prob_strs = collect(keys(answers))
    elseif should_skip
        skip_set = Set(prob_nums)
        prob_strs = String[]
        for prob_str in keys(answers)
            if parse(Int, prob_str) ∉ skip_set
                push!(prob_strs, prob_str)
            end
        end
    else
        prob_strs = map(prob_nums) do prob_num
            @sprintf "%03d" prob_num
        end
    end
    return sort(prob_strs)
end

"Tests a specific problem. Returns `true` on success or `false` on failure."
function test_problem(answers::Dict{String, Int}, prob_str::String)
    if !haskey(answers, prob_str)
        println(stderr, "Problem $prob_str: no answer in input file")
        return false
    end

    print("Testing Problem $prob_str...")
    flush(stdout)

    try_or_fail(prob_str, "no file found for this problem") do
        include(joinpath("..", "src", "Problem$(prob_str).jl"))
    end || return false

    try_or_fail(prob_str, "solution failed with an error", show_error = true) do
        solve_expr = Meta.parse("Problem$(prob_str).solve()")
        (answer, time_sec) = @timed eval(solve_expr)
        correct = answers[prob_str]
        @assert (answer == correct) "expected $correct, but got $answer"
        @printf "PASSED (%6.2f s)\n" time_sec
    end || return false

    return true
end

"Runs `f` or prints an error. Returns `true` on success or `false` on failure."
function try_or_fail(
    f::Function,
    prob_str::String,
    error_msg::String;
    show_error::Bool = false
)
    try
        f()
        return true
    catch e
        println("FAILED")
        println(stderr, "Problem $prob_str: $error_msg")
        if (show_error)
            showerror(stderr, e)
            println(stderr)
            print_trace(stacktrace(catch_backtrace()))
        end
        return false
    end
end

"Prints a stack trace to `stderr` in a human-readable format."
function print_trace(trace::StackTraces.StackTrace)
    for (index, frame) in enumerate(trace)
        println(stderr, " [$index] $frame")
    end
end

main()

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
using ArgMacros

let
    @inlinearguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end
end


text:
"""
    xy06(date1, date2)

X,Y coordinates of celestial intermediate pole from series based
on IAU 2006 precession and IAU 2000A nutation.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- `x`, `y`: CIP X,Y coordinates (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The X,Y coordinates are those of the unit vector towards the
   celestial intermediate pole.  They represent the combined effects
   of frame bias, precession and nutation.

3. The fundamental arguments used are as adopted in IERS Conventions
   (2003) and are from Simon et al. (1994) and Souchay et al.
   (1999).

4. This is an alternative to the angles-based method, via the ERFA
   function [`fw2xy`](@ref) and as used in [`xys06a`](@ref) for example.  The two
   methods agree at the 1 microarcsecond level (at present), a
   negligible amount compared with the intrinsic accuracy of the
   models.  However, it would be unwise to mix the two methods
   (angles-based and series-based) in a single application.

### Called ###

- [`fal03`](@ref): mean anomaly of the Moon
- [`falp03`](@ref): mean anomaly of the Sun
- [`faf03`](@ref): mean argument of the latitude of the Moon
- [`fad03`](@ref): mean elongation of the Moon from the Sun
- [`faom03`](@ref): mean longitude of the Moon's ascending node
- [`fame03`](@ref): mean longitude of Mercury
- [`fave03`](@ref): mean longitude of Venus
- [`fae03`](@ref): mean longitude of Earth
- [`fama03`](@ref): mean longitude of Mars
- [`faju03`](@ref): mean longitude of Jupiter
- [`fasa03`](@ref): mean longitude of Saturn
- [`faur03`](@ref): mean longitude of Uranus
- [`fane03`](@ref): mean longitude of Neptune
- [`fapa03`](@ref): general accumulated precession in longitude

### References ###

- Capitaine, N., Wallace, P.T. & Chapront, J., 2003,
    Astron.Astrophys., 412, 567

- Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

- McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
    IERS Technical Note No. 32, BKG

- Simon, J.L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
    Francou, G. & Laskar, J., Astron.Astrophys., 1994, 282, 663

- Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M., 1999,
    Astron.Astrophys.Supp.Ser. 135, 111

- Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
function xy06(date1, date2)
    x = Ref{Cdouble}()
    y = Ref{Cdouble}()
    ccall((:eraXy06, liberfa), Cvoid,
          (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          date1, date2, x, y)
    return x[], y[]
end

"""
    xys00a(date1, date2)

For a given TT date, compute the X,Y coordinates of the Celestial
Intermediate Pole and the CIO locator s, using the IAU 2000A
precession-nutation model.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- `x`, `y`: Celestial Intermediate Pole (Note 2)
- `s`: The CIO locator s (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The Celestial Intermediate Pole coordinates are the x,y
   components of the unit vector in the Geocentric Celestial
   Reference System.

3. The CIO locator s (in radians) positions the Celestial
   Intermediate Origin on the equator of the CIP.

4. A faster, but slightly less accurate result (about 1 mas for
   X,Y), can be obtained by using instead the [`xys00b`](@ref) function.

### Called ###

- [`pnm00a`](@ref): classical NPB matrix, IAU 2000A
- [`bpn2xy`](@ref): extract CIP X,Y coordinates from NPB matrix
- [`s00`](@ref): the CIO locator s, given X,Y, IAU 2000A

### Reference ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
xys00a

"""
    xys00b(date1, date2)

For a given TT date, compute the X,Y coordinates of the Celestial
Intermediate Pole and the CIO locator s, using the IAU 2000B
precession-nutation model.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- `x`, `y`: Celestial Intermediate Pole (Note 2)
- `s`: The CIO locator s (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The Celestial Intermediate Pole coordinates are the x,y
   components of the unit vector in the Geocentric Celestial
   Reference System.

3. The CIO locator s (in radians) positions the Celestial
   Intermediate Origin on the equator of the CIP.

4. The present function is faster, but slightly less accurate (about
   1 mas in X,Y), than the [`xys00a`](@ref) function.

### Called ###

- [`pnm00b`](@ref): classical NPB matrix, IAU 2000B
- [`bpn2xy`](@ref): extract CIP X,Y coordinates from NPB matrix
- [`s00`](@ref): the CIO locator s, given X,Y, IAU 2000A

### Reference ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
xys00b

"""
    xys06a(date1, date2)

For a given TT date, compute the X,Y coordinates of the Celestial
Intermediate Pole and the CIO locator s, using the IAU 2006
precession and IAU 2000A nutation models.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- `x`, `y`: Celestial Intermediate Pole (Note 2)
- `s`: The CIO locator s (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The Celestial Intermediate Pole coordinates are the x,y components
   of the unit vector in the Geocentric Celestial Reference System.

3. The CIO locator s (in radians) positions the Celestial
   Intermediate Origin on the equator of the CIP.

4. Series-based solutions for generating X and Y are also available:
   see Capitaine & Wallace (2006) and [`xy06`](@ref).

### Called ###

- [`pnm06a`](@ref): classical NPB matrix, IAU 2006/2000A
- [`bpn2xy`](@ref): extract CIP X,Y coordinates from NPB matrix
- [`s06`](@ref): the CIO locator s, given X,Y, IAU 2006

### References ###

- Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

- Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
xys06a

for name in ("xys00a",
             "xys00b",
             "xys06a")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(date1, date2)
            x = Ref{Cdouble}()
            y = Ref{Cdouble}()
            s = Ref{Cdouble}()
            ccall(($fc, liberfa), Cvoid,
                  (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
                  date1, date2, x, y, s)
            return x[], y[], s[]
        end
    end
end


text:
# WARNING this file is only loaded if GBIF.jl is also active
# This all happens thanks to the Requires.jl package

import Base: getindex
import Base: setindex!
import SimpleSDMLayers: clip, latitudes, longitudes, mask!, mask

SimpleSDMLayers.longitudes(record::GBIF.GBIFRecord) = record.longitude
SimpleSDMLayers.latitudes(record::GBIF.GBIFRecord) = record.latitude

"""
    latitudes(records::GBIFRecords)

Returns the non-missing latitudes.
"""
SimpleSDMLayers.latitudes(records::GBIF.GBIFRecords) = filter(!ismissing, [latitudes(record) for record in records])

"""
    longitudes(records::GBIFRecords)

Returns the non-missing longitudes.
"""
SimpleSDMLayers.longitudes(records::GBIF.GBIFRecords) = filter(!ismissing, [longitudes(record) for record in records])

"""
    Base.getindex(p::T, occurrence::GBIF.GBIFRecord) where {T <: SimpleSDMLayer}

Extracts the value of a layer at a given position for a `GBIFRecord`. If the
`GBIFRecord` has no latitude or longitude, this will return `nothing`.
"""
function Base.getindex(layer::T, record::GBIF.GBIFRecord) where {T <: SimpleSDMLayer}
   ismissing(record.latitude) && return nothing
   ismissing(record.longitude) && return nothing
   return layer[Point(record.longitude, record.latitude)]
end

"""
    Base.setindex!(layer::T, v, record::GBIFRecord) where {T <: SimpleSDMResponse}

Changes the values of the cell including the point at the requested latitude and
longitude. **Be careful**, this function will not update a cell that has
`nothing`.
"""
function Base.setindex!(layer::SimpleSDMResponse{T}, v::T, record::GBIF.GBIFRecord) where {T}
   ismissing(record.latitude) && return nothing
   ismissing(record.longitude) && return nothing
   isnothing(layer[record]) && return nothing
   setindex!(layer, v, Point(record.longitude, record.latitude))
end

"""
    clip(layer::T, records::GBIF.GBIFRecords)

Returns a clipped version (with a 10% margin) around all occurences in a
GBIFRecords collection.
"""
function SimpleSDMLayers.clip(layer::T, records::GBIF.GBIFRecords) where {T <: SimpleSDMLayer}
   occ_latitudes = latitudes(records)
   occ_longitudes = longitudes(records)

   lat_min, lat_max = extrema(occ_latitudes)
   lon_min, lon_max = extrema(occ_longitudes)

   lat_Δ = abs(lat_max - lat_min)
   lon_Δ = abs(lon_max - lon_min)

   scaling = 0.1
   lon_s = scaling*lon_Δ
   lat_s = scaling*lat_Δ

   lat_max = min(layer.top, lat_max+lat_s)
   lat_min = max(layer.bottom, lat_min-lat_s)
   lon_max = min(layer.right, lon_max+lon_s)
   lon_min = max(layer.left, lon_min-lon_s)

   return clip(layer; left=lon_min, right=lon_max, bottom=lat_min, top=lat_max)
end

"""
    Base.getindex(layer::T, records::GBIF.GBIFRecords) where {T <: SimpleSDMLayer}

Returns the values of a layer at all occurrences in a `GBIFRecords` collection.
"""
function Base.getindex(layer::T, records::GBIF.GBIFRecords) where {T <: SimpleSDMLayer}
   return convert(Vector{eltype(layer)}, filter(!isnothing, [layer[records[i]] for i in 1:length(records)]))
end

"""
    Base.getindex(layer::T, records::Vector{GBIF.GBIFRecord}) where {T <: SimpleSDMLayer}

Returns the values of a layer at all occurrences in a `GBIFRecord` array.
"""
function Base.getindex(layer::T, records::Vector{GBIF.GBIFRecord}) where {T <: SimpleSDMLayer}
   return [layer[record] for record in records]
end

"""
    mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: AbstractBool}

Fills a layer (most likely created with `similar`) so that the values are `true`
if an occurrence is found in the cell, `false` if not.
"""
function mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: Bool}
    for record in records
        if !isnothing(layer[record])
            layer[record] = true
        end
    end
    return layer
end

"""
    mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: Number}

Fills a layer (most likely created with `similar`) so that the values reflect
the number of occurrences in the cell.
"""
function mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: Number}
    for record in records
        if !isnothing(layer[record])
            layer[record] = layer[record] + one(T)
        end
    end
    return layer
end

"""
    mask(layer::SimpleSDMLayer, records::GBIF.GBIFRecords, element_type::Type=Bool)

Create a new layer storing information about the presence of occurrences in the
cells, either counting (numeric types) or presence-absence-ing (boolean types)
them.
"""
function mask(layer::SimpleSDMLayer, records::GBIF.GBIFRecords, element_type::Type=Bool)
    returnlayer = similar(layer, element_type)
    mask!(returnlayer, records)
    return returnlayer
end


text:
# ECOS orders differently than MOI the second and third dimension of the
# exponential cone

struct PermutedExponentialCone <: MOI.AbstractVectorSet end

Base.copy(set::PermutedExponentialCone) = set

struct PermutedExponentialBridge{T,F} <: MOI.Bridges.Constraint.SetMapBridge{
    T,
    PermutedExponentialCone,
    MOI.ExponentialCone,
    F,
    F,
}
    constraint::MOI.ConstraintIndex{F,PermutedExponentialCone}
end

function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{<:PermutedExponentialBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.ExponentialCone},
) where {T}
    return PermutedExponentialBridge{T,F}
end

function MOI.Bridges.map_set(
    ::Type{<:PermutedExponentialBridge},
    ::MOI.ExponentialCone,
)
    return PermutedExponentialCone()
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:PermutedExponentialBridge},
    ::PermutedExponentialCone,
)
    return MOI.ExponentialCone()
end

function MOI.Bridges.map_function(
    ::Type{<:PermutedExponentialBridge{T}},
    func,
) where {T}
    scalars = MOI.Utilities.eachscalar(func)
    return scalars[[1, 3, 2]]
end

# Swapping indices 2 <-> 3 is an involution (it is its own inverse)

function MOI.Bridges.inverse_map_function(
    BT::Type{<:PermutedExponentialBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end

# It is also symmetric (it is its own transpose)

function MOI.Bridges.adjoint_map_function(
    BT::Type{<:PermutedExponentialBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end

function MOI.Bridges.inverse_adjoint_map_function(
    BT::Type{<:PermutedExponentialBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end


text:
"""
    Chain <: Imputor

Runs multiple `Imputor`s on the same data in the order they're provided.

# Fields
* `imputors::Array{Imputor}`
"""
struct Chain <: Imputor
    imputors::Vector{Imputor}
end

"""
    Chain(imputors::Imputor...) -> Chain

Creates a Chain using the `Imputor`s provided (ordering matters).
"""
Chain(imputors::Imputor...) = Chain(collect(imputors))

"""
    impute!(imp::Chain, missing::Function, data::Dataset; limit::Float64=0.1)

Creates a `Context` and runs the `Imputor`s on the supplied data.

# Arguments
* `imp::Chain`: the chain to run
* `missing::Function`: the missing function to use in the `Context` to pass to the `Imputor`s
* `data::Dataset`: our data to impute
* `limit::Float64`: the missing data ration limit/threshold

# Returns
* `Dataset`: our imputed data
"""
function impute!(imp::Chain, missing::Function, data::Dataset; limit::Float64=0.1)
    ctx = Context(*(size(data)...), 0, limit, missing)

    for imputor in imp.imputors
        impute!(imputor, copy(ctx), data)
    end

    return data
end

"""
    impute!(imp::Chain, data::Dataset; limit::Float64=0.1)

Infers the missing data function from the `data` and passes that to
`impute!(imp::Chain, missing::Function, data::Dataset; limit::Float64=0.1)`.

# Arguments
* `imp::Chain`: the chain to run
* `data::Dataset`: our data to impute
* `limit::Float64`: the missing data ration limit/threshold

# Returns
* `Dataset`: our imputed data
"""
function impute!(imp::Chain, data::Dataset; limit::Float64=0.1)
    impute!(imp, ismissing, data; limit=limit)
end


text:
module FluxGogh

include(joinpath(@__DIR__, "config.jl"))
import .Config: Configuration

using Flux: Flux, Chain, Conv, relu
using LinearAlgebra
using Metalhead: Colorant, Metalhead, VGG19, channelview
using NNlib
using Images
# using Random
using Statistics: mean

loadmodel() = VGG19()

function preprocess(image::AbstractMatrix{<:Colorant})
    # Convert from CHW (Image.jl's channel ordering) to WHCN (Flux.jl's ordering)
    # and enforce Float32, as that seems important to Flux
    Float32.(permutedims(channelview(RGB.(image)), (3, 2, 1))[:, :, :, :] .* 255 .- 120)
end
preprocess(x) = x

loadimage(image) = preprocess(Metalhead.load(image))

saveimage(file, image::AbstractArray{<:Float32}) = save(file, clamp.(permutedims(image[:, :, :], (2, 1, 3)) .+ 120, 0, 255) ./ 255f0)

function choosemids(model::VGG19, x)
    layers1 = Chain(model.layers[1], Conv(weight=model.layers[2].weight, bias=model.layers[2].bias))
    x1 = layers1(x)
    layers2 = Chain(x -> relu.(x), model.layers[3:4], Conv(weight=model.layers[5].weight, bias=model.layers[5].bias))
    x2 = layers2(x1)
    layers3 = Chain(x -> relu.(x), model.layers[6:9], Conv(weight=model.layers[10].weight, bias=model.layers[10].bias))
    x3 = layers3(x2)
    layers4 = Chain(x -> relu.(x), model.layers[11:14], Conv(weight=model.layers[15].weight, bias=model.layers[15].bias))
    x4 = layers4(x3)
    [x1, x2, x3, x4]
end

function get_matrix(y)
    width, height, ch_num, batch_size = size(y)
    y_reshaped = reshape(y, (width * height, ch_num, batch_size))
    result = NNlib.batched_mul(permutedims(y_reshaped, (2,1,3)), y_reshaped)
    result ./ Float32(width * height * ch_num)
end

function generate_image(model, img_orig, img_style, config::Configuration)
    img_gen = rand(Float32, config.shape) .* 40 .- 20
    generate_image(model, img_orig, img_style, img_gen, config)
end

function generate_image(model, img_orig, img_style, img_gen, config::Configuration)
    _model = model .|> Flux.gpu
    mids_orig = choosemids(model, img_orig) .|> Flux.gpu
    mids_style = choosemids(model, img_style)  # .|> Flux.gpu
    style_mats = get_matrix.(mids_style) .|> Flux.gpu  # == [get_matrix(y) for y in mids_orig]

    function loss(x, (y1, y2)::NTuple{2, Any})
        ŷ = choosemids(_model, x)
        l = length(ŷ)
        s1 = l - l ÷ 2 + 1
        L1 = sum(Flux.Losses.mse.(ŷ[s1:end], y1[s1:end]))
        L2 = mean(Flux.Losses.mse.(get_matrix.(ŷ), y2))
        config.lam * L1 + L2
    end

    _img_gen = img_gen |> Flux.gpu
    ps = Flux.params(_img_gen)
    opt = Flux.ADAM(config.lr, (0.9f0, 0.999f0))
    for i in 1:config.iteration
        # Flux.train!(loss, ps, [(_img_gen, (mids_orig, style_mats))], opt)
        gs = Flux.gradient(ps) do
            loss(_img_gen, (mids_orig, style_mats))
        end
        Flux.update!(opt, ps, gs)
        if i % 100 == 0
            ŷ = choosemids(_model, _img_gen)
            l = length(ŷ)
            s1 = l - l ÷ 2 + 1
            L1 = Flux.Losses.mse.(ŷ[s1:end], mids_orig[s1:end])
            L2 = Flux.Losses.mse.(get_matrix.(ŷ), style_mats)
            L = config.lam * sum(L1) + mean(L2)
            @info i, L1, L2, L
            saveimage(joinpath(config.outputdir, "$(lpad(i, 5, '0')).png"), _img_gen |> Flux.cpu)
        end
    end
    _img_gen |> Flux.cpu
end

end # module


text:
@testset "Basic" begin

    d = Flux.Dense(10, 5)

    c = Flux.Chain(d, Flux.softmax)
    p = Positionwise(d, Flux.softmax)

    x = randn(10, 5, 3)

    @test p(x)[:,:,2] == c(x[:,:,2])
    @test p(x[:,:,2]) == c(x[:,:,2])
end


text:
function animate_gif(xs, y, t, x_str; x_label=["" for i in 1:length(xs)], filename=x_str, header="", directory="Output")
    filepath = pwd() * "/" * directory * "/"
    isdir(dirname(filepath)) || mkpath(filepath)

    anim = @animate for n in 1:length(t)
        x_max = maximum([maximum(x) for x in xs])
        x_min = minimum([minimum(x) for x in xs])

        fig = plot(xlim=(x_min, x_max), ylim=(minimum(y), maximum(y)), legend=:bottom, size=(400,400))
        for i in 1:length(xs)
            plot!(fig, xs[i][:,n], y, label=x_label[i], title=header*"t = $(round(t[n]/86400, digits=1)) days", linewidth=3)
        end

        xlabel!(fig, "$x_str")
        ylabel!(fig, "z")
    end

    gif(anim, pwd() * "/$(directory)/$(filename).gif", fps=20)
end


text:
mutable struct Entity
    genes:: Array{Float64,1}
    fitness:: Float64
end

function generateEntity(genesLength)
    return Entity(rand(-1:0.001:1,genesLength), 0)
end

function calculateFitness!(entity)
    rez = 4*(entity.genes[1]^2) - 6*entity.genes[1] - 3*(entity.genes[2]^3)
    + 0.5*entity.genes[2] + 3*entity.genes[3] + 8*entity.genes[4] - 6.1*entity.genes[5]
    + 2*entity.genes[6] - 10
    entity.fitness = abs(rez)
end

function printEntity(entity)
    for i in 1:length(entity.genes)
        print("$(entity.genes[i]) ")
    end
    print(" $(entity.fitness) \n");
end

function crossover!(entity1,entity2,crossoverPoint)
    for i in 1:crossoverPoint
        x = entity1.genes[i]
        entity1.genes[i] = entity2.genes[i]
        entity2.genes[i] = x
    end
end

function mutate!(entity,mutatePercentage)
    if rand(Float64) < mutatePercentage
        mutationPoint = rand(1:genesLength)
        entity.genes[mutationPoint] = 1 - entity.genes[mutationPoint]
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using ..Grids
using ..MPIStateArrays
using ..DGBalanceLawDiscretizations
using ..AtmosDycore.VanillaAtmosDiscretizations

"""
    writevtk(prefix, Q::MPIStateArray, disc::DGBalanceLaw [, fieldnames])

Write a vtk file for all the fields in the state array `Q` using geometry and
connectivity information from `disc.grid`. The filename will start with `prefix`
which may also contain a directory path. The names used for each of the fields
in the vtk file can be specified through the collection of strings `fieldnames`;
if not specified the fields names will be `"Q1"` through `"Qk"` where `k` is the
number of states in `Q`, i.e., `k = size(Q,2)`.

"""
function writevtk(prefix, Q::MPIStateArray, disc::DGBalanceLaw,
                  fieldnames=nothing)
  vgeo = disc.grid.vgeo
  host_array = Array ∈ typeof(Q).parameters
  (h_vgeo, h_Q) = host_array ? (vgeo, Q.Q) : (Array(vgeo), Array(Q))
  writevtk_helper(prefix, h_vgeo, h_Q, disc.grid, fieldnames)
end

"""
    writevtk(prefix, Q::MPIStateArray, disc::DGBalanceLaw, fieldnames,
             auxstate::MPIStateArray, auxfieldnames)

Write a vtk file for all the fields in the state array `Q` and auxiliary state
`auxstate` using geometry and connectivity information from `disc.grid`. The
filename will start with `prefix` which may also contain a directory path. The
names used for each of the fields in the vtk file can be specified through the
collection of strings `fieldnames` and `auxfieldnames`.

If `fieldnames === nothing` then the fields names will be `"Q1"` through `"Qk"`
where `k` is the number of states in `Q`, i.e., `k = size(Q,2)`.

If `auxfieldnames === nothing` then the fields names will be `"aux1"` through
`"auxk"` where `k` is the number of states in `auxstate`, i.e., `k =
size(auxstate,2)`.

"""
function writevtk(prefix, Q::MPIStateArray, disc::DGBalanceLaw,
                  fieldnames, auxstate, auxfieldnames)
  vgeo = disc.grid.vgeo
  host_array = Array ∈ typeof(Q).parameters
  (h_vgeo, h_Q, h_aux) = host_array ? (vgeo, Q.Q, auxstate.Q) :
                                      (Array(vgeo), Array(Q), Array(auxstate))
  writevtk_helper(prefix, h_vgeo, h_Q, disc.grid, fieldnames, h_aux,
                  auxfieldnames)
end

"""
    writevtk_helper(prefix, vgeo::Array, Q::Array, grid, fieldnames)

Internal helper function for `writevtk`
"""
function writevtk_helper(prefix, vgeo::Array, Q::Array, grid,
                         fieldnames, auxstate=nothing, auxfieldnames=nothing)

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq  = N+1

  nelem = size(Q)[end]
  Xid = (grid.xid, grid.yid, grid.zid)
  X = ntuple(j->reshape((@view vgeo[:, Xid[j], :]),
                        ntuple(j->Nq, dim)...,
                        nelem), dim)
  if fieldnames == nothing
    fields = ntuple(i->("Q$i", reshape((@view Q[:, i, :]),
                                       ntuple(j->Nq, dim)..., nelem)),
                    size(Q, 2))
  else
    fields = ntuple(i->(fieldnames[i], reshape((@view Q[:, i, :]),
                                               ntuple(j->Nq, dim)..., nelem)),
                    size(Q, 2))
  end
  if auxstate !== nothing
    if auxfieldnames === nothing
      auxfields = ntuple(i->("aux$i", reshape((@view auxstate[:, i, :]),
                                              ntuple(j->Nq, dim)..., nelem)),
                         size(auxstate, 2))
    else
      auxfields = ntuple(i->(auxfieldnames[i], reshape((@view auxstate[:, i, :]),
                                                       ntuple(j->Nq, dim)...,
                                                       nelem)),
                         size(auxstate, 2))
    end
    fields = (fields..., auxfields...)
  end
  writemesh(prefix, X...; fields=fields, realelems=grid.topology.realelems)
end

function writevtk(prefix, Q::MPIStateArray, disc::VanillaAtmosDiscretization)
  vgeo = disc.grid.vgeo
  host_array = Array ∈ typeof(Q).parameters
  (vgeo, Q) = host_array ? (vgeo, Q.Q) : (Array(vgeo), Array(Q))
  writevtk_VanillaAtmosDiscretization(prefix, vgeo, Q, disc.grid)
end

function writevtk_VanillaAtmosDiscretization(prefix, vgeo::Array, Q::Array,
                                             G::Grids.AbstractGrid{T, dim, N}
                                            ) where {T, dim, N}

  Nq  = N+1

  nelem = size(Q)[end]
  Xid = (G.xid, G.yid, G.zid)
  X = ntuple(j->reshape((@view vgeo[:, Xid[j], :]),
                        ntuple(j->Nq, dim)...,
                        nelem), dim)

  _ρ = VanillaAtmosDiscretizations._ρ
  _U = VanillaAtmosDiscretizations._U
  _V = VanillaAtmosDiscretizations._V
  _W = VanillaAtmosDiscretizations._W
  _E = VanillaAtmosDiscretizations._E

  ρ = reshape((@view Q[:, _ρ, :]), ntuple(j->Nq, dim)..., nelem)
  U = reshape((@view Q[:, _U, :]), ntuple(j->Nq, dim)..., nelem)
  V = reshape((@view Q[:, _V, :]), ntuple(j->Nq, dim)..., nelem)
  W = reshape((@view Q[:, _W, :]), ntuple(j->Nq, dim)..., nelem)
  E = reshape((@view Q[:, _E, :]), ntuple(j->Nq, dim)..., nelem)
  writemesh(prefix, X...;
            fields=(("ρ", ρ), ("U", U), ("V", V), ("W", W), ("E", E)),
            realelems=G.topology.realelems)
end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
# Custom Collections

"""
Base type for all custom db collections.
"""
abstract type AbstractSimStringDB end

"""
Abstract type for feature extraction structs
"""
abstract type FeatureExtractor end

# Feature Extraction Definitions

"""
Feature extraction on character-level ngrams
"""
struct CharacterNGrams{T1<:Int, T2<:AbstractString} <: FeatureExtractor
    n::T1         # number of n-grams to extract
    padder::T2    # string to use to pad n-grams
end

"""
Feature extraction based on word-level ngrams
"""
struct WordNGrams{T1<:Int, T2<:AbstractString} <: FeatureExtractor
    n::T1           # number of n-grams to extract
    padder::T2      # string to use to pad n-grams
    splitter::T2    # string to use to split words
end


text:
image = open("password.sif") do f
    reshape([parse(Int, s) for s in readlines(f)[1]], 25, 6, :)
end

# image = reshape([parse(Int, s) for s in "0222112222120000"], 2, 2, :)

function checksum(image::Array{Int,3})::Int
    min = size(image, 1)*size(image, 2)
    prod = 0
    for i in 1:size(image, 3)
        zeros = sum(x -> x == 0, image[:,:,i])
        ones = sum(x -> x == 1, image[:,:,i])
        twos = sum(x -> x == 2, image[:,:,i])
        if zeros < min
            min = zeros
            prod = ones * twos
        end
    end

    return prod
end

function visible_pixel(pixels::Array{Int,1})::Int
    for p in pixels
        # println(p)
        if p == 0 || p == 1
            return p
        end
    end
    return -1
end

function parse_image(image::Array{Int,3})::Array{Int,2}
    res = zeros(Int, size(image, 1), size(image, 2))
    for x in 1:size(image, 1), y in 1:size(image,2)
        # println(image[x, y, :])
        res[x, y] = visible_pixel(image[x, y, :])
    end
    for row in 1:size(res, 2)
        println(replace(join(res[:,row]), "0" => " "))
    end
    return res
end

println(checksum(image))
parse_image(image)


text:
# Tests for the PSet types.
# Author: Noah C. Benson <n@nben.net>

function compare_test(p::SIMM, s::SMUT, ks::Vector{T}, n::Int) where {T,SIMM<:AbstractSet{T},SMUT<:AbstractSet{T}}
    let k, q, pd = [:push, :delete], ks = collect(ks)
        for i in 1:n
            k = rand(ks)
            q = rand(pd)
            if q == :push
                p = Air.push(p, k)
                push!(s, k)
            else
                p = Air.delete(p, k)
                delete!(s, k)
            end
            @test length(p) == length(s)
            @test (k in p) == (k in s)
            @test isequal(p, s)
        end
    end
end
@testset "PSet" begin
    # First, do a standard test:
    numops = 100
    numits = 10
    for it in 1:numits
        mut = Set{Real}()
        imm = Air.PSet{Real}()
        ops = [:push, :pop, :get]
        for ii in 1:numops
            q = rand(ops)
            if q == :push
                x = rand(Float64)
                push!(mut, x)
                imm = Air.push(imm, x)
            elseif q == :pop && length(mut) > 0
                x = rand(mut)
                delete!(mut, x)
                imm = Air.delete(imm, x)
            elseif q == :get && length(mut) > 0
                x = rand(mut)
                @test in(x, imm)
            end
            @test length(mut) == length(imm)
            @test mut == imm
        end
    end
    # Next, using the above function
    let syms = [:a, :b, :c, :d, :e, :f, :g, :h, :i, :j], n = 100
        @testset "PIdSet" begin
            compare_test(Air.PIdSet{Symbol}(), Base.IdSet{Symbol}(), syms, n)
            idset = Base.IdSet{Symbol}()
            for k in [:b, :d, :e]
                push!(idset, k)
            end
            compare_test(Air.PIdSet{Symbol}([:b, :d, :e]), idset, syms, n)
        end
        @testset "PSet" begin
            compare_test(Air.PSet{Symbol}(), Set{Symbol}(), syms, n)
            compare_test(Air.PSet{Symbol}([:b, :d, :e]), Set{Symbol}([:b, :d, :e]), syms, n)
        end
    end
end


text:
function my_f(x::Float64, y::Float64)
        return 2x + 3y
end

max_stars_repo_name: Mikejmnez/MyExample

max_issues_repo_name: Mikejmnez/MyExample

max_forks_repo_name: Mikejmnez/MyExample


text:
function make_dirname(savelocation::String;
                      stub::String="unnamed",
                      tag::String="")
    println("Saving in: $savelocation")
    if isempty(tag)
        print("Enter a directory tag "*
              "(or hit Enter to generate from the date): ")
        tag = readline()
        if isempty(tag)
            tag = Libc.strftime("%d-%b-%y-at-%H-%M-%S", time())
        end
    end

    if isempty(stub)
        dirname = tag
    else
        dirname = stub*"-"*tag
    end
    pathstr = joinpath(savelocation, dirname)
    return pathstr
end

function make_dir(pathstr::String)
    if isdir(pathstr)
        @info "Directory already exists, so not copying source files."
    else
        mkpath(pathstr)
        cplist = ["LICENSE.md", "Manifest.toml", "Project.toml", "README.md",
                  "src", "tex"]
        sourcepathstr = joinpath(@__DIR__, "../")
        for c in cplist
            cp(joinpath(sourcepathstr, c),
               joinpath(pathstr, c),
               follow_symlinks=true)
        end
    end

    # Make a directory for holding output
    resultsdir = joinpath(pathstr, "results")
    already_existed = false
    if isdir(resultsdir)
        @info "Results directory already exists, not modifying."
        already_existed = true
    else
        mkdir(resultsdir)
    end

    # Copy .tex files to results directory
    texdir = joinpath(resultsdir, "../tex")
    cplist = readdir(texdir)
    for c in cplist
        cp(joinpath(texdir, c), joinpath(resultsdir, c))
    end

    return resultsdir, already_existed
end

function setup(savelocation::String;
               stub::String="unnamed",
               tag::String="")
    pathstr = make_dirname(savelocation, stub = stub, tag = tag)
    resultsdir, already_existed = make_dir(pathstr)
    return resultsdir, already_existed
end


text:
module DefaultHeadersRequest

export defaultheaderslayer, setuseragent!

using ..Messages
using ..IOExtras
import ..Form, ..content_type

"""
    defaultheaderslayer(req) -> Response

Sets default expected headers.
"""
function defaultheaderslayer(handler)
    return function(req; iofunction=nothing, kw...)
        headers = req.headers
        if isempty(req.url.port) ||
            (req.url.scheme == "http" && req.url.port == "80") ||
            (req.url.scheme == "https" && req.url.port == "443")
            hostheader = req.url.host
        else
            hostheader = req.url.host * ":" * req.url.port
        end
        defaultheader!(headers, "Host" => hostheader)
        defaultheader!(headers, "Accept" => "*/*")
        if USER_AGENT[] !== nothing
            defaultheader!(headers, "User-Agent" => USER_AGENT[])
        end

        if !hasheader(headers, "Content-Length") &&
        !hasheader(headers, "Transfer-Encoding") &&
        !hasheader(headers, "Upgrade")
            l = nbytes(req.body)
            if l !== nothing
                setheader(headers, "Content-Length" => string(l))
            elseif req.method == "GET" && iofunction isa Function
                setheader(headers, "Content-Length" => "0")
            end
        end
        if !hasheader(headers, "Content-Type") && req.body isa Form && req.method in ("POST", "PUT")
            # "Content-Type" => "multipart/form-data; boundary=..."
            setheader(headers, content_type(req.body))
        end
        return handler(req; iofunction=iofunction, kw...)
    end
end

const USER_AGENT = Ref{Union{String, Nothing}}("HTTP.jl/$VERSION")

"""
    setuseragent!(x::Union{String, Nothing})

Set the default User-Agent string to be used in each HTTP request.
Can be manually overridden by passing an explicit `User-Agent` header.
Setting `nothing` will prevent the default `User-Agent` header from being passed.
"""
function setuseragent!(x::Union{String, Nothing})
    USER_AGENT[] = x
    return
end

end # module


text:
using Distributions
using LinearAlgebra
using StatsBase

num_bins = 20

function get_hist_edges(input_dict, nbins, population)
    all_outputs = reduce(vcat, [input_dict[i] for i in keys(input_dict)])

    max_output = maximum(all_outputs)
    min_output = minimum(all_outputs)

    if (population) # percent off
        # println("Population CC bin edges!")
        bin_edges = collect(-1:102.0/nbins:101.0)
    else # FITC
        # println("Single cell CC bin edges!")
        bin_edges = collect(min_output-1:(max_output+2-min_output)/nbins:max_output+1)
    end

    h = fit(Histogram, all_outputs, bin_edges)
    return collect(h.edges[1])
end

function find_bin_index(value, bin_edges) # wheeee no bounds checking here
    if debug
        print("Bin edges: ")
        println(bin_edges)
        print("Value: ")
        println(value)
    end
    for i in 1:length(bin_edges)-1
        if value >= bin_edges[i] && value < bin_edges[i+1]
            if debug
                print("Found index at: ")
                println(i)
            end
            return i
        else
            if debug
                print("Did not find index at ")
                println(i)
            end
        end
    end
end

# each row of matrix is input, each column is output
function build_conditional_matrix(input_dict, population)
    inputs = collect(keys(input_dict))
    output_bin_edges = get_hist_edges(input_dict, num_bins, population) # lower edges, bin into 20 things
    freqs = zeros(length(inputs), length(output_bin_edges)-1)

    if debug
        println("Frequencies: ")
        println(freqs)
    end

    for i in 1:length(inputs)
        input = inputs[i]
        for j in 1:length(input_dict[input])
            output = input_dict[input][j]
            bin_idx = find_bin_index(output, output_bin_edges)

            if isnothing(bin_idx)
                println(output_bin_edges)
                println(output)
            end
            if debug
                print("Adding count to ")
                print(i)
                print(", ")
                println(bin_idx)
            end
            freqs[i, bin_idx] += 1
        end
    end

    if debug
        println("Frequencies: ")
        println(freqs)
    end

    # now we have to normalize each row
    for i in 1:size(freqs)[1]
        tot = sum(freqs[i,:])
        for j in 1:length(output_bin_edges)-1
            freqs[i, j] = freqs[i, j]/tot
        end
    end

    if debug
        println("Frequencies: ")
        println(freqs)
    end

    return freqs
end

function blahut_arimoto(x)
    numrows, numcols = size(x)

    # check for negative entries
    for i in 1:numrows
        for j in 1:numcols
            if (x[i,j]<0)
                println("ERROR: negative entry in input matrix")
                return
            end
        end
    end

    # check for zero columns
    for i in 1:numcols
        arr = x[:,i]
        if arr == zeros(1, numrows)
            println("ERROR: there is a zero column in the matrix")
            return;
        end
    end

    # check for zero rows and make sure they sum to 1
    for i in 1:numrows
        arr = x[i,:]
        if arr == zeros(1, numcols)
            println("ERROR: there is a zero row in the matrix")
            return;
        end
        x[i,:] = x[i,:]/sum(x[i,:]) # normalize rows
    end

    r = ones(numrows)/numrows
    q = zeros(numrows, numcols)
    error_tolerance = 1e-5/numrows
    r1 = zeros(numrows)

    for i in 1:numrows
        x[i,:] = x[i,:]/sum(x[i,:])
    end
    for iter in 1:10000
        for j in 1:numcols
            q[:,j] = r .* x[:,j]
            q[:,j] = q[:,j]/sum(q[:,j])
        end
        for i in 1:numrows
            r1[i] = prod(q[i,:] .^ x[i,:])
        end
        r1 = r1/sum(r1)
        r = r1

        #  if (norm(r1-r) < error_tolerance)
        #      break
        #  else
        #      r = r1
        #  end
    end

    C = 0
    for i in 1:numrows
        for j in 1:numcols
            if r[i] > 0 && q[i,j] > 0
                C = C + r[i]*x[i,j]*log(q[i,j]/r[i])
            end
        end
    end

    C = C / log(2)
    return (C, r)
end

#= Computes the channel capacity given a dict of input => outputs
=#
function compute_channel_capacity(input_dict, population=true)
    m = build_conditional_matrix(input_dict, population)
    C, r = blahut_arimoto(m)
    return C
end


text:
module TestExamples1DMHD

using Test
using Trixi

include("test_trixi.jl")

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

@testset "MHD" begin

  @testset "elixir_mhd_alfven_wave.jl with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.440611823425164e-15, 1.1373567770134494e-14, 3.024482376149653e-15, 2.0553143516814395e-15, 3.9938347410210535e-14, 3.984545392098788e-16, 2.4782402104201577e-15, 1.551737464879987e-15],
      linf = [1.9984014443252818e-15, 1.3405943022348765e-14, 3.3584246494910985e-15, 3.164135620181696e-15, 7.815970093361102e-14, 8.881784197001252e-16, 2.886579864025407e-15, 2.942091015256665e-15],
      initial_condition = initial_condition_constant,
      tspan = (0.0,1.0))
  end

  @testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.4717982305884098e-5, 1.182263100132273e-5, 2.3571203951257584e-5, 2.3571203951305505e-5, 1.8950292712340346e-6, 1.196322416573199e-16, 2.3661336207200768e-5, 2.366133620713721e-5],
      linf = [6.156423568670633e-5, 3.474803933150424e-5, 0.00011882557236665703, 0.0001188255723667403, 6.861887285491974e-6, 2.220446049250313e-16, 0.00012195115526943134, 0.00012195115526922318])
  end

  @testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.15286156171121895, 0.16926639610967018, 0.30495838128766584, 0.0, 0.345735219586204, 6.608258910237594e-17, 0.31766258436717354, 0.0],
      linf = [0.7012110703585162, 0.508655909370987, 1.251766990108334, 0.0, 1.9957985139014158, 1.1102230246251565e-16, 1.6877793846395974, 0.0])
  end

  @testset "elixir_mhd_briowu_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_briowu_shock_tube.jl"),
      l2   = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0],
      linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0])
  end

  @testset "elixir_mhd_torrilhon_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
      l2   = [0.45702524713047815, 0.4792453322653316, 0.3406561762729464, 0.4478037841848423, 0.9204646319222438, 1.3216517820475193e-16, 0.2889772516628848, 0.2552071770333965],
      linf = [1.221507578527215, 0.8927062606782392, 0.8536666880170606, 0.9739151339925485, 1.6602343391033632, 2.220446049250313e-16, 0.686743631371773, 0.6555428379163893])
  end

  @testset "elixir_mhd_ryujones_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"),
      l2   = [0.23480132779180388, 0.3922573469657149, 0.08236145370312341, 0.17559850888696307, 0.9615951543069574, 6.608258910237605e-17, 0.21534816808200136, 0.10699642480726586],
      linf = [0.6397374584516875, 0.9457640306380323, 0.3546254702548819, 0.8532865350683194, 2.078210620993694, 1.1102230246251565e-16, 0.49334335395420337, 0.2504525070653655],
      tspan = (0.0, 0.1))
  end

  @testset "elixir_mhd_shu_osher_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"),
      l2   = [1.011173124884092, 8.271206982561576, 1.3087914055352385, 0.0, 52.19213714506333, 6.743357635607637e-16, 1.0102652381088522, 0.0],
      linf = [2.8649230108993886, 22.600699637929466, 4.166209106166631, 0.0, 135.11479022821882, 3.4416913763379853e-15, 2.83396622829752, 0.0],
      tspan = (0.0, 0.2))
  end

  @testset "elixir_mhd_shu_osher_shock_tube.jl with flipped shock direction" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"),
      l2   = [1.0150653343383023, 8.293796006442639, 1.2953342215936525, 0.0, 52.34290101591115, 3.352946970627566e-16, 1.0046607148062305, 0.0],
      linf = [2.9091039190988806, 22.734543670960285, 4.124890055211936, 0.0, 136.89057030285923, 1.5543122344752192e-15, 2.806658251289774, 0.0],
      initial_condition = initial_condition_shu_osher_shock_tube_flipped,
      boundary_conditions=boundary_condition_shu_osher_shock_tube_flipped,
      tspan = (0.0, 0.2))
  end
end

end # module


text:
module ScalingWithExtrapTests

using Interpolations, Base.Test

xs = linspace(-5, 5, 10)
ys = sin(xs)

function run_tests{T,N,IT}(sut::Interpolations.AbstractInterpolation{T,N,IT,OnGrid}, itp)
    for x in xs
        @test_approx_eq_eps sut[x] sin(x) sqrt(eps(sin(x)))
    end
    @test sut[-5] == sut[-5.1] == sut[-15.8] == sut[-Inf] == itp[1]
    @test sut[5] == sut[5.1] == sut[15.8] == sut[Inf] == itp[end]
end

function run_tests{T,N,IT}(sut::Interpolations.AbstractInterpolation{T,N,IT,OnCell}, itp)
    halfcell = (xs[2] - xs[1]) / 2

    for x in (5 + halfcell, 5 + 1.1halfcell, 15.8, Inf)
        @test sut[-x] == itp[.5]
        @test sut[x] == itp[end+.5]
    end
end

for GT in (OnGrid, OnCell)
    itp = interpolate(ys, BSpline(Quadratic(Flat())), GT())

    # Test extrapolating, then scaling
    eitp = extrapolate(itp, Flat())
    seitp = scale(eitp, xs)
    run_tests(seitp, itp)

    # Test scaling, then extrapolating
    sitp = scale(itp, xs)
    esitp = extrapolate(sitp, Flat())
    run_tests(esitp, itp)
end

end


text:
InputFolder = './Images/NucleiNoisy/';
OutputFolder = './Results/Images/NucleiNoisy2/';

@iA = '*.tif';

@fxg_gDenoiseDCT [iA] > [D, An];
params.WinSize = 32;
params.Step = 2;
params.Thresh = 0.45;
params.TopOpen = 0;
params.AddNoiseVar = 0;
/endf

/show iA >;
/show D >;
/keep D > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
## Number of observations

nobs(x::AbstractVector) = length(x)

nobs(x::AbstractArray, dim::Integer = 1) = size(x, dim)

nobs(x::DataFrames.AbstractDataFrame) = nrow(x)


text:
using NumericalMethodsforEngineers, Test

#=
@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1
          den = 1
          For( j=1, j <= N, Increment(j),
            begin
              If(j != i, num *= (x-xi[j]))
              If(j != i, den *= (xi[i]-xi[j]))
            end
          )
          sum += yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
  xi = [0, 1//2, 1]
  yi = [0, a, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  Println(Y(x))
end
=#

@sym begin
  ClearAll(xi, yi, Y, R, sol, ytilde)
  xi = [0, 1//2, 1]
  yi = [0, a, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  #
  # Can be formulated as ytile(x) = F(x) + C1(a) * Ψ(x)
  #
  # F(x_) := 2*x^2 - x
  # C(a_) := -4a
  # Ψ(x_) := x^2 - x
  #
  #Y(x_) := F(x) + C(a)*Ψ(x)
  #
  R(x_) = Simplify(D(Y(x), x, 2) - 3*x - 4*Y(x))
  R(x_) = Simplify(R(x) ./ (a => -C1/4))
  SetJ(r, ToString(Simplify(R(x))))
end

println("\n\nExample 7.13: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("Residual for Weighted Residual Method using 1 point Lagragian Polynomial")
@sym Println("\nY(x) = ", Y(x), "\n")
@sym Println("R(x) = ", R(x), "\n")
println("( Example 7.13 gives: R = -4x^2*(2 + C1) + x*(1 + 4C1) + 2*(2 + C1) )", "\n")

println(r)

@eval f713(x, C1) = $(Meta.parse(r))

@test f713(1, 2) == 4.0 + 2.0*2 - 3 - 4*(-1.0 + 2.0*1 + 2*(-1 + 1))


text:
function Segment(arg0::Vector3D, arg1::Vector3D, arg2::Line)
    return Segment((Vector3D, Vector3D, Line), arg0, arg1, arg2)
end

function get_end(obj::Segment)
    return jcall(obj, "getEnd", Vector3D, ())
end

function get_line(obj::Segment)
    return jcall(obj, "getLine", Line, ())
end

function get_start(obj::Segment)
    return jcall(obj, "getStart", Vector3D, ())
end


text:
using DynamicalSystems
using Base.Test

# Test 1: harmonicOscillator.

omega0_test1 = 2.0;
beta_test1 = 0.5;
gamma_test1 = 1.5;
omega_test1 = 1.0;

t_test1 = 3.0;
x_test1 = [1.0, 2.0];

functions_test1 = harmonicOscillator(omega0_test1, beta=beta_test1,
    gamma=gamma_test1, omega=omega_test1);

# Validate the results of the functions.
truth_test1 = [2.0, -4.788319987910199];
answer_test1 = map(fe -> fe(t_test1, x_test1), functions_test1);

@test_approx_eq answer_test1 truth_test1

# Test 2: duffingOscillator.

delta_test2 = 0.5;
gamma_test2 = 1.5;
alpha_test2 = 0.25;
beta_test2 = -0.75;
omega_test2 = 2.5;
t_test2 = 3.0;
x_test2 = [1.5, 3.5];

functions_test2 = duffingOscillator(delta_test2, gamma_test2,
    alpha=alpha_test2, beta=beta_test2, omega=omega_test2);

# Validate the results of the functions.
truth_test2 = [3.5, 1.8132499651621083];
answer_test2 = map(fe -> fe(t_test2, x_test2), functions_test2);

@test_approx_eq answer_test2 truth_test2

# Test 3: lorenzAttractor.

sigma_test3 = 0.25;
rho_test3 = 0.5;
beta_test3 = 0.75;
t_test3 = 0.0;
x_test3 = [0.55; 0.35; 0.15];

functions_test3 = lorenzAttractor(sigma_test3, rho_test3, beta_test3);

# Validate the results of the functions.
truth_test3 = [-0.05, -0.15749999999999, 0.08];
answer_test3 = map(fe -> fe(t_test3, x_test3), functions_test3);

@test_approx_eq answer_test3 truth_test3


text:
using HDF5

Base.@kwdef mutable struct TelInfo <: Info
    antenna_positions::AbstractArray{Float64,2} = Array{Float64,2}(undef, 3, 0)
    antenna_position_frame::Union{AbstractString,Nothing} = nothing
    antenna_names::AbstractVector{AbstractString} = Vector{String}(undef, 0)
    antenna_numbers::AbstractVector{Int} = Vector{Int}(undef, 0)
    antenna_diameters::AbstractVector{Float64} = Vector{Float64}(undef, 0)
    latitude::Float64 = NaN
    longitude::Float64 = NaN
    altitude::Float64 = NaN
    telescope_name::AbstractString = "Unknown"
end # mutable struct TelInfo

function to_hdf5(group::HDF5.Group, info::TelInfo)
    invoke(to_hdf5, Tuple{HDF5.Group, Info}, group, info;
           antenna_frame=(layout=HDF5.H5D_COMPACT,),
           latitude=(layout=HDF5.H5D_COMPACT,),
           longitude=(layout=HDF5.H5D_COMPACT,),
           altitude=(layout=HDF5.H5D_COMPACT,),
           telescope_name=(layout=HDF5.H5D_COMPACT,)
          )
end