text:
module Arsenic

    using Gallium
    using Cxx
    using DataStructures
    using DWARF
    using DWARF: CallFrameInfo
    using ObjFileBase
    using ObjFileBase: handle
    using ELF
    using COFF

    function get_insts(session, modules, ip)
        base, mod = Gallium.find_module(session, modules, UInt(ip))
        modrel = UInt(UInt(ip)-base)
        if isa(mod, Gallium.SyntheticModule)
            bounds = mod.get_proc_bounds(session, ip)
            insts = Gallium.load(session, Gallium.RemotePtr{UInt8}(base+first(bounds)),
              length(bounds))
            return base, first(bounds), insts
        end
        if isnull(mod.xpdata)
            loc, fde = Gallium.Unwinder.find_fde(mod, modrel)
            seekloc = loc
            cie = CallFrameInfo.realize_cie(fde)
            nbytes = UInt(CallFrameInfo.fde_range(fde, cie))
        else
            entry = Gallium.Unwinder.find_seh_entry(mod, modrel)
            loc = entry.start
            # Need to translate from virtual to file addresses. Hardcode 0xa00 for
            # now.
            seekloc = loc - 0xa00
            nbytes = entry.stop - entry.start
        end
        if ObjFileBase.isrelocatable(handle(mod))
            # For JIT frames, base is the start of .text, so we need to add that
            # offset back
            text = first(filter(x->sectionname(x)==
                ObjFileBase.mangle_sname(handle(mod),"text"),Sections(handle(mod))))
            seekloc += sectionoffset(text)
        elseif isa(handle(mod), ELF.ELFHandle)
            # In weird executables, the executable segment may not be at offset 0
            phs = ELF.ProgramHeaders(handle(mod))
            idx = findfirst(p->p.p_type==ELF.PT_LOAD &&
                               ((p.p_flags & ELF.PF_X) != 0), phs)
            seekloc += phs[idx].p_offset
        elseif isa(handle(mod), COFF.COFFHandle)
            text = ObjFileBase.deref(first(filter(x->sectionname(x)==
                ObjFileBase.mangle_sname(handle(mod),"text"),Sections(handle(mod)))))
            seekloc -= text.VirtualAddress
            seekloc += text.PointerToRawData
        end
        seek(handle(mod), seekloc)
        insts = read(handle(mod), UInt8, nbytes)
        base, loc, insts
    end

    function compute_stack(modules, session::Gallium.Ptrace.Session)
        regs = Gallium.getregs(session)
        stack, RCs = Gallium.stackwalk(regs, session, modules, rich_c = true, collectRCs = true)
        Gallium.NativeStack(stack,RCs,modules,session)
    end

    function compute_stack(modules, session::Gallium.FakeMemorySession, RC)
        stack, RCs = Gallium.stackwalk(RC, session, modules, rich_c = true, collectRCs = true)
        Gallium.NativeStack(stack,RCs,modules,session)
    end

    # We should be in the same frame, fast path stack unwinding
    function update_stack_same_frame!(state, session = state.top_interp.session)
        regs = Gallium.getregs(session)
        modules = state.top_interp.modules
        oldstack, oldRCs = state.top_interp.stack, state.top_interp.RCs
        (length(oldRCs) == 1) &&
            return update_stack!(state, session)
        (ok, oneupRC) = try
            Gallium.Unwinder.unwind_step(session, modules, regs; stacktop = true, ip_only = false)
        catch err
            @show err
            return update_stack!(state, session)
        end
        parentsmatch = Gallium.ip(oldRCs[end-1]) == Gallium.ip(oneupRC) &&
            Gallium.get_dwarf(oldRCs[end-1], :rsp) == Gallium.get_dwarf(oneupRC, :rsp)
        (!ok || !parentsmatch) && return update_stack!(state, session)
        stack, RCs = copy(oldstack), copy(oldRCs)
        ok, frame = Gallium.frameinfo(regs, session, modules; rich_c = true)
        !ok && return update_stack!(state, session)
        stack[end] = get(frame)
        RCs[end] = regs
        state.interp = state.top_interp = Gallium.NativeStack(stack,RCs,modules,session)
    end

    function update_stack!(state, session = state.top_interp.session)
        state.interp = state.top_interp = compute_stack(state.top_interp.modules, session)
    end
    update_stack!(state::Void, _ = nothing) = nothing

    include("disassembler.jl")
    include("stepping.jl")
    include("remoteexec.jl")
    is_linux() && isfile(Pkg.dir("RR","src","RR.jl")) && include("arrsenic.jl")
    include("interface.jl")
end # module


text:
# Source: https://www.lsac.org/lsat/taking-lsat/test-format/analytical-reasoning/analytical-reasoning-sample-questions
#
# A university library budget committee must reduce exactly five of eight areas of
# expenditure—G, L, M, N, P, R, S, and W—in accordance with the following conditions:
#
# 1 If both G and S are reduced, W is also reduced.
# 2 If N is reduced, neither R nor S is reduced.
# 3 If P is reduced, L is not reduced.
# 4 Of the three areas L, M, and R, exactly two are reduced.

julia> using AbstractLogic
Start the repl in command prompt by typing `=`.

abstractlogic> G, L, M, N, P, R, S, W ∈ red, not [clear]

# 1 If both G and S are reduced, W is also reduced.
abstractlogic> G, S = 'red' ==> W = 'red'

# 2 If N is reduced, neither R nor S is reduced.
abstractlogic> N = 'red' ==> R, S = 'not'

# 3 If P is reduced, L is not reduced.
abstractlogic> P = 'red' ==> L = 'not'

# 4 Of the three areas L, M, and R, exactly two are reduced.
abstractlogic> L|M|R = 'red' {2}

# Question 1
# If both M and R are reduced,
# Save state
abstractlogic> M, R = 'red'

# which one of the following is a pair of areas neither of which could be reduced?
# A. G, L
abstractlogic> all: G, L = 'not'

# B. G, N
abstractlogic> all: G, N = 'not'

# C. L, N
abstractlogic> all: L, N = 'not'
# D. L, P
abstractlogic> all: L, P = 'not'

# E. P, S
abstractlogic> all: P, S = 'not'

####
# Passage for Questions 2 and 3
# Seven piano students—T, U, V, W, X, Y, and Z—are to give a recital, and their
# instructor is deciding the order in which they will perform. Each student will
# play exactly one piece, a piano solo. In deciding the order of performance, the
# instructor must observe the following restrictions:

# 1. X cannot play first or second.
# 2. W cannot play until X has played.
# 3. Neither T nor Y can play seventh.
# 4. Either Y or Z must play immediately after W plays.
# 5. V must play either immediately after or immediately before U plays.

# Each students plays exactly one solo
abstractlogic> T, U, V, W, X, Y, Z ∈ unique [clear]

# 1. X cannot play first or second.
abstractlogic> X != 1,2

# 2. W cannot play until X has played.
abstractlogic> W > X

# 3. Neither T nor Y can play seventh.
abstractlogic> T,Y != 7

# 4. Either Y or Z must play immediately after W plays.
abstractlogic> Y|Z = W + 1

# 5. V must play either immediately after or immediately before U plays.
abstractlogic> V = U + 1 ||| V = U - 1

# Question 2
# If V plays first, which one of the following must be true?

# First lets save our state since Question 1 creates an addition constraint
# which we would like to reverse.
abstractlogic> export preserver

abstractlogic> V = 1

# A. T plays sixth.
abstractlogic> all: T = 6

# B. X plays third.
abstractlogic> all: X = 3

# C. Z plays seventh.
abstractlogic> all: Z = 7

# D. T plays immediately after Y.
abstractlogic> all: T = Y + 1

# E. W plays immediately after X.
abstractlogic> all: W = X + 1

# lets restore the previous state (we could also use the back command but that is more likely to lead to errors)
abstractlogic> import preserver

# Question 3
# If U plays third, what is the latest position in which Y can play?
abstractlogic> U = 3

# The fast answer is:
abstractlogic> range Y

# first
# second
# fifth
# sixth - true
# seventh

# Passage for Question 4
# A charitable foundation awards grants in exactly four areas—medical services,
# theater arts, wildlife preservation, and youth services—each grant being in
# one of these areas.
abstractlogic> q1.m, q1.t, q1.w, q1.y ∈ 0:1 [clear]
abstractlogic> q2.m, q2.t, q2.w, q2.y ∈ 0:1
abstractlogic> q3.m, q3.t, q3.w, q3.y ∈ 0:1
abstractlogic> q4.m, q4.t, q4.w, q4.y ∈ 0:1
# One or more grants are awarded in each of the four quarters of a calendar year.

abstractlogic> {{>=:q1.m,<=:q1.y}} == 1 {{1,}}
abstractlogic> {{>=:q2.m,<=:q2.y}} == 1 {{1,}}
abstractlogic> {{>=:q3.m,<=:q3.y}} == 1 {{1,}}
abstractlogic> {{>=:q4.m,<=:q4.y}} == 1 {{1,}}

# Additionally, over the course of a calendar year, the following must obtain:

# Grants are awarded in all four areas.
abstractlogic> {{j}}.m = 1 {{1,}}; {{j}}.t = 1 {{1,}}; {{j}}.w = 1 {{1,}}; {{j}}.y = 1 {{1,}}

# No more than six grants are awarded.
abstractlogic> {{i}} = 1 {{,6}}

# No grants in the same area are awarded in the same quarter or in consecutive quarters.
# Medical and theater arts
abstractlogic> {{j}}.m = 1 ==> {{j+1}}.m != 1; {{j}}.t = 1 ==> {{j+1}}.t != 1

# wildlife and youthservice
abstractlogic> {{j}}.w = 1 ==> {{j+1}}.w != 1; {{j}}.y = 1 ==> {{j+1}}.y != 1

# Exactly two medical services grants are awarded.
abstractlogic> {{j}}.m = 1 {{2}}

# A wildlife preservation grant is awarded in the second quarter.
abstractlogic> q2.w = 1

# Question 4
# If a wildlife preservation grant and a youth services grant are awarded in the
# same quarter of a particular calendar year, then any of the following could be
# true that year EXCEPT:
abstractlogic> export preserver

abstractlogic> {{j}}.w = 1 &&& {{j}}.y = 1 {{1,}}

# A. A medical services grant is awarded in the second quarter.
abstractlogic> any: q2.m = 1

# B. A theater arts grant is awarded in the first quarter.
abstractlogic> any: q1.t = 1

# C. A theater arts grant is awarded in the second quarter.
abstractlogic> any: q2.t = 1

# D. A wildlife preservation grant is awarded in the fourth quarter.
abstractlogic> any: q4.w = 1

# E. A youth services grant is awarded in the third quarter.
abstractlogic> any: q3.y = 1

# Passage for Questions 5 and 6
# From a group of seven people—J, K, L, M, N, P, and Q—exactly four will be selected
# to attend a diplomat’s retirement dinner. Selection conforms to the following conditions:
abstractlogic> J, K, L, M, N, P, Q ∈ 0:1 [clear]

abstractlogic> {{i}} = 1 {{4}}

# Either J or K must be selected, but J and K cannot both be selected.
abstractlogic> J ^ K

# Either N or P must be selected, but N and P cannot both be selected.
abstractlogic> N ^ P

# N cannot be selected unless L is selected.
abstractlogic> L! ==> N!

# Q cannot be selected unless K is selected.
abstractlogic> K! ==> !Q

# Question 5
abstractlogic> export preserver
# If P is not selected to attend the retirement dinner, then exactly how many
# different groups of four are there each of which would be an acceptable selection?
abstractlogic> P = 1
abstractlogic> show
# one
# two
# three
# four - we can see from the output that 4 is available
# five
abstractlogic> import preserver

# Question 6
# There is only one acceptable group of four that can be selected to attend the
# retirement dinner if which one of the following pairs of people is selected?

# J and L
abstractlogic> check: J == 1 &&& L == 1

# K and M
abstractlogic> check: K == 1 &&& M == 1

# L and N
abstractlogic> check: L == 1 &&& N == 1

# L and Q
abstractlogic> check: L == 1 &&& Q == 1

# M and Q
abstractlogic> check: M == 1 &&& Q == 1

# Explanation for Question 6
#
# Passage for Questions 7 and 8
# On a particular Saturday, a student will perform six activities—grocery
# shopping, hedge trimming, jogging, kitchen cleaning, laundry, and motorbike servicing.
# Each activity will be performed once, one at a time.
abstractlogic> s, h, j, k, l, m ∈ unique [clear]
# The order in which the activities are performed is subject to the following conditions:

# Grocery shopping has to be immediately after hedge trimming.
abstractlogic> s = h + 1

# Kitchen cleaning has to be earlier than grocery shopping.
abstractlogic> k < s

# Motorbike servicing has to be earlier than laundry.
abstractlogic> m < l

# Motorbike servicing has to be either immediately before or immediately after jogging.
abstractlogic> m = j + 1 ||| m = j - 1

# Question 7
# If laundry is earlier than kitchen cleaning, then hedge trimming must be

abstractlogic> l < k

abstractlogic> range h
# fifth
# fourth
# third
# second
# first

# Question 8
# Which one of the following, if substituted for the condition that motorbike
# servicing has to be earlier than laundry,

# remove `m < l`
abstractlogic> s, h, j, k, l, m ∈ unique [clear]

abstractlogic> s = h + 1; k < s; m = j + 1 ||| m = j - 1

# would have the same effect in
# determining the order of the student’s activities?

abstractlogic> export preserver

abstractlogic> m < l

abstractlogic> export refset

abstractlogic> import preserver

# Laundry has to be one of the last three activities.
abstractlogic> l >=4

abstractlogic> compare refset

# Laundry has to be either immediately before or immediately after jogging.
abstractlogic> import preserver

abstractlogic> l = j + 1 ||| l = j - 1

abstractlogic> compare refset

abstractlogic> import preserver

# Jogging has to be earlier than laundry.
abstractlogic> j < l

abstractlogic> compare refset

abstractlogic> import preserver

# Laundry has to be earlier than hedge trimming.
abstractlogic> l < h

abstractlogic> compare refset

abstractlogic> import preserver

# Laundry has to be earlier than jogging.
abstractlogic> l < j

abstractlogic> compare refset

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Knapsack problem/Bounded

using MathProgBase

immutable KPDSupply{S<:String, T<:Integer}
    item::S
    weight::T
    value::T
    quant::T
end

function solve{S<:String, T<:Integer}(gear::Array{KPDSupply{S,T},1},
                                      capacity::T)
    w = map(x->x.weight, gear)
    v = map(x->x.value, gear)
    q = map(x->x.quant, gear)
    sol = mixintprog(-v, w', '<', capacity, :Int, 0, q)
    sol.status == :Optimal || error("This Problem could not be solved")
    if all(q .== 1)
        return gear[sol.sol .== 1.0]
    else
        pack = KPDSupply[]
        s = int(sol.sol)
        for (i, g) in enumerate(gear)
            s[i] != 0 || continue
            push!(pack, KPDSupply(g.item, g.weight, g.value, s[i]))
        end
        return pack
    end
end

gear = [KPDSupply("map", 9, 150, 1),
        KPDSupply("compass", 13, 35, 1),
        KPDSupply("water", 153, 200, 2),
        KPDSupply("sandwich", 50, 60, 2),
        KPDSupply("glucose", 15, 60, 2),
        KPDSupply("tin", 68, 45, 3),
        KPDSupply("banana", 27, 60, 3),
        KPDSupply("apple", 39, 40, 3),
        KPDSupply("cheese", 23, 30, 1),
        KPDSupply("beer", 52, 10, 3),
        KPDSupply("suntan cream", 11, 70, 1),
        KPDSupply("camera", 32, 30, 1),
        KPDSupply("T-shirt", 24, 15, 2),
        KPDSupply("trousers", 48, 10, 2),
        KPDSupply("umbrella", 73, 40, 1),
        KPDSupply("waterproof trousers", 42, 70, 1),
        KPDSupply("waterproof overclothes", 43, 75, 1),
        KPDSupply("note-case", 22, 80, 1),
        KPDSupply("sunglasses", 7, 20, 1),
        KPDSupply("towel", 18, 12, 2),
        KPDSupply("socks", 4, 50, 1),
        KPDSupply("book", 30, 10, 2)]

pack = solve(gear, 400)

println("The hiker should pack:")
for s in pack
    println("    ", s.quant, " ", s.item)
end
println()
println("Packed Weight:  ", mapreduce(x->x.weight*x.quant, +, pack))
println("Packed Value:  ", mapreduce(x->x.value*x.quant, +, pack))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
const global_application = Ref{Application}()

struct DisplayInline
    dom_function
end

"""
    with_session(f)::DisplayInline

calls f with the session, that will become active when displaying the result
of with_session. f is expected to return a valid DOM.
"""
function with_session(f)
    return DisplayInline(f)
end

const WebMimes = (
    MIME"text/html",
    MIME"application/prs.juno.plotpane+html",
    # MIME"application/vnd.webio.application+html"
)

function get_global_app()
    if !isassigned(global_application) || istaskdone(global_application[].server_task[])
        global_application[] = Application(
            (ctx, request)-> "Nothing to see",
            get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1"),
            parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")),
            verbose = get(ENV, "JSCALL_VERBOSITY_LEVEL", "false") == "true"
        )
    end
    global_application[]
end

for M in WebMimes
    @eval function Base.show(io::IO, m::$M, dom::DisplayInline)
        application = get_global_app()
        session_url = "/show"
        route!(application, session_url) do context
            # Serve the actual content
            return serve_dom(context, dom.dom_function)
        end
        # Display the route we just added in an iframe inline:
        println(io, "<iframe src=$(repr(url(session_url))) frameborder=\"0\" width = '100%' height = '100%'>")
        println(io, "</iframe>")
    end
end
# function Base.show(io::IO, m::MIME"text/html", dom::Hyperscript.Node)
#     inline_display = with_session() do session
#         dom
#     end
#     show(io, m, inline_display)
# end
# function Base.show(io::IO, m::MIME"text/html", dom::Markdown.MD)
#     inline_display = with_session() do session
#         dom
#     end
#     show(io, m, inline_display)
# end
function Base.show(io::IO, m::MIME"application/vnd.webio.application+html", dom::DisplayInline)
    application = get_global_app()
    session = Session()
    sessionid = string(uuid4())
    application.sessions[sessionid] = Dict("base" => session)
    dom2html(io, session, sessionid, dom.dom_function(session, (target = "/show",)))
end


text:
function name2uuid(name::String)
    u4 = UUIDs.UUID("07f2c1e0-90b0-56cf-bda7-b44b56e34eed")  # Modia3D package UUID
    uuid = string(UUIDs.uuid5(u4, name))
end

function exportObject(object, shapes, obj::Modia3D.Composition.Object3D, sphere::Modia3D.Shapes.Sphere, initPos, initRot)
    r_obj = Modia3D.ZeroVector3D
    R_obj = Modia3D.NullRotation
    name = String(Modia3D.fullName(obj)) * ".geometry"
    shape = (; name=name, uuid=name2uuid(name), type="SphereBufferGeometry", radius=sphere.diameter/2, heightSegments=16, widthSegments=32)
    material = printVisuMaterialToJSON(obj, obj.visualMaterial)
    objectInfo = getObjectInfo(obj, shape, material, initPos, initRot)
    printInfoToFile(object, shapes, shape, material, objectInfo)
    return (r_obj, R_obj)
end

function exportObject(object, shapes, obj::Modia3D.Composition.Object3D, ellipsoid::Modia3D.Shapes.Ellipsoid, initPos, initRot)
    r_obj = Modia3D.ZeroVector3D
    R_obj = Modia3D.NullRotation
    name = String(Modia3D.fullName(obj)) * ".geometry"
    shape = (; name=name, uuid=name2uuid(name), type="SphereBufferGeometry", radius=ellipsoid.lengthX/2, heightSegments=16, widthSegments=32)
    material = printVisuMaterialToJSON(obj, obj.visualMaterial)
    objectInfo = getObjectInfo(obj, shape, material, initPos, initRot, scale=[1.0, ellipsoid.lengthY/ellipsoid.lengthX, ellipsoid.lengthZ/ellipsoid.lengthX])
    printInfoToFile(object, shapes, shape, material, objectInfo)
    return (r_obj, R_obj)
end

function exportObject(object, shapes, obj::Modia3D.Composition.Object3D, box::Modia3D.Shapes.Box, initPos, initRot)
    r_obj = Modia3D.ZeroVector3D
    R_obj = Modia3D.NullRotation
    name = String(Modia3D.fullName(obj)) * ".geometry"
    shape = (; name=name, uuid=name2uuid(name), type="BoxBufferGeometry", width=box.lengthX, height=box.lengthY, depth=box.lengthZ)
    material = printVisuMaterialToJSON(obj, obj.visualMaterial)
    objectInfo = getObjectInfo(obj, shape, material, initPos, initRot)
    printInfoToFile(object, shapes, shape, material, objectInfo)
    return (r_obj, R_obj)
end

function exportObject(object, shapes, obj::Modia3D.Composition.Object3D, cylinder::Modia3D.Shapes.Cylinder, initPos, initRot)
    r_obj = Modia3D.ZeroVector3D
    if cylinder.axis == 1
        R_obj = @SMatrix[0 -1 0; 1 0 0; 0 0 1]
    elseif cylinder.axis == 2
        R_obj = Modia3D.NullRotation
    else
        R_obj = @SMatrix[1 0 0; 0 0 1; 0 -1 0]
    end
    name = String(Modia3D.fullName(obj)) * ".geometry"
    shape = (; name=name, uuid=name2uuid(name), type="CylinderBufferGeometry", radiusBottom=cylinder.diameter/2, radiusTop=cylinder.diameter/2, height=cylinder.length, radialSegments=32, heightSegments=1)
    material = printVisuMaterialToJSON(obj, obj.visualMaterial)
    objectInfo = getObjectInfo(obj, shape, material, initPos, initRot, R_obj=R_obj)
    printInfoToFile(object, shapes, shape, material, objectInfo)
    return (r_obj, R_obj)
end

function exportObject(object, shapes, obj::Modia3D.Composition.Object3D, cone::Modia3D.Shapes.Cone, initPos, initRot)
    if cone.axis == 1
        r_obj = @SVector[cone.length/2, 0.0, 0.0]
        R_obj = @SMatrix[0 -1 0; 1 0 0; 0 0 1]
    elseif cone.axis == 2
        r_obj = @SVector[0.0, cone.length/2, 0.0]
        R_obj = Modia3D.NullRotation
    else
        r_obj = @SVector[0.0, 0.0, cone.length/2]
        R_obj = @SMatrix[1 0 0; 0 0 1; 0 -1 0]
    end
    name = String(Modia3D.fullName(obj)) * ".geometry"
    shape = (; name=name, uuid=name2uuid(name), type="CylinderBufferGeometry", radiusBottom=cone.diameter/2, radiusTop=cone.topDiameter/2, height=cone.length, radialSegments=32, heightSegments=1)
    material = printVisuMaterialToJSON(obj, obj.visualMaterial)
    objectInfo = getObjectInfo(obj, shape, material, initPos, initRot, r_obj=r_obj, R_obj=R_obj)
    printInfoToFile(object, shapes, shape, material, objectInfo)
    return (r_obj, R_obj)
end

function printInfoToFile(object, shapes, shape, material, objectInfo)
    push!(shapes.geometries, shape)
    push!(shapes.materials, material)
    push!(object.children, objectInfo)
end

function exportObject(object, shapes, obj::Modia3D.Composition.Object3D, fileMesh::Modia3D.Shapes.FileMesh, initPos, initRot)
    r_obj = Modia3D.ZeroVector3D
    R_obj = Modia3D.NullRotation
    filenameOld = fileMesh.filename
    filename = filenameOld
    (head,ext) = splitext(filenameOld)
    if ext == ".json"
        filename = filenameOld
    elseif ext == ".obj"
        filename = head * ".json"
        if !isfile(filename)
            @warn("Please convert $filenameOld into a .json file.")
            return nothing
        end
    else
        @warn("Please convert $filenameOld into a .json file.")
        return nothing
    end

    io = open(filename, "r")
    meshObject = JSON.parse(io)
    close(io)

    push!(shapes.geometries, meshObject["geometries"]...)
    push!(shapes.materials, meshObject["materials"]...)
    objectInfo = getObjectInfoMesh(obj, initPos, initRot, fileMesh.scaleFactor, R_obj, meshObject["object"])
    push!(object.children, objectInfo)

    return (r_obj, R_obj)
end

function getObjectInfoMesh(obj, initPos, initRot, scale, R_obj, meshInfo)
    push!(meshInfo, "name" => Modia3D.fullName(obj) )
    delete!(meshInfo, "matrix")
    push!(meshInfo, "position" => initPos )
    push!(meshInfo, "rotation" => Modia3D.rot123fromR(R_obj*initRot) )
    push!(meshInfo, "scale" => scale )

    return meshInfo
end

function getObjectInfo(obj, shape, material, initPos, initRot; r_obj=Modia3D.ZeroVector3D, R_obj=Modia3D.NullRotation, scale=ones(3))
    name = String(Modia3D.fullName(obj))
    return (; name=name, uuid=name2uuid(name), type=:Mesh, geometry=get(shape, :uuid, nothing), material=get(material, :uuid, ""), position=initPos+initRot'*r_obj, rotation=Modia3D.rot123fromR(R_obj*initRot), scale=scale )
end

printVisuMaterialToJSON(obj, visuMaterial) = nothing
function printVisuMaterialToJSON(obj, visuMaterial::Modia3D.Shapes.VisualMaterial)
    name = String(Modia3D.fullName(obj)) * ".material"

    transparent = visuMaterial.transparency > 0

    material = (; name=name, uuid=name2uuid(name), type="MeshPhongMaterial", color=((visuMaterial.color[1]*256 + visuMaterial.color[2])*256 + visuMaterial.color[3]), opacity=1-visuMaterial.transparency, transparent=transparent,
    shininess=visuMaterial.shininess*100, side="DoubleSide")

    return material
end

function printObjectToJSON(object, shapes, obj; initPos=obj.r_abs, initRot=obj.R_abs)
    r_obj = Modia3D.ZeroVector3D
    R_obj = Modia3D.NullRotation
    shapeKind = obj.shapeKind

    if shapeKind == Modia3D.SphereKind
        sphere::Modia3D.Shapes.Sphere = obj.shape
        (r_obj, R_obj) = exportObject(object, shapes, obj, sphere, initPos, initRot)

    elseif shapeKind == Modia3D.EllipsoidKind
        ellipsoid::Modia3D.Shapes.Ellipsoid = obj.shape
        (r_obj, R_obj) = exportObject(object, shapes, obj, ellipsoid, initPos, initRot)

    elseif shapeKind == Modia3D.BoxKind
        box::Modia3D.Shapes.Box = obj.shape
        (r_obj, R_obj) = exportObject(object, shapes, obj, box, initPos, initRot)

    elseif shapeKind == Modia3D.CylinderKind
        cylinder::Modia3D.Shapes.Cylinder = obj.shape
        (r_obj, R_obj) = exportObject(object, shapes, obj, cylinder, initPos, initRot)

    elseif shapeKind == Modia3D.ConeKind
        cone::Modia3D.Shapes.Cone = obj.shape
        (r_obj, R_obj) = exportObject(object, shapes, obj, cone, initPos, initRot)

    elseif shapeKind == Modia3D.FileMeshKind
        fileMesh::Modia3D.Shapes.FileMesh = obj.shape
        (r_obj, R_obj) = exportObject(object, shapes, obj, fileMesh, initPos, initRot)

    else
        @warn("$shapeKind is not supported by animation export.")
        return (nothing, nothing)
    end
    return (r_obj, R_obj)
end

function addLights!(object, uuid4, position)
    name = "Directional Light"
    uuid = string(UUIDs.uuid5(UUIDs.UUID(uuid4), name))
    directionalLight = (; name=name, uuid=uuid, type="DirectionalLight", position=position)
    push!(object.children, directionalLight)

    name="Ambient Light"
    uuid = string(UUIDs.uuid5(UUIDs.UUID(uuid4), name))
    ambientLight = (; name=name, uuid=uuid, type="AmbientLight", intensity=0.5)
    push!(object.children, ambientLight)
end

function addCameras!(object, uuid4, position, orientation)
    radius = LinearAlgebra.norm(position)

    name = "Perspective Camera"
    uuid = string(UUIDs.uuid5(UUIDs.UUID(uuid4), name))
    camera = (; name=name, uuid=uuid, type="PerspectiveCamera", position=position, rotation=Modia3D.rot123fromR(SMatrix{3,3,Float64}(orientation)))
    push!(object.children, camera)

    name="Orthographic Camera XY"
    uuid = string(UUIDs.uuid5(UUIDs.UUID(uuid4), name))
    camera = (; name=name, uuid=uuid, type="OrthographicCamera", position=[0, 0, radius], rotation=[0, 0, 0], bottom=-radius/2, top=radius/2, left=-radius/2, right=radius/2)
    push!(object.children, camera)

    name="Orthographic Camera YZ"
    uuid = string(UUIDs.uuid5(UUIDs.UUID(uuid4), name))
    camera = (; name=name, uuid=uuid, type="OrthographicCamera", position=[radius, 0, 0], rotation=[0, deg2rad(90), deg2rad(90)], bottom=-radius/2, top=radius/2, left=-radius/2, right=radius/2)
    push!(object.children, camera)

    name="Orthographic Camera ZX"
    uuid = string(UUIDs.uuid5(UUIDs.UUID(uuid4), name))
    camera = (; name=name, uuid=uuid, type="OrthographicCamera", position=[0, radius, 0], rotation=[deg2rad(-90), 0, deg2rad(-90)], bottom=-radius/2, top=radius/2, left=-radius/2, right=radius/2)
    push!(object.children, camera)
end

function createAnimationPositionTrack(tracks, animation, obj, iobj, r_obj)
    keys = []
    for istep in 1:length(animation)
        r_all = animation[istep].objectData[iobj].position + (Modia3D.from_q(SVector{4,Float64}(animation[istep].objectData[iobj].quaternion)))'*r_obj
        push!(keys, (; time=animation[istep].time, value=r_all) )
    end
    name = string(Modia3D.fullName(obj), ".position")
    push!(tracks, (; name=name, uuid=name2uuid(name), type="vector3", keys ) )
end

function createAnimationQuaternionTrack(tracks, animation, obj, iobj, R_obj)
    keys = []
    for istep in 1:length(animation)
        R_all = Modia3D.from_R(R_obj*Modia3D.from_q(SVector{4,Float64}(animation[istep].objectData[iobj].quaternion)) )
        push!(keys, (; time=animation[istep].time, value=R_all) )
    end
    name = string(Modia3D.fullName(obj), ".quaternion")
    push!(tracks, (; name=name, uuid=name2uuid(name), type="quaternion", keys ) )

end

function createAnimationQuaternionTrack(object, animation, obj, iobj, R_obj::Nothing)
    return nothing
end

function exportAnimation(scene)
    animationFile = scene.options.animationFile
    allVisuElements = scene.allVisuElements
    animation = scene.animation
    if length(allVisuElements) > 0 && !isnothing(animationFile)
        (head,ext) = splitext(animationFile)
        if ext != ".json"
            @warn("The path of an animationFile=$animationFile must end with .json.")
            return nothing
        end
        print("Export animation to $animationFile ... ")
        shapes = (; geometries=[], materials=[])
        metadata = (; generator = "Modia3D", type = "Object")
        name = scene.name
        uuid = name2uuid(name)
        object = (; name, uuid, animations=[uuid], children=[])

        options = scene.options
        cameraUpDir = 2
        sceneUpDir = Modia3D.upwardsDirection(options.gravityField)
        (r_Light, R_Light) = Modia3D.cameraPosition(options.lightDistance, options.lightLongitude, options.lightLatitude, cameraUpDir, sceneUpDir)
        addLights!(object, uuid, r_Light)
        (r_Camera, R_Camera) = Modia3D.cameraPosition(options.cameraDistance, options.cameraLongitude, options.cameraLatitude, cameraUpDir, sceneUpDir)
        addCameras!(object, uuid, r_Camera, R_Camera)

        if !isnothing(animation) && length(animation) != 0
            iobj = 0
            tracks = []
            for obj in allVisuElements
                iobj = iobj + 1
                (r_obj, R_obj) = printObjectToJSON(object, shapes, obj, initPos=animation[1].objectData[iobj].position, initRot=Modia3D.from_q(SVector{4,Float64}(animation[1].objectData[iobj].quaternion)) )
                if !isnothing(R_obj)
                    createAnimationPositionTrack(tracks, animation, obj, iobj, r_obj)
                    createAnimationQuaternionTrack(tracks, animation, obj, iobj, R_obj)
                end
            end
            animations = [(; name="Simulation", uuid=uuid, tracks)]
            scene = (; metadata, shapes.geometries, shapes.materials, object, animations)
        else
            for obj in allVisuElements
                (r_obj, R_obj) = printObjectToJSON(object, shapes, obj)
            end
            scene = (; metadata, shapes.geometries, shapes.materials, object)
        end

        io = open(animationFile, "w")
        JSON.print(io, scene, 0)  # 3rd parameter: nothing = not human readable; 0 = no indentation; n = number of indentation spaces
        close(io)

        println("done.")
    end
end


text:
"""
    - nonnegative orthant (no)
        x >= 0

    A Semismooth Newton Method for Fast, Generic Convex Programming
        https://arxiv.org/abs/1705.00772
"""

# nonnegative orthant cone
function κ_no(z)
    max.(0.0, z)
end

# nonnegative orthant cone Jacobian
function Jκ_no(z)
    p = zero(z)
    for (i, pp) in enumerate(z)
        if pp >= 0.0
            p[i] = 1.0
        end
    end
    return Diagonal(p)
end


text:
__precompile__()

module FullRankModels

import LowRankModels: ConvergenceHistory, update_ch!, fit!

include("regularizers.jl")
include("gfrm.jl")
include("evaluate_fit.jl")

include("algorithms/prisma.jl")
# include("algorithms/fw_sketch.jl")
include("algorithms/fw_sketch_multithread.jl")

include("utilities.jl")
include("cross_validate.jl")

end # module


text:
struct PSOSolution{M,T}

    Swarm::Swarm # Final Swarm
    fbest::T
    xbest::SVector{M,T}

    Iters::Int
    Time::T
    ExitFlag::Int

    # If a hybrid function is also used. (May decide to remove this in the future)
    HybridFcnRes::Union{Nothing, Optim.OptimizationResults}
    fbestPSO::Union{Nothing, T}
    xbestPSO::Union{Nothing, SVector{M,T}}

end

function PSOSolution(S::Swarm, fbest::T, xbest::SVector{M,T},
                     Iters::Int, Time::T, ExitFlag::Int) where {M,T}
    PSOSolution(S, fbest, xbest, Iters, Time, ExitFlag, nothing, nothing, nothing)
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
using MLStyle

macro cond(cases)
    @match cases begin
        quote
        $(cases...)
        end =>
        let default = Expr(:call, throw, "None of the branches have satisfied conditions.")
            foldr(cases, init = default) do case, last
                last_lnode = @__LINE__
                @match case begin
                    ::LineNumberNode => begin
                        last_lnode = case
                        Expr(:block, case, last)
                    end
                    :(_ => $b) => b
                    :($a => $b) => Expr(:if, a, b, last)
                    _ => throw("Invalid syntax at $last_lnode.")
                end
            end
        end |> esc
    end
end

function check_more_than1(x :: Int)
    @cond begin
        x > 1 => true
        _ => false
    end
end

@assert check_more_than1(1) === false
@assert check_more_than1(0) === false
@assert check_more_than1(2) === true

max_stars_repo_name: cherichy/MLStyle-Playground

max_issues_repo_name: cherichy/MLStyle-Playground

max_forks_repo_name: cherichy/MLStyle-Playground


text:
module two_tanks #jl

using ReachabilityAnalysis, Symbolics

function _two_tanks()
    var = @variables x1 x2

    HA = LightAutomaton(4)

    B = hcat([1.0, 1])
    U = Interval(-0.1, 0.1)

    ## off/off
    X = HPolyhedron([x1 ≥ -1, x2 ≤ 1], var)
    A = [-1.0 0; 1 0]
    c = [-2.0, 0]
    m1 = @system(x' = A*x + B*u + c, x ∈ X, u ∈ U)

    ## on/off
    X = HPolyhedron([x2 ≤ 1], var)
    A = [-1.0 0; 1 0]
    c = [3.0, 0]
    m2 = @system(x' = A*x + B*u + c, x ∈ X, u ∈ U)

    ## off/on
    X = HPolyhedron([x1 ≥ -1, x2 ≥ 0], var)
    A = [-1.0 0; 1 -1]
    c = [-2.0, -5]
    m3 = @system(x' = A*x + B*u + c, x ∈ X, u ∈ U)

    ## on/on
    X = HPolyhedron([x1 ≤ -1, x2 ≥ 0], var)
    A = [-1.0 0; 1 -1]
    c = [3.0, -5]
    m4 = @system(x' = A*x + B*u + c, x ∈ X, u ∈ U)

    ## transition off/off -> on/off
    add_transition!(HA, 1, 2, 1)
    G = Hyperplane(x1 == -1, var)
    t1 = ConstrainedIdentityMap(2, G)

    ## transition off/off -> off/on
    add_transition!(HA, 1, 3, 2)
    G = Hyperplane(x2 == 1, var)
    t2 = ConstrainedIdentityMap(2, G)

    ## transition on/off -> off/on
    add_transition!(HA, 2, 3, 3)
    G = Hyperplane(x2 == 1, var)
    t3 = ConstrainedIdentityMap(2, G)

    ## transition off/on -> off/off
    add_transition!(HA, 3, 1, 4)
    G = Hyperplane(x2 == 0, var)
    t4 = ConstrainedIdentityMap(2, G)

    ## transition off/on -> on/on
    add_transition!(HA, 3, 4, 5)
    G = Hyperplane(x1 == -1, var)
    t5 = ConstrainedIdentityMap(2, G)

    ## transition on/on -> on/off
    add_transition!(HA, 4, 2, 6)
    G = Hyperplane(x2 == 0, var)
    t6 = ConstrainedIdentityMap(2, G)

    ## transition on/on -> off/on
    add_transition!(HA, 4, 3, 7)
    G = Hyperplane(x1 == 1, var)
    t7 = ConstrainedIdentityMap(2, G)

    ## hybrid system
    S = HybridSystem(HA, [m1, m2, m3, m4], [t1, t2, t3, t4, t5, t6, t7])

    return S
end

function model(X0)
    H = _two_tanks()
    return IVP(H, X0)
end

end #jl

# ## References

# https://ths.rwth-aachen.de/research/projects/hypro/two-tank/

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# Find the linear regression of the form: f = (x) -> ∑βixi+ϵi

# Load CSV data
data = readcsv("gasoline.csv")

# Dissect data
x = data[:, 2:4]
y = data[:, 6]

# Call linreg
coefs = linreg(x, y);

max_stars_repo_name: neilpanchal/Julia-Algorithms

max_issues_repo_name: neilpanchal/Julia-Algorithms

max_forks_repo_name: neilpanchal/Julia-Algorithms


text:
using SubstitutionTilings
using SubstitutionTilings.Chair
using SubstitutionTilings.Collaring
using Test

using Luxor

cross = [
    chair(0,0,0),
    chair(2,-2,-2),
    chair(1,0,-2),
    chair(3,-2,0),
]

tiling = substitute(chair_system, [chair(0,0,0)], 8, Chair.in_bounds, (w=80, h=80))
tiling_dict = Dict(tiling)
@time empirical_frequency(cross, tiling_dict)
@draw begin
    colors = ["#3CD0E6", "#CA7EE6", "#E6873C", "#B4E647"]
    sc = 10
    tiling = substitute(chair_system, [chair(0,0,0)], 8, Chair.in_bounds, (w=80, h=80))
    collar = substitute(chair_system, collar_class(tiling, ChairElem(0,1,1)), 0)

    setline(0.1)

    for tile in collar
        draw(tile, sc, colors[1 + tile[1].angle % 4], :fill)
        draw(tile, sc, "black", :stroke)
    end
end 800 800

frequency(chair_system, collar, cross, 4)
@time sizeof(substitute(chair_system, [chair(0,0,0)], 20, Chair.in_bounds, (w=800, h=800)))


text:
#=##############################################################################
# DESCRIPTION
    Unsteady simulation of eVTOL transition maneuver on tilt-wing, tandem,
    distributed propulsion aircraft. Geometry of a modified Vahana aircraft with
    tilt rotors.

    REFERENCES
    * Vahana geometry: Droandi, G., Syal, M., and Bower, G., “Tiltwing Multi-Rotor
    Aerodynamic Modeling in Hover, Transition and Cruise Flight Conditions,” AHS
    International 74th Annual Forum & Technology Display, 2018, p. 2018.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Feb 2021
  * License   : MIT
=###############################################################################

# ------------ MODULES ---------------------------------------------------------
using Revise
import QuadGK
import LinearAlgebra: I
import FLOWVPM
import FLOWUnsteady

uns = FLOWUnsteady
vpm = FLOWVPM
vlm = uns.vlm
gt = uns.gt

# ------------ GLOBAL VARIABLES ------------------------------------------------
# Default path where to save data
# extdrive_path = "/media/edoalvar/Samsung_T5/simulationdata202102/"
extdrive_path = "/media/flowlab/Storage/ealvarez/simulations/"

# Default data path where to find rotor and airfoil data
data_path = uns.def_data_path

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["geometry", "kinematics", "monitor",
                                                   "misc", "postprocessing"]
    include("vahana2_"*header_name*".jl")
end

# ------------ DRIVERS ---------------------------------------------------------

function run_simulation_vahana(;    save_path=extdrive_path*"vahana2_sim19",
                                    prompt=true,
                                    run_name="vahana2",
                                    verbose=true, v_lvl=1)

    # ----------------- PARAMETERS ---------------------------------------------

    # Geometry options
    n_factor = 1                            # Refinement factor
    add_rotors = true                       # Whether to include rotors

    ## 72 steps per rev settings
    Vcruise = 0.25 * 125*0.44704            # Cruise speed
    Vinf(x,t) = 1e-5*[1,0,-1]               # (m/s) freestream velocity, if 0 the simulation might crash
    RPMh_w = 600.0                          # RPM of main wing rotors in hover
    telapsed = 30.0                         # Total time to perform maneuver
    nsteps = 4*5400                         # Time steps for complete maneuver
    lambda = 2.125                          # Target minimum core overlap
    # p_per_step = 2                          # Particle sheds per time step
    vlm_rlx = 0.2                           # VLM relaxation (deactivated with -1)

    p_per_step = 2*2                          # Particle sheds per time step

    # # Maneuver to perform
    # ## 18 steps per rev settings
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # Vinf(x,t) = 1e-5*[1,0,-1]               # (m/s) freestream velocity, if 0 the simulation might crash
    # RPMh_w = 600.0                          # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 5400                           # Time steps
    # lambda = 2.125                          # Target minimum core overlap
    # p_per_step = 4                          # Particle sheds per time step
    # vlm_rlx = 0.75                          # VLM relaxation (deactivated with -1)

    # # Maneuver to perform
    # Vcruise = 0.125 * 125*0.44704         # Cruise speed
    # RPMh_w = 600                          # RPM of main wing rotors in hover
    # telapsed = 60.0                       # Total time to perform maneuver
    # nsteps = 9000                         # Time steps

    # # Maneuver to perform
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # Vinf(x,t) = 1e-5*[1,0,-1]               # (m/s) freestream velocity, if 0 the simulation might crash
    # # RPMh_w = 200                          # RPM of main wing rotors in hover
    # # RPMh_w = 20
    # RPMh_w = 100
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 1500                           # Time steps
    # # nsteps = 100

    dt = telapsed/nsteps

    # Simulation options
    rho = 1.225                             # (kg/m^3) air density
    mu = 1.81e-5                            # Air dynamic viscosity

    # Solver options
    R = 0.75                                # (m) Reference blade radius
    # lambda = 2.125                          # Target minimum core overlap
    # p_per_step = 4                          # Particle sheds per time step
    overwrite_sigma = lambda * (2*pi*RPMh_w/60*R + Vcruise)*dt / p_per_step
    surf_sigma = R/10                       # Surface regularization
    # vlm_sigma = R/25                      # VLM regularization
    # vlm_sigma = R
    vlm_sigma = -1
    # vlm_sigma = surf_sigma # MAKE THIS -1 if unstable
    # vlm_rlx = 0.75                          # VLM relaxation (deactivated with -1)
    # shed_unsteady = false                   # Shed unsteady-loading particles
    shed_unsteady = true
    # VehicleType = uns.QVLMVehicle           # Type of vehicle to generate
    VehicleType = uns.UVLMVehicle           # Type of vehicle to generate

    vpm_formulation = vpm.formulation_sphere # VPM formulation
    # vpm_formulation = vpm.formulation_classic
    # vpm_sgsmodel    = vpm.sgs_stretching1_fmm
    vpm_sgsmodel    = vpm.generate_sgs_directionfiltered(vpm.generate_sgs_lowfiltered(vpm.sgs_stretching1_fmm))
    # vpm_sgsscaling(args...) = 0.75
    vpm_sgsscaling(args...) = 1.0
    # vpm_sgsscaling = vpm.sgs_scaling_none
    # vpm_relaxation  = vpm.norelaxation
    # vpm_relaxation  = vpm.pedrizzetti
    vpm_relaxation  = vpm.correctedpedrizzetti

    # ----------------- MANEUVER DEFINITION ------------------------------------
    # Generate maneuver
    maneuver = generate_maneuver_vahana(; add_rotors=add_rotors)

    # Plot maneuver path and controls
    gt.create_path(save_path, prompt)
    uns.plot_maneuver(maneuver; tstages=[0.2, 0.3, 0.5, 0.6], save_path=save_path)

    # ----------------- GEOMETRY GENERATION ------------------------------------
    # Generate geometry
    (vehicle, grounds) = generate_geometry_vahana(; n_factor=n_factor,
                                                    xfoil=false,
                                                    data_path=data_path,
                                                    run_name=run_name,
                                                    add_rotors=add_rotors,
                                                    VehicleType=VehicleType)

    # Move landing pad to landing area
    gt.lintransform!(grounds[2], Array(1.0I, 3, 3), Vcruise*telapsed*[-0.25, 0, -0.0025])

    # Save ground
    for (i, ground) in enumerate(grounds)
        gt.save(ground, run_name*"_Ground$i"; path=save_path)
    end

    # ----------------- SIMULATION SETUP ---------------------------------------
    Vref = Vcruise
    RPMref = RPMh_w
    ttot = telapsed
    # max_particles = ceil(Int, (nsteps+2)*(2*vlm.get_m(vehicle.vlm_system)+1)*p_per_step)
    max_particles = 2*400000

    Vinit = Vref*maneuver.Vvehicle(0)       # (m/s) initial vehicle velocity
                                            # (rad/s) initial vehicle angular velocity
    Winit = pi/180 * (maneuver.anglevehicle(0+1e-12)-
                                          maneuver.anglevehicle(0))/(ttot*1e-12)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                                    Vinit=Vinit, Winit=Winit)

    # ----------------- SIMULATION MONITOR -------------------------------------
    monitor_vahana = generate_monitor_vahana(vehicle, rho, RPMref, nsteps, save_path, Vinf)
    monitor_vpm = uns.generate_monitor_enstrophy(; save_path=save_path)

    monitor(args...; optargs...) = monitor_vahana(args...; optargs...) || monitor_vpm(args...; optargs...)

    # ----------------- WAKE TREATMENT FUNCTION --------------------------------
    remove_particles(args...; optargs...) = (remove_particles_lowstrength(args...; optargs...)
                                          || remove_particles_sphere(args...; optargs...))

    # ----------------- RUNTIME FUNCTION ---------------------------------------
    runtime_function(args...; optargs...) = remove_particles(args...; optargs...) || monitor(args...; optargs...)

    # ----------------- RUN SIMULATION -----------------------------------------
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      rho=rho,
                                      mu=mu,
                                      # SOLVERS OPTIONS
                                      vpm_formulation=vpm_formulation,
                                      vpm_sgsmodel=vpm_sgsmodel,
                                      vpm_sgsscaling=vpm_sgsscaling,
                                      vpm_relaxation=vpm_relaxation,
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      vlm_rlx=vlm_rlx,
                                      surf_sigma=surf_sigma,
                                      max_particles=max_particles,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=runtime_function,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      create_savepath=false,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose, v_lvl=v_lvl,
                                      save_code=splitdir(@__FILE__)[1]
                                      )

    return simulation, pfield
end

function visualize_maneuver_vahana(; save_path=extdrive_path*"vahana2_maneuver04/",
                                        prompt=true,
                                        run_name="vahana2",
                                        verbose=true, v_lvl=0,
                                        paraview=true,
                                        optargs...)

    # Maneuver to perform
    # ## 72 steps per rev settings
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # RPMh_w = 600.0                          # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 21600                          # Time steps

    ## 4 steps per rev settings
    Vcruise = 0.25 * 125*0.44704            # Cruise speed
    RPMh_w = 600.0                          # RPM of main wing rotors in hover
    telapsed = 30.0                         # Total time to perform maneuver
    nsteps = Int(21600/18)                  # Time steps

    # # Maneuver to perform
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # RPMh_w = 400.0                          # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 100                            # Time steps

    # # Maneuver to perform
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # RPMh_w = 1200.0                         # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 9000                           # Time steps

    # # Maneuver to perform
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # RPMh_w = 1200.0                         # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 27000                           # Time steps

    # # Maneuver to perform
    # Vcruise = 0.125 * 125*0.44704            # Cruise speed
    # RPMh_w = 600.0                          # RPM of main wing rotors in hover
    # telapsed = 60.0                         # Total time to perform maneuver
    # nsteps = 9000                           # Time steps

    dt = telapsed/nsteps

    # Generate maneuver
    maneuver = generate_maneuver_vahana(; optargs...)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; tstages=[0.2, 0.3, 0.5, 0.6])

    # Generate geometry
    (vehicle, grounds) = generate_geometry_vahana(; n_factor=1,
                                                    xfoil=false,
                                                    data_path=data_path,
                                                    run_name=run_name,
                                                    optargs...)

    # Simulation setup
    Vref = Vcruise
    RPMref = RPMh_w
    ttot = telapsed

    Vinit = Vref*maneuver.Vvehicle(0)       # (m/s) initial vehicle velocity
                                            # (rad/s) initial vehicle angular velocity
    Winit = pi/180 * (maneuver.anglevehicle(0+1e-12)-
                                          maneuver.anglevehicle(0))/(ttot*1e-12)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                                    Vinit=Vinit, Winit=Winit)

    save_vtk_optsargs = [(:save_horseshoes, false)]

    # Visualize maneuver
    strn = uns.visualize_kinematics(simulation, nsteps, save_path;
                                    run_name=run_name,
                                    save_vtk_optsargs=save_vtk_optsargs,
                                    prompt=prompt, verbose=verbose, v_lvl=v_lvl,
                                    paraview=false
                                    )

    # Move landing pad to landing area
    gt.lintransform!(grounds[2], Array(1.0I, 3, 3), Vcruise*telapsed*[-0.25, 0, -0.0025])

    # Save ground
    for (i, ground) in enumerate(grounds)
        gt.save(ground, run_name*"_Ground$i"; path=save_path)
        strn *= run_name*"_Ground$i.vtk;"
    end

    # Call paraview
    if paraview
        run(`paraview --data="$save_path/$strn"`)
    end

    return strn
end

"""
    Generates geometry of Vahana aircraft, saves it as vtk files, and calls
    Paraview visualizing the VTK geometry.
"""
function visualize_geometry_vahana(; save_path=extdrive_path*"vahana2_geometry01/",
                                     prompt=true, run_name="vahana2", optargs...)

    (vehicle, grounds) = generate_geometry_vahana(; n_factor=1,
                                                     xfoil=false,
                                                     data_path=data_path,
                                                     run_name=run_name,
                                                     optargs...)

    # Setup dummy freestream and RPMs for horseshoe visualization
    Vinf(x,t) = [1,0,0]
    vlm.setVinf(vehicle.system, Vinf)

    for rotor_system in vehicle.rotor_systems
        for rotor in rotor_system
            vlm.setRPM(rotor, 6000)
        end
    end

    # Create save path
    gt.create_path(save_path, prompt)

    # Save vehicle
    strn = uns.save_vtk(vehicle, run_name; path=save_path, save_horseshoes=false)

    # Save ground
    for (i, ground) in enumerate(grounds)
        gt.save(ground, run_name*"_Ground$i"; path=save_path)
        strn *= run_name*"_Ground$i.vtk;"
    end

    # Call Paraview
    run(`paraview --data="$save_path/$strn"`)
end

function vahana2_postprocess(; ite=1)

    range = collect(4:4:6480)                   # Time steps to read
    nchunks = 4                                 # Number of chunks
    chunk = ceil(Int, length(range)/nchunks)    # Chunk length
    # ite = 1                                     # This chunk to iterate over

    nums = range[chunk*(ite-1) + 1 : (chunk*ite > length(range) ? length(range) : chunk*ite)]

    read_path = extdrive_path*"vahana2_sim16/"
    save_path = extdrive_path*"vahana2_sim16-fdom00-$(ite)/"

    postprocessing_fluiddomain(read_path, save_path, nums;
                                        # PROCESSING OPTIONS
                                        fdx=1/80,                # Scaling of cell length
                                        static_particles=false,  # Add static particles
                                        # OUTPUT OPTIONS
                                        run_name="vahana2",
                                        prompt=true, paraview=false)
end


text:
# Let's try incommensurate system
using FdeSolver
using Plots

# Intputs
tSpan = [0, 50]
h = 2^(-5)
y0 = [1.2, 2.8]
β = [0.8, 0.7]
A = 1 ; B = 3
par = [ A , B ]
# Definition of the System
function F(t, y, par)

    F1 = par[1] .- (par[2]+1) .* y[1] .+ y[1] .^ 2 .* y[2]
    F2 = par[2] .* y[1] .- y[1] .^ 2 .* y[2]

    return [F1, F2]

end

t, Yapp = FDEsolver(F, tSpan, y0, β, par)
plot(t, Yapp)


text:
@vectorize_1arg Real HaarScaling

function HaarScaling(xi::AbstractArray, J::Int, k::Int)
	y = Array{Float64}(size(xi))
	for idx in eachindex(xi)
		@inbounds y[idx] = HaarScaling( xi[idx], J, k )
	end
	return y
end


text:
using XGPaint
using Healpix

halo_pos, halo_mass = read_halo_catalog_hdf5(
    "/home/zequnl/Data/websky_halos-light.hdf5")
## Load halos from HDF5 files, establish a CIB model and cosmology
cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0)
radio_model = Radio_Sehgal2009{Float32}()

##

@time begin
    sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass);
end;

##
@time begin
m = Map{Float64,RingOrder}(radio_model.nside)
flux_I, flux_II, redshift_I, redshift_II = paint!(
    m, 151f6, radio_model, sources, return_fluxes=true)
end
##

using DelimitedFiles
using Unitful
using UnitfulAstro
using Cosmology
using PyPlot
using StatsBase

figure(figsize=(7,7))

sehgal_red = readdlm("examples/data/sehgal_figure8_red.txt", ',', Float64, '\n')
sehgal_blue = readdlm("examples/data/sehgal_figure8_blue.txt", ',', Float64, '\n')
sehgal_green = readdlm("examples/data/sehgal_figure8_green.txt", ',', Float64, '\n')

vol = ustrip(u"Mpc^3",comoving_volume(cosmo, 0.3))
bins = (range(23, stop=29, step=0.5));
mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0
count_I = fit(Histogram, log10.(sources.L_I_151[redshift_I .< 0.3]), bins).weights

count_II = fit(Histogram, log10.(sources.L_II_151[redshift_II .< 0.3]), bins).weights

plot(mids, (count_I .+ count_II) ./ diff((bins)) ./ vol, label="websky total" )
plot(mids, count_I ./ diff((bins)) ./ vol, label="websky FR I" )
plot(mids, count_II ./ diff((bins)) ./ vol, label="websky FR II" )

plot(log10.(sehgal_red[:, 1]), sehgal_red[:, 2], "rp",
    label = "sehgal total", alpha = 0.6)
plot(log10.(sehgal_green[:, 1]), sehgal_green[:, 2], color = "green",
    "^", label = "sehgal FR I", alpha = 0.6)
plot(log10.(sehgal_blue[:, 1]), sehgal_blue[:, 2], color = "blue", "s",
    label = "sehgal FR II", alpha = 0.6)

ylabel("counts / log bin / Mpc\$^3\$")
xlabel("log L [W/Hz]")
yscale("log")

legend()
tight_layout()
savefig("sehgal_LF.pdf")
gcf()

##

plt.figure(figsize=(6,6))

bins = (range(-4.1, stop=2.5, step=0.2));
mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0

flux_I, flux_II, redshift_I, redshift_II = paint!(
    m, 1.4f9, radio_model, sources, return_fluxes=true)
count_I = fit(Histogram, log10.(flux_I), bins).weights
count_II = fit(Histogram, log10.(flux_II), bins).weights

Sehgal = readdlm("examples/data/sehgal_fig9_1.4GHz.txt", ',', Float64, '\n');
plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 1.4 GHz")

plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total")
plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I")
plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II")

yscale("log")
xscale("log")
xlabel("\$ S \$ [Jy]")
ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]")
xticks(10 .^ [-3.,-2.,-1.,0.,1.])
legend()
tight_layout()
savefig("sehgal_1.4GHz.pdf")
gcf()

##
plt.figure(figsize=(6,6))

bins = (range(-4.1, stop=2.5, step=0.2));
mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0

flux_I, flux_II, redshift_I, redshift_II = paint!(
    m, 33f9, radio_model, sources, return_fluxes=true)
count_I = fit(Histogram, log10.(flux_I), bins).weights
count_II = fit(Histogram, log10.(flux_II), bins).weights

Sehgal = readdlm("examples/data/sehgal_fig9_33GHz.txt", ',', Float64, '\n');
plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 33 GHz")

plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total")
plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I")
plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II")

yscale("log")
xscale("log")
xlabel("\$ S \$ [Jy]")
ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]")
xticks(10 .^ [-3.,-2.,-1.,0.,1.])
legend()
tight_layout()
savefig("sehgal_33GHz.pdf")
gcf()
##
plt.figure(figsize=(6,6))

bins = (range(-4.1, stop=2.5, step=0.2));
mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0

flux_I, flux_II, redshift_I, redshift_II = paint!(
    m, 95f9, radio_model, sources, return_fluxes=true)
count_I = fit(Histogram, log10.(flux_I), bins).weights
count_II = fit(Histogram, log10.(flux_II), bins).weights

Sehgal = readdlm("examples/data/sehgal_fig9_95GHz.txt", ',', Float64, '\n');
plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 95 GHz")

plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total")
plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I")
plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II")

yscale("log")
xscale("log")
xlabel("\$ S \$ [Jy]")
ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]")
xticks(10 .^ [-3.,-2.,-1.,0.,1.])
legend()
tight_layout()
savefig("sehgal_95GHz.pdf")
gcf()

##
plt.figure(figsize=(6,6))

bins = (range(-4.1, stop=2.5, step=0.2));
mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0

flux_I, flux_II, redshift_I, redshift_II = paint!(
    m, 145f9, radio_model, sources, return_fluxes=true)
count_I = fit(Histogram, log10.(flux_I), bins).weights
count_II = fit(Histogram, log10.(flux_II), bins).weights

Sehgal = readdlm("examples/data/sehgal_fig9_145GHz.txt", ',', Float64, '\n');
plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 145 GHz")

plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total")
plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I")
plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II")

yscale("log")
xscale("log")
xlabel("\$ S \$ [Jy]")
ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]")
xticks(10 .^ [-3.,-2.,-1.,0.,1.])
legend()
tight_layout()
savefig("sehgal_145GHz.pdf")
gcf()

##


text:
# ModiaLang interface to Modia3D
Object3D(        ; kwargs...) = Par(; _constructor = :(Modia3D.Object3D), _path = true, kwargs...)
Scene(           ; kwargs...) = Par(; _constructor = :(Modia3D.SceneOptions)          , kwargs...)
Visual(          ; kwargs...) = Par(; _constructor = :(Modia3D.Visual)                , kwargs...)
Solid(           ; kwargs...) = Par(; _constructor = :(Modia3D.Solid)                 , kwargs...)
Box(             ; kwargs...) = Par(; _constructor = :(Modia3D.Box)                   , kwargs...)
Beam(            ; kwargs...) = Par(; _constructor = :(Modia3D.Beam)                  , kwargs...)
Cylinder(        ; kwargs...) = Par(; _constructor = :(Modia3D.Cylinder)              , kwargs...)
Sphere(          ; kwargs...) = Par(; _constructor = :(Modia3D.Sphere)                , kwargs...)
Ellipsoid(       ; kwargs...) = Par(; _constructor = :(Modia3D.Ellipsoid)             , kwargs...)
Capsule(         ; kwargs...) = Par(; _constructor = :(Modia3D.Capsule)               , kwargs...)
Cone(            ; kwargs...) = Par(; _constructor = :(Modia3D.Cone)                  , kwargs...)
SpringShape(     ; kwargs...) = Par(; _constructor = :(Modia3D.Spring)                , kwargs...)
GearWheel(       ; kwargs...) = Par(; _constructor = :(Modia3D.GearWheel)             , kwargs...)
Grid(            ; kwargs...) = Par(; _constructor = :(Modia3D.Grid)                  , kwargs...)
VisualMaterial(  ; kwargs...) = Par(; _constructor = :(Modia3D.VisualMaterial)        , kwargs...)
MassProperties(  ; kwargs...) = Par(; _constructor = :(Modia3D.MassProperties)        , kwargs...)
CoordinateSystem(; kwargs...) = Par(; _constructor = :(Modia3D.CoordinateSystem)      , kwargs...)
FileMesh(        ; kwargs...) = Par(; _constructor = :(Modia3D.FileMesh)              , kwargs...)
Font(            ; kwargs...) = Par(; _constructor = :(Modia3D.Font)                  , kwargs...)
TextShape(       ; kwargs...) = Par(; _constructor = :(Modia3D.TextShape)             , kwargs...)
Fix(             ; kwargs...) = Par(; _constructor = :(Modia3D.Fix)                   , kwargs...)

MassPropertiesFromShape()  = Par(; _constructor = :(Modia3D.MassPropertiesFromShape))
MassPropertiesFromShapeAndMass(;mass) = Par(; _constructor = :(Modia3D.MassPropertiesFromShapeAndMass), mass = mass)
UniformGravityField(; kwargs...) = Par(; _constructor = :(Modia3D.UniformGravityField), kwargs...)

RefPath(; kwargs...) = Modia3D.ReferencePath(; kwargs...)

ptpJointSpace(; kwargs...) = Modia3D.ptpJointSpace(; kwargs...)
scheduleReferenceMotion(; kwargs...) = Modia3D.scheduleReferenceMotion(; kwargs...)

calculateRobotMovement(args...) = Modia3D.calculateRobotMovement(args...)
getRefPathPosition(args...) = Modia3D.getRefPathPosition(args...)
getRefPathInitPosition(args...) = Modia3D.getRefPathInitPosition(args...)

getVariables(args...) = (args...,)

multibodyResiduals(args...) = Modia3D.multibodyResiduals(args...)

Revolute(; obj1, obj2, axis=3, phi=Var(init=0.0), w=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia3D.Revolute), _path = true, ndof = 1),
    obj1 = Par(value = obj1),
    obj2 = Par(value = obj2),
    axis = Par(value = axis),
    canCollide = Par(value = canCollide),
    phi  = phi,
    w    = w,
    equations = :[
        w   = der(phi)
        qdd = der(w)   # standardized name for the generalized joint accelerations
        variables = getVariables(phi, w, 0.0) # standardized name for the generalized joint position, velocity, force
        ]
)

RevoluteWithFlange(; obj1, obj2, axis=3, phi=Var(init=0.0), w=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia3D.Revolute), _path = true, ndof = 1),
    obj1   = Par(value = obj1),
    obj2   = Par(value = obj2),
    axis   = Par(value = axis),
    canCollide = Par(value = canCollide),
    flange = Flange,
    phi    = phi,
    w      = w,
    equations = :[
        phi = flange.phi
        w   = der(phi)
        qdd = der(w)   # standardized name for the generalized joint accelerations
        variables = getVariables(phi, w, flange.tau) # standardized name for the generalized joint position, velocity, force
        ]
)

Prismatic(; obj1, obj2, axis=1, s=Var(init=0.0), v=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia3D.Prismatic), _path = true, ndof = 1),
    obj1 = Par(value = obj1),
    obj2 = Par(value = obj2),
    axis = Par(value = axis),
    canCollide = Par(value = canCollide),
    s    = s,
    v    = v,
    equations = :[
        v   = der(s)
        qdd = der(v)   # standardized name for the generalized joint accelerations
        variables = getVariables(s, v, 0.0) # standardized name for the generalized joint position, velocity, force
        ]
)

PrismaticWithFlange(; obj1, obj2, axis=1, s=Var(init=0.0), v=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia3D.Prismatic), _path = true, ndof = 1),
    obj1   = Par(value = obj1),
    obj2   = Par(value = obj2),
    axis   = Par(value = axis),
    canCollide = Par(value = canCollide),
    flange = TranslationalFlange,
    s      = s,
    v      = v,
    equations = :[
        s   = flange.s
        v   = der(s)
        qdd = der(v)   # standardized name for the generalized joint accelerations
        variables = getVariables(s, v, flange.f) # standardized name for the generalized joint position, velocity, force
        ]
)

"""
    J = J123(rot123::AbstractVector)

Return joint rot. kinematics matrix `J` for Cardan angles `rot123` (rotation sequence x-y-z).
"""
function J123(rot123::AbstractVector)

    (sbe, cbe) = sincos(rot123[2])
    (sga, cga) = sincos(rot123[3])
    return [     cga     -sga  0.0 ;
             cbe*sga  cbe*cga  0.0 ;
            -sbe*cga  sbe*sga  cbe ] / cbe

end

"""
    J = J132(rot132::AbstractVector)

Return joint rot. kinematics matrix `J` for Cardan angles `rot132` (rotation sequence x-z-y).
"""
function J132(rot132::AbstractVector)

    (sga, cga) = sincos(rot132[2])
    (sbe, cbe) = sincos(rot132[3])
    return [ cbe      0.0  sbe     ;
            -sbe*cga  0.0  cbe*cga ;
             cbe*sga  cga  sbe*sga ] / cga

end

FreeMotion(; obj1, obj2, r=Var(init=zeros(3)), rot=Var(init=zeros(3)), v=Var(init=zeros(3)), w=Var(init=zeros(3))) = Model(; _constructor = Par(value = :(Modia3D.FreeMotion), _path = true, ndof = 6),
    obj1 = Par(value = obj1),
    obj2 = Par(value = obj2),
    r    = r,
    rot  = rot,
    v    = v,
    w    = w,
    equations = :[
        der(r)   = v
        der(rot) = J123(rot) * w
        der(v) = qdd[1:3]
        der(w) = qdd[4:6]
        variables = getVariables(r, rot, v, w)
        ]
)


text:
include("comm_manager.jl")
include("execute_request.jl")

using IJulia.CommManager

function complete_request(socket, msg)
    text = msg.content["text"]
    line = msg.content["line"]
    cursorpos = chr2ind(line, msg.content["cursor_pos"])

    comps, positions = completions(line,cursorpos)
    if sizeof(text) > length(positions)
        comps = [line[(cursorpos-sizeof(text)+1):(cursorpos-length(positions))]*s for s in comps]
    end
    send_ipython(requests, msg_reply(msg, "complete_reply", [
        "status" => "ok",
        "matches" => comps,
        "matched_text" => line[positions],
    ]))
end

function kernel_info_request(socket, msg)
    send_ipython(requests,
                 msg_reply(msg, "kernel_info_reply",
                           ["protocol_version" => [4, 0],
                            "language_version" => [VERSION.major,
                                                   VERSION.minor,
                                                   VERSION.patch],
                            "language" => "julia" ]))
end

function connect_request(socket, msg)
    send_ipython(requests,
                 msg_reply(msg, "connect_reply",
                           ["shell_port" => profile["shell_port"],
                            "iopub_port" => profile["iopub_port"],
                            "stdin_port" => profile["stdin_port"],
                            "hb_port" => profile["hb_port"]]))
end

function shutdown_request(socket, msg)
    send_ipython(requests, msg_reply(msg, "shutdown_reply",
                                     msg.content))
    exit()
end

# TODO: better Julia help integration (issue #13)
docstring(o) = ""
docstring(o::Union(Function,DataType)) = sprint(show, methods(o))

function object_info_request(socket, msg)
    try
        s = parse(msg.content["oname"])
        o = eval(Main, s)
        content = ["oname" => msg.content["oname"],
                   "found" => true,
                   "ismagic" => false,
                   "isalias" => false,
                   "docstring" => docstring(o),
                   "type_name" => string(typeof(o)),
                   "base_class" => string(typeof(o).super),
                   "string_form" => get(msg.content,"detail_level",0) == 0 ?
                   sprint(16384, show, o) : repr(o) ]
        if method_exists(length, (typeof(o),))
            content["length"] = length(o)
        end
        send_ipython(requests, msg_reply(msg, "object_info_reply", content))
    catch e
        @verror_show e catch_backtrace()
        send_ipython(requests,
                     msg_reply(msg, "object_info_reply",
                               ["oname" => msg.content["oname"],
                                "found" => false ]))
    end
end

function history_request(socket, msg)
    # we will just send back empty history for now, pending clarification
    # as requested in ipython/ipython#3806
    send_ipython(requests,
                 msg_reply(msg, "history_reply",
                           ["history" => []]))

end

const handlers = (String=>Function)[
    "execute_request" => execute_request_0x535c5df2,
    "complete_request" => complete_request,
    "kernel_info_request" => kernel_info_request,
    "object_info_request" => object_info_request,
    "connect_request" => connect_request,
    "shutdown_request" => shutdown_request,
    "history_request" => history_request,
    "comm_open" => comm_open,
    "comm_msg" => comm_msg,
    "comm_close" => comm_close
]


text:
function calculate_tgrad(sys::AbstractODESystem)
  isempty(sys.tgrad[]) || return sys.tgrad[]  # use cached tgrad, if possible
  rhs = [detime_dvs(eq.rhs) for eq ∈ equations(sys)]
  iv = sys.iv()
  notime_tgrad = [expand_derivatives(ModelingToolkit.Differential(iv)(r)) for r in rhs]
  tgrad = retime_dvs.(notime_tgrad,(states(sys),),iv)
  sys.tgrad[] = tgrad
  return tgrad
end

function calculate_jacobian(sys::AbstractODESystem)
    isempty(sys.jac[]) || return sys.jac[]  # use cached Jacobian, if possible
    rhs = [eq.rhs for eq ∈ equations(sys)]

    iv = sys.iv()
    dvs = [dv(iv) for dv ∈ states(sys)]

    jac = expand_derivatives.(calculate_jacobian(rhs, dvs))
    sys.jac[] = jac  # cache Jacobian
    return jac
end

struct ODEToExpr
    sys::AbstractODESystem
    states::Vector{Variable}
end
ODEToExpr(sys) = ODEToExpr(sys,states(sys))
function (f::ODEToExpr)(O::Operation)
    if isa(O.op, Variable)
        isequal(O.op, f.sys.iv) && return O.op.name  # independent variable
        O.op ∈ f.states         && return O.op.name  # dependent variables
        isempty(O.args)         && return O.op.name  # 0-ary parameters
        return build_expr(:call, Any[O.op.name; f.(O.args)])
    end
    return build_expr(:call, Any[Symbol(O.op); f.(O.args)])
end
(f::ODEToExpr)(x) = convert(Expr, x)

function generate_tgrad(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); kwargs...)
    tgrad = calculate_tgrad(sys)
    return build_function(tgrad, dvs, ps, sys.iv;
                          conv = ODEToExpr(sys), kwargs...)
end

function generate_jacobian(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); sparse = false, kwargs...)
    jac = calculate_jacobian(sys)
    if sparse
        jac = SparseArrays.sparse(jac)
    end
    return build_function(jac, dvs, ps, sys.iv;
                          conv = ODEToExpr(sys), kwargs...)
end

function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); kwargs...)
    rhss = [deq.rhs for deq ∈ equations(sys)]
    dvs′ = convert.(Variable,dvs)
    ps′ = convert.(Variable,ps)
    return build_function(rhss, dvs′, ps′, sys.iv;
                          conv = ODEToExpr(sys),kwargs...)
end

function calculate_factorized_W(sys::AbstractODESystem, simplify=true)
    isempty(sys.Wfact[]) || return (sys.Wfact[],sys.Wfact_t[])

    jac = calculate_jacobian(sys)
    gam = Variable(:__MTKWgamma)()

    W = - LinearAlgebra.I + gam*jac
    Wfact = lu(W, Val(false), check=false).factors

    if simplify
        Wfact = simplify_constants.(Wfact)
    end

    W_t = - LinearAlgebra.I/gam + jac
    Wfact_t = lu(W_t, Val(false), check=false).factors
    if simplify
        Wfact_t = simplify_constants.(Wfact_t)
    end
    sys.Wfact[] = Wfact
    sys.Wfact_t[] = Wfact_t

    (Wfact,Wfact_t)
end

function generate_factorized_W(sys::AbstractODESystem, vs = states(sys), ps = parameters(sys), simplify=true, expression = Val{true}; kwargs...)
    Wfact,Wfact_t = calculate_factorized_W(sys,simplify)
    siz = size(Wfact)
    constructor = :(x -> begin
                        A = SMatrix{$siz...}(x)
                        StaticArrays.LU(LowerTriangular( SMatrix{$siz...}(UnitLowerTriangular(A)) ), UpperTriangular(A), SVector(ntuple(n->n, max($siz...))))
                    end)

    Wfact_func   = build_function(Wfact  , vs, ps, Variable(:__MTKWgamma), sys.iv;
                                  conv = ODEToExpr(sys), expression = expression, constructor=constructor,kwargs...)
    Wfact_t_func = build_function(Wfact_t, vs, ps, Variable(:__MTKWgamma), sys.iv;
                                  conv = ODEToExpr(sys), expression = expression, constructor=constructor,kwargs...)

    return (Wfact_func, Wfact_t_func)
end

function calculate_massmatrix(sys::AbstractODESystem, simplify=true)
    eqs = equations(sys)
    dvs = states(sys)
    M = zeros(length(eqs),length(eqs))
    for (i,eq) in enumerate(eqs)
        if eq.lhs isa Constant
            @assert eq.lhs.value == 0
        elseif eq.lhs.op isa Differential
            j = findfirst(x->isequal(x.name,var_from_nested_derivative(eq.lhs)[1].name),dvs)
            M[i,j] = 1
        else
            error("Only semi-explicit constant mass matrices are currently supported")
        end
    end
    M = simplify ? simplify_constants.(M) : M
    M == I ? I : M
end

function DiffEqBase.ODEFunction(sys::AbstractODESystem, args...; kwargs...)
    ODEFunction{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false, Wfact = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create an `ODEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and `ps`
are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false, Wfact = false,
                                     sparse = false,
                                     kwargs...) where {iip}
    f_oop,f_iip = generate_function(sys, dvs, ps; expression=Val{false}, kwargs...)

    f(u,p,t) = f_oop(u,p,t)
    f(du,u,p,t) = f_iip(du,u,p,t)

    if tgrad
        tgrad_oop,tgrad_iip = generate_tgrad(sys, dvs, ps; expression=Val{false}, kwargs...)
        _tgrad(u,p,t) = tgrad_oop(u,p,t)
        _tgrad(J,u,p,t) = tgrad_iip(J,u,p,t)
    else
        _tgrad = nothing
    end

    if jac
        jac_oop,jac_iip = generate_jacobian(sys, dvs, ps; sparse = sparse, expression=Val{false}, kwargs...)
        _jac(u,p,t) = jac_oop(u,p,t)
        _jac(J,u,p,t) = jac_iip(J,u,p,t)
    else
        _jac = nothing
    end

    if Wfact
        tmp_Wfact,tmp_Wfact_t = generate_factorized_W(sys, dvs, ps; expression=Val{false}, kwargs...)
        Wfact_oop, Wfact_iip = tmp_Wfact
        Wfact_oop_t, Wfact_iip_t = tmp_Wfact_t
        _Wfact(u,p,dtgamma,t) = Wfact_oop(u,p,dtgamma,t)
        _Wfact(W,u,p,dtgamma,t) = Wfact_iip(W,u,p,dtgamma,t)
        _Wfact_t(u,p,dtgamma,t) = Wfact_oop_t(u,p,dtgamma,t)
        _Wfact_t(W,u,p,dtgamma,t) = Wfact_iip_t(W,u,p,dtgamma,t)
    else
        _Wfact,_Wfact_t = nothing,nothing
    end

    M = calculate_massmatrix(sys)

    ODEFunction{iip}(f,jac=_jac,
                      tgrad = _tgrad,
                      Wfact = _Wfact,
                      Wfact_t = _Wfact_t,
                      mass_matrix = M,
                      syms = Symbol.(states(sys)))
end

function DiffEqBase.ODEProblem(sys::AbstractODESystem, args...; kwargs...)
    ODEProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false, Wfact = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, multithread=false,
                                    kwargs...) where iip
```

Generates an ODEProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false, Wfact = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, multithread=false,
                                    kwargs...) where iip
    f = ODEFunction(sys;tgrad=tgrad,jac=jac,Wfact=Wfact,checkbounds=checkbounds,
                        linenumbers=linenumbers,multithread=multithread,
                        sparse=sparse)
    u0 = varmap_to_vars(u0map,states(sys))
    p = varmap_to_vars(parammap,parameters(sys))
    ODEProblem(f,u0,tspan,p;kwargs...)
end