text:
# using Base.Test
# using MIToS.Utils

print("""

Test General Utils
==================
""")

print("""
deleteitems!
""")
vec = collect("Hola mundo\n")
deleteitems!(vec, Char[' ', '\n'])
@test vec == Char['H', 'o', 'l', 'a', 'm', 'u', 'n', 'd', 'o']

vec = Int[ char for char in "Hola mundo\n" ]
deleteitems!(vec, IntSet(Int[' ', '\n']) )
@test vec == Int['H', 'o', 'l', 'a', 'm', 'u', 'n', 'd', 'o']

print("""
get_n_words
""")
line = "#=GF AC PF00571"
line_n = "#=GF AC PF00571\n"
@test get_n_words(line, 1) == ASCIIString[line]
@test get_n_words(line, 2) == ASCIIString["#=GF", "AC PF00571"]
@test get_n_words(line, 3) == ASCIIString["#=GF", "AC", "PF00571"]
@test get_n_words(line, 4) == ASCIIString["#=GF", "AC", "PF00571"]
for i in 1:4
  @test get_n_words(line, i) == get_n_words(line_n, i)
end

@test get_n_words("\n",1) == ASCIIString[]
@test get_n_words("#", 1) == ASCIIString["#"]

let str = "#=GF CC   (Römling U.  and Galperin M.Y. “Bacterial cellulose\n"

  @test get_n_words(str, 3) == UTF8String["#=GF", "CC", "(Römling U.  and Galperin M.Y. “Bacterial cellulose"]
end

let str = "#=GF CC   not present in all SecA2–SecY2 systems. This family of Asp5 is\n"

  @test get_n_words(str, 3) == UTF8String["#=GF", "CC", "not present in all SecA2–SecY2 systems. This family of Asp5 is"]
end

print("""
hascoordinates
""")
@test hascoordinates("O83071/192-246")
@test !hascoordinates("O83071")

print("""
select_element
""")
@test select_element([1]) == 1
# print("Test warn: ")
# @test select_element([1,2]) == 1
@test_throws ErrorException select_element([])

print("""

matrix2list
-----------
""")

let mat = [ 1 2 3
            4 5 6
            7 8 9 ]

  @test matrix2list(mat) == [2, 3, 6]
  @test matrix2list(mat, diagonal=true) == [1, 2, 3, 5, 6, 9]
  @test matrix2list(mat, part="lower") == [4, 7, 8]
  @test matrix2list(mat, part="lower", diagonal=true) == [1, 4, 7, 5, 8, 9]
end

print("""

list2matrix
-----------
""")

let mat = [ 1 2 3
            2 5 6
            3 6 9 ]

  @test triu(list2matrix([2, 3, 6], 3), 1) == triu(mat, 1)
  @test list2matrix([1, 2, 3, 5, 6, 9], 3, diagonal=true) == mat
end

print("""

Test findobjects
================
""")

import Base: ==, hash

using AutoHashEquals

@auto_hash_equals immutable Dummy
  string::ASCIIString
  int::Int
end

print("""
findobjects, isobject and AbstractTest
""")

@test findobjects([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], Not(Is(:int,4))) == [1, 2]
@test findobjects([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], Is(:int,4)) == [3]

@test findobjects([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], Is(:int, x -> x > 2 )) == [2, 3]
@test findobjects([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], Not(Is(:int, x -> x > 2 ))) == [1]

@test findobjects([ Dummy("abc", 2), Dummy("abcd", 2), Dummy("bc", 2) ], Is(:string, r"^ab")) == [1, 2]
@test findobjects([ Dummy("abc", 2), Dummy("abcd", 2), Dummy("bc", 2) ], Not(Is(:string, r"^ab"))) == [3]
@test findobjects([ Dummy("abc", 2), Dummy("abcd", 2), Dummy("bc", 2) ], Is(:string, r"^nothere$")) == []

@test findobjects([ Dummy("H", 2), Dummy("C", 2), Dummy("O", 2) ], Not(Is(:string, "H"))) == [2, 3]

@test findobjects([ Dummy("H", 2), Dummy("C", 2), Dummy("O", 2) ], In(:string, ["C", "O"])) == [2, 3]
@test findobjects([ Dummy("H", 2), Dummy("C", 2), Dummy("O", 2) ], Not(In(:string, ["C", "O"]))) == [1]

print("""
collectobjects & isobject
""")

@test collectobjects([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], Is(:int,4)) == [Dummy("c", 4)]
@test collectobjects([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], Is(:int, x -> x > 2 )) == [Dummy("b", 3), Dummy("c", 4)]
@test collectobjects([ Dummy("abc", 2), Dummy("abcd", 2), Dummy("bc", 2) ], Not(Is(:string, r"^ab"))) == [Dummy("bc", 2)]
@test collectobjects([ Dummy("abc", 2), Dummy("abcd", 2), Dummy("bc", 2) ], Is(:string, r"^nothere$")) == []

print("""
collectcaptures & capture
""")

@test map(isnull, collectcaptures([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], :string, Is( :int, x -> x > 2 )))  == [true, false, false]
@test map(isnull, collectcaptures([ "NotDummy", 42, Dummy("c", 4) ], :string, Is( :int, x -> x > 2 )))   == [true, true, false]

print("""
collectcaptures & guess_type
""")

@test isa(collectcaptures([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], :string, Is( :int, 3 )), Array{Nullable{ASCIIString},1})
@test isa(collectcaptures([ Dummy("a", 2), Dummy("b", 3), Dummy("c", 4) ], :int, Is( :int, 3 )), Array{Nullable{Int},1})


text:
using GLAbstraction, Makie, GeometryTypes
using GLAbstraction: orthographicprojection, translationmatrix
using StaticArrays

const AbstractCamera = Scene
@default function camera2d(scene, kw_args)
    area = SimpleRectangle(0, 0, 500, 500)
    projection = eye(Mat4f0)
    view = eye(Mat4f0)

    translationspeed = to_float(1)
    eyeposition = Vec3f0(3)
    lookat = Vec3f0(0)
    upvector = Vec3f0((0, 0, 1))
    near = to_float(10_000)
    far = to_float(-10_000)
end

wscale(screenrect, viewrect) = widths(viewrect) ./ widths(screenrect)

function update_cam!(cam::AbstractCamera, area)
    x, y = minimum(area)
    w, h = widths(area) ./ 2f0
    cam[:area] = FRect(area)
    cam[:projection] = orthographicprojection(-w, w, -h, h, -10_000f0, 10_000f0)
    cam[:view] = translationmatrix(Vec3f0(-x - w, -y - h, 0))
    return
end

function add_mousebuttons(scene::Scene)
    nw = GLWindow.nativewindow(scene[:screen])
    scene[:mousebuttons] = Set(Mouse.Button[])
    GLFW.SetMouseButtonCallback(nw, (native_window, button, action, mods) -> begin
        set = to_value(scene, :mousebuttons)
        button_enum = Mouse.Button(button)
        if button != GLFW.KEY_UNKNOWN
            if action == GLFW.PRESS
                push!(set, button_enum)
            elseif action == GLFW.RELEASE
                delete!(set, button_enum)
            elseif action == GLFW.REPEAT
                # nothing needs to be done, besides returning the same set of keys
            else
                error("Unrecognized enum value for GLFW button press action: $action")
            end
        end
        scene[:mousebuttons] = set # trigger setfield event!
        return
    end)
    return
end

function camspace(cam_area, screen_area, point)
    point = point .* wscale(screen_area, cam_area)
    point .+ minimum(cam_area)
end

function selection_rect(
        scene, cam,
        key = Mouse.left,
        button = Set([Keyboard.left_control, Keyboard.space])
    )
    rect = RefValue(FRect(0, 0, 0, 0))
    lw = 2f0
    rect_vis = lines(
        scene,
        rect[],
        linestyle = :dot,
        thickness = 1f0,
        color = (:black, 0.4),
        drawover = true,
        camera = :pixel
    )
    waspressed = RefValue(false)
    dragged_rect = lift_node(scene, :mousedrag) do drag
        if ispressed(scene, key) && ispressed(scene, button)
            screen_area = to_value(scene, :window_area)
            cam_area = to_value(cam, :area)
            mp = to_value(scene, :mouseposition)
            mp = camspace(cam_area, screen_area, mp)

            if drag == Mouse.down
                waspressed[] = true
                rect_vis[:visible] = true # start displaying
                rect[] = FRect(mp, 0, 0)
                rect_vis[:positions] = rect[]
            elseif drag == Mouse.pressed
                mini = minimum(rect[])
                rect[] = FRect(mini, mp - mini)
                # mini, maxi = min(mini, mp), max(mini, mp)
                rect_vis[:positions] = rect[]
            end
        else
            if drag == Mouse.up && waspressed[]
                waspressed[] = false
                update_cam!(cam, rect[])
            end
            rect[] = FRect(0, 0, 0, 0)
            rect_vis[:positions] = rect[]
            # always hide if not the right key is pressed
            rect_vis[:visible] = false # hide
        end
        return rect[]
    end
    rect_vis, dragged_rect
end

function add_mousedrag(scene::Scene)
    indrag = RefValue(false)
    tracked_mousebutton = RefValue(Mouse.left)
    drag = RefValue(Mouse.notpressed)
    scene[:mousedrag] = lift_node(scene, :mouseposition, :mousebuttons) do mp, buttons
        # only track if still the same button is pressed
        if length(buttons) == 1 &&
                (!indrag[] || tracked_mousebutton[] == first(buttons))

            if !indrag[]
                tracked_mousebutton[] = first(buttons); indrag[] = true
                drag[] = Mouse.down # just started, so dragging is still false
                return drag[]
            else
                drag[] = Mouse.pressed
                return drag[]
            end
        end
        # already on notpressed, no need for update
        if drag[] != Mouse.notpressed
            drag[] = indrag[] ? Mouse.up : Mouse.notpressed
        end
        indrag[] = false
        return drag[]
    end
end

function add_pan(cam::AbstractCamera, scene)
    startpos = RefValue(Vec(0.0, 0.0))
    lift_node(scene, :mouseposition, :mousedrag) do mp, dragging
        if ispressed(scene, Mouse.middle)
            screen_area = to_value(scene, :window_area)
            cam_area = to_value(cam, :area)
            if dragging == Mouse.down
                startpos[] = mp
            elseif dragging == Mouse.pressed && ispressed(scene, Mouse.middle)
                diff = startpos[] .- mp
                startpos[] = mp
                diff = diff .* wscale(screen_area, cam_area)
                update_cam!(cam, FRect(minimum(cam_area) .+ diff, widths(cam_area)))
            end
        end
        return
    end
end

function add_zoom(cam::AbstractCamera, scene)
    lift_node(scene[:scroll]) do x
        zoom = Float32(x[2])
        if zoom != 0
            a = to_value(cam, :area)
            z = 1 + (zoom * 0.10)
            mp = Vec2f0(to_value(scene, :mouseposition))
            mp = (mp .* wscale(to_value(scene, :window_area), a)) + minimum(a)
            p1, p2 = minimum(a), maximum(a)
            p1, p2 = p1 - mp, p2 - mp # translate to mouse position
            p1, p2 = z * p1, z * p2
            p1, p2 = p1 + mp, p2 + mp
            update_cam!(cam, FRect(p1, p2 - p1))
            z
        end
        return
    end
end

function reset!(cam::AbstractCamera, boundingbox, preserveratio = true)
    w1 = widths(boundingbox)
    if preserveratio
        w2 = widths(cam[Screen][Area])
        ratio = w2 ./ w1
        w1 = if ratio[1] > ratio[2]
            s = w2[1] ./ w2[2]
            Vec2f0(s * w1[2], w1[2])
        else
            s = w2[2] ./ w2[1]
            Vec2f0(w1[1], s * w1[1])
        end
    end
    p = minimum(w1) .* 0.001 # 2mm padding
    update_cam!(cam, FRect(-p, -p, w1 .+ 2p))
    return
end

lerp{T}(a::T, b::T, val::AbstractFloat) = (a .+ (val * (b .- a)))

function add_restriction!(cam::AbstractCamera, window, rarea::SimpleRectangle, minwidths::Vec)
    area_ref = Base.RefValue(cam[Area])
    restrict_action = paused_action(1.0) do t
        o = lerp(origin(area_ref[]), origin(cam[Area]), t)
        wh = lerp(widths(area_ref[]), widths(cam[Area]), t)
        update_cam!(cam, FRect(o, wh))
    end
    on(window, Mouse.Drag) do drag
        if drag == Mouse.up && !isplaying(restrict_action)
            area = cam[Area]
            o = origin(area)
            maxi = maximum(area)
            newo = max.(o, origin(rarea))
            newmax = min.(maxi, maximum(rarea))
            maxi = maxi - newmax
            newo = newo - maxi
            newwh = newmax - newo
            scale = 1f0
            for (w1, w2) in zip(minwidths, newwh)
                stmp = w1 > w2 ? w1 / w2 : 1f0
                scale = max(scale, stmp)
            end
            newwh = newwh * scale
            area_ref[] = FRect(newo, newwh)
            if area_ref[] != cam[Area]
                play!(restrict_action)
            end
        end
        return
    end
    restrict_action
end

ext: jl


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

const LOG_STD_MAX = 2
const LOG_STD_MIN = -20

"""
Create multilayer perceptron with given parameters.
"""
function mlp(sizes::Vector{Int}, activation::Function, output_activation::Function=identity)
    layers = []
    for j = 1:length(sizes) - 1
        act = j < length(sizes) - 1 ? activation : output_activation
        push!(layers, Dense(sizes[j], sizes[j + 1], act))
    end
    return Chain(layers...) |> dev
end

"""
Gaussian MLP policy.
"""
mutable struct SquashedGaussianMLPActor
    net::Chain                              # primary network
    mu_layer::Dense                         # mean layer
    log_std_layer::Dense                    # log standard deviation layer
    act_mins::AbstractVector{Float32}       # minimum values of actions
    act_maxs::AbstractVector{Float32}       # maximum values of actions
    rng::AbstractRNG                        # internal RNG
    rng_gpu::Union{AbstractRNG, Nothing}    # internal RNG for GPU calculations
    linearized::Bool                        # linearized squashing (allows post-analysis)
end

function SquashedGaussianMLPActor(
	obs_dim::Int,
	act_dim::Int,
	hidden_sizes::Vector{Int},
	activation::Function,
	act_mins::Vector{Float64},
	act_maxs::Vector{Float64},
	rng::AbstractRNG,
    linearized::Bool
)
    net = mlp(vcat(obs_dim, hidden_sizes), activation, activation) |> dev
    mu_layer = Dense(hidden_sizes[end], act_dim) |> dev
    log_std_layer = Dense(hidden_sizes[end], act_dim) |> dev
    act_mins, act_maxs = dev(Float32.(act_mins)), dev(Float32.(act_maxs))
    if WITH_GPU[]
        rng_gpu = CURAND.RNG()
        seed = rand(rng, UInt32)
        Random.seed!(rng_gpu, seed)
    else
        rng_gpu = nothing
    end
    return SquashedGaussianMLPActor(net, mu_layer, log_std_layer, act_mins, act_maxs, rng, rng_gpu, linearized)
end

"""
Calculate log pdf of normal distribution.
"""
function normal_logpdf(μ, σ, x)
    lz = sum(-0.5f0 .* ((x .- μ) ./ σ).^2; dims=1)
    lden = (size(μ, 1) * Float32(log(2π)) / 2) .+ sum(log.(σ); dims=1)
    lpdf = lz .- lden
    lpdf = dropdims(lpdf; dims=1)
    return lpdf
end

"""
Calculate random normal deviates on GPU (in-place).
"""
function normal_deviates(rng::AbstractRNG, sz::NTuple{N, Int64}) where N
    z = CuArray{Float32}(undef, sz)
    randn!(rng, z)
    return z
end
Zygote.@nograd normal_deviates

"""
Retrieve action (and optional log probability) from policy.
"""
function (pi::SquashedGaussianMLPActor)(obs, deterministic::Bool=false, with_logprob::Bool=true)
    net_out = pi.net(obs)
    mu = pi.mu_layer(net_out)
    log_std = pi.log_std_layer(net_out)
    log_std = clamp.(log_std, LOG_STD_MIN, LOG_STD_MAX)
    std = exp.(log_std)

    # Pre-squash distribution and sample
    if deterministic
        pi_action = mu
    else
        if obs isa CuArray
            z = normal_deviates(pi.rng_gpu, size(mu))
        else
            z = randn(pi.rng, Float32, size(mu))
        end
        pi_action = mu .+ std .* Zygote.dropgrad(z)
    end

    if with_logprob
        logp_pi = normal_logpdf(mu, std, pi_action)
        logp_pi = logp_pi .- dropdims(sum((2.0f0 .* (log(2.0f0) .- pi_action .- softplus.(-2.0f0 .* pi_action))); dims=1); dims=1)
    else
        logp_pi = NaN
    end

    if pi.linearized # probabilities still correspond to non-linear squashing
        pi_action = clamp.(pi_action, pi.act_mins, pi.act_maxs)
    else
        pi_action = tanh.(pi_action)
        pi_action = @. pi.act_mins + (pi.act_maxs - pi.act_mins) * (pi_action / 2 + 0.5f0)
    end

    return pi_action, logp_pi
end

"""
Q-value function.
"""
mutable struct MLPQFunction
    q::Chain
end

function MLPQFunction(obs_dim::Int, act_dim::Int, hidden_sizes::Vector{Int}, activation::Function)
    q = mlp(vcat(obs_dim + act_dim, hidden_sizes, 1), activation)
    return MLPQFunction(q)
end

"""
Determine Q-value of observation and action.
"""
function (qf::MLPQFunction)(obs, act)
    q = qf.q(cat(obs, act; dims=1))
    q = dropdims(q; dims=1)
    return q
end

"""
Actor-critic agent.
"""
mutable struct MLPActorCritic <: GlobalResult
    pi::SquashedGaussianMLPActor
    qs::Vector{MLPQFunction}
end

function MLPActorCritic(
	obs_dim::Int,
	act_dim::Int,
	act_mins::Vector{Float64},
	act_maxs::Vector{Float64},
	hidden_sizes::Vector{Int}=[100,100,100],
	num_q::Int=2,
	activation::Function=SoftActorCritic.relu,
	rng::AbstractRNG=Random.GLOBAL_RNG,
    linearized::Bool=false
)
    pi = SquashedGaussianMLPActor(obs_dim, act_dim, hidden_sizes, activation, act_mins, act_maxs, rng, linearized)
    qs = [MLPQFunction(obs_dim, act_dim, hidden_sizes, activation) for _ in 1:num_q]
    return MLPActorCritic(pi, qs)
end

"""
Retrieve action from policy.
"""
function (ac::MLPActorCritic)(obs, deterministic::Bool=true)
    a, _ = ac.pi(obs, deterministic, false)
    return a
end

"""
Define native softplus function to avoid CUDA bugs. #TODO: removable?
"""
softplus(x) = x > 0 ? x + log(oftype(x, 1) + exp(-x)) : log(oftype(x, 1) + exp(x))

"""
Define native relu function to avoid Flux bugs.
"""
relu(x) = max(x, oftype(x, 0))


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

mutable struct SecurityRulePropertiesFormat2 <: SwaggerModel
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    sourcePortRange::Any # spec type: Union{ Nothing, String } # spec name: sourcePortRange
    destinationPortRange::Any # spec type: Union{ Nothing, String } # spec name: destinationPortRange
    sourceAddressPrefix::Any # spec type: Union{ Nothing, String } # spec name: sourceAddressPrefix
    sourceAddressPrefixes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: sourceAddressPrefixes
    sourceApplicationSecurityGroups::Any # spec type: Union{ Nothing, Vector{ApplicationSecurityGroup} } # spec name: sourceApplicationSecurityGroups
    destinationAddressPrefix::Any # spec type: Union{ Nothing, String } # spec name: destinationAddressPrefix
    destinationAddressPrefixes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: destinationAddressPrefixes
    destinationApplicationSecurityGroups::Any # spec type: Union{ Nothing, Vector{ApplicationSecurityGroup} } # spec name: destinationApplicationSecurityGroups
    sourcePortRanges::Any # spec type: Union{ Nothing, Vector{String} } # spec name: sourcePortRanges
    destinationPortRanges::Any # spec type: Union{ Nothing, Vector{String} } # spec name: destinationPortRanges
    access::Any # spec type: Union{ Nothing, SecurityRuleAccess2 } # spec name: access
    priority::Any # spec type: Union{ Nothing, Int32 } # spec name: priority
    direction::Any # spec type: Union{ Nothing, SecurityRuleDirection2 } # spec name: direction
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function SecurityRulePropertiesFormat2(;description=nothing, protocol=nothing, sourcePortRange=nothing, destinationPortRange=nothing, sourceAddressPrefix=nothing, sourceAddressPrefixes=nothing, sourceApplicationSecurityGroups=nothing, destinationAddressPrefix=nothing, destinationAddressPrefixes=nothing, destinationApplicationSecurityGroups=nothing, sourcePortRanges=nothing, destinationPortRanges=nothing, access=nothing, priority=nothing, direction=nothing, provisioningState=nothing)
        o = new()
        validate_property(SecurityRulePropertiesFormat2, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(SecurityRulePropertiesFormat2, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(SecurityRulePropertiesFormat2, Symbol("sourcePortRange"), sourcePortRange)
        setfield!(o, Symbol("sourcePortRange"), sourcePortRange)
        validate_property(SecurityRulePropertiesFormat2, Symbol("destinationPortRange"), destinationPortRange)
        setfield!(o, Symbol("destinationPortRange"), destinationPortRange)
        validate_property(SecurityRulePropertiesFormat2, Symbol("sourceAddressPrefix"), sourceAddressPrefix)
        setfield!(o, Symbol("sourceAddressPrefix"), sourceAddressPrefix)
        validate_property(SecurityRulePropertiesFormat2, Symbol("sourceAddressPrefixes"), sourceAddressPrefixes)
        setfield!(o, Symbol("sourceAddressPrefixes"), sourceAddressPrefixes)
        validate_property(SecurityRulePropertiesFormat2, Symbol("sourceApplicationSecurityGroups"), sourceApplicationSecurityGroups)
        setfield!(o, Symbol("sourceApplicationSecurityGroups"), sourceApplicationSecurityGroups)
        validate_property(SecurityRulePropertiesFormat2, Symbol("destinationAddressPrefix"), destinationAddressPrefix)
        setfield!(o, Symbol("destinationAddressPrefix"), destinationAddressPrefix)
        validate_property(SecurityRulePropertiesFormat2, Symbol("destinationAddressPrefixes"), destinationAddressPrefixes)
        setfield!(o, Symbol("destinationAddressPrefixes"), destinationAddressPrefixes)
        validate_property(SecurityRulePropertiesFormat2, Symbol("destinationApplicationSecurityGroups"), destinationApplicationSecurityGroups)
        setfield!(o, Symbol("destinationApplicationSecurityGroups"), destinationApplicationSecurityGroups)
        validate_property(SecurityRulePropertiesFormat2, Symbol("sourcePortRanges"), sourcePortRanges)
        setfield!(o, Symbol("sourcePortRanges"), sourcePortRanges)
        validate_property(SecurityRulePropertiesFormat2, Symbol("destinationPortRanges"), destinationPortRanges)
        setfield!(o, Symbol("destinationPortRanges"), destinationPortRanges)
        validate_property(SecurityRulePropertiesFormat2, Symbol("access"), access)
        setfield!(o, Symbol("access"), access)
        validate_property(SecurityRulePropertiesFormat2, Symbol("priority"), priority)
        setfield!(o, Symbol("priority"), priority)
        validate_property(SecurityRulePropertiesFormat2, Symbol("direction"), direction)
        setfield!(o, Symbol("direction"), direction)
        validate_property(SecurityRulePropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type SecurityRulePropertiesFormat2

const _property_map_SecurityRulePropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("description")=>Symbol("description"), Symbol("protocol")=>Symbol("protocol"), Symbol("sourcePortRange")=>Symbol("sourcePortRange"), Symbol("destinationPortRange")=>Symbol("destinationPortRange"), Symbol("sourceAddressPrefix")=>Symbol("sourceAddressPrefix"), Symbol("sourceAddressPrefixes")=>Symbol("sourceAddressPrefixes"), Symbol("sourceApplicationSecurityGroups")=>Symbol("sourceApplicationSecurityGroups"), Symbol("destinationAddressPrefix")=>Symbol("destinationAddressPrefix"), Symbol("destinationAddressPrefixes")=>Symbol("destinationAddressPrefixes"), Symbol("destinationApplicationSecurityGroups")=>Symbol("destinationApplicationSecurityGroups"), Symbol("sourcePortRanges")=>Symbol("sourcePortRanges"), Symbol("destinationPortRanges")=>Symbol("destinationPortRanges"), Symbol("access")=>Symbol("access"), Symbol("priority")=>Symbol("priority"), Symbol("direction")=>Symbol("direction"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_SecurityRulePropertiesFormat2 = Dict{Symbol,String}(Symbol("description")=>"String", Symbol("protocol")=>"String", Symbol("sourcePortRange")=>"String", Symbol("destinationPortRange")=>"String", Symbol("sourceAddressPrefix")=>"String", Symbol("sourceAddressPrefixes")=>"Vector{String}", Symbol("sourceApplicationSecurityGroups")=>"Vector{ApplicationSecurityGroup}", Symbol("destinationAddressPrefix")=>"String", Symbol("destinationAddressPrefixes")=>"Vector{String}", Symbol("destinationApplicationSecurityGroups")=>"Vector{ApplicationSecurityGroup}", Symbol("sourcePortRanges")=>"Vector{String}", Symbol("destinationPortRanges")=>"Vector{String}", Symbol("access")=>"SecurityRuleAccess2", Symbol("priority")=>"Int32", Symbol("direction")=>"SecurityRuleDirection2", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ SecurityRulePropertiesFormat2 }) = collect(keys(_property_map_SecurityRulePropertiesFormat2))
Swagger.property_type(::Type{ SecurityRulePropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SecurityRulePropertiesFormat2[name]))}
Swagger.field_name(::Type{ SecurityRulePropertiesFormat2 }, property_name::Symbol) =  _property_map_SecurityRulePropertiesFormat2[property_name]

const _allowed_SecurityRulePropertiesFormat2_protocol = ["Tcp", "Udp", "Icmp", "Esp", "*", "Ah"]

function check_required(o::SecurityRulePropertiesFormat2)
    (getproperty(o, Symbol("access")) === nothing) && (return false)
    (getproperty(o, Symbol("direction")) === nothing) && (return false)
    true
end

function validate_property(::Type{ SecurityRulePropertiesFormat2 }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "SecurityRulePropertiesFormat2", :enum, val, _allowed_SecurityRulePropertiesFormat2_protocol)
    end
end


text:
# # Plot Comparison Between Model Runs

# This plots zonal and time mean u anomalies from a control run (assumed to be the first in names below). Used for sensitivity analysis of parameter sweeps.

# Load modules
using Plots
using NCDatasets
using Statistics: mean

# Select the variable to plot
var_name = "u";

# Specify needed directories and filenames: these can be replaced automatically by a bash run script or can be modified manually
CLIMA_ANALYSIS = "to_be_specified";
CLIMA_NETCDF = "to_be_specified";
RUNNAME = "to_be_specified";

# Get all netcdf file names with a specified string in the CLIMA_NETCDF directory
fnames = filter(x -> occursin(".nc", x), readdir( CLIMA_NETCDF ) );

# Aggregate data in output into a multi-file dataset
nexp = size(fnames);

# Extract and postprocess data
function get_zonal_mean(file_name, var_name, t_spinup);
    ds = NCDataset(file_name, "r");
    lon = ds["long"][:];
    lat = ds["lat"][:];
    z = ds["level"][:] / 1e3; # height in kilometers
    time = ds["time"][:];  # time
    var = ds[var_name][:];
    var=var[t_spinup:end, :, :, :]; # cut out the spinup;
    data_time_zonal_mean = mean( mean(var, dims=1) , dims=4)[1,:,:,1]; # create statistics
    return lat, z, data_time_zonal_mean
    none
end

# ## Plotting

# Control
lat, z, zm_c= get_zonal_mean( "$CLIMA_NETCDF/"fnames[1], var_name, 1);
p_ctrl = contourf( lat, z, zm_c , title="$var_name (ctrl)", xlabel="lat (deg N)", ylabel="z (km)");

# Make anomaly plots and save them in an array
plot_array = Any[p_ctrl];
for i in 2:nexp[1]
  lat, z, zm = get_zonal_mean( "$CLIMA_NETCDF/"fnames[i], var_name, 1);
  one_plot = contourf( lat, z, zm-zm_c, title="$var_name (exp$i-ctrl)", xlabel="lat (deg N)", ylabel="z (km)");
  push!(plot_array,one_plot); # make a plot and add it to the plot_array
end

# Plot and save the array of plots
fig=plot(plot_array...);
savefig(fig, string("$CLIMA_ANALYSIS/plot_$var_name","_ctrl_anomaly.pdf"));


text:
## If Rois are given as an array of indices
function FluorescentSerie(img::AxisArray,rois::Array{Array{Int64,1},1},summaryFunc::Function=sum)
    nt = nimages(img)
    ax = timeaxis(img)
    results = zeros(nt,length(rois))
    for i in eachindex(rois)
        for j in 1:nt
            results[j,i] = summaryFunc(img[ax(j)][rois[i]])
        end
    end
    AxisArray(results,ax,Axis{:ROI}(1:length(rois)))
end

# Constructing from a ROI image and the raw data
function FluorescentSerie(rawImage::AxisArray,roiIm::AbstractArray{Int64},summaryFunc::Function=sum)
    size_spatial(rawImage) != size(roiIm) ? error("ROI image has a different size than the data") :
    rois = Array{Array{Int64,1},1}(maximum(roiIm))
    for i in 1:maximum(roiIm)
        rois[i] = LinearIndices(findall(roiIm.==i))
    end
    FluorescentSerie(rawImage,rois,summaryFunc)
end

# Constructing from an Array and  associated metadata
function FluorescentSerie(img::AbstractArray,metadata::Dict,rois::Array{Array{CartesianIndex{3},1},1},summaryFunc::Function=sum)

    framesPerTrial = metadata["framesPerTrial"]
    nMax = maximum(framesPerTrial)
    results = zeros(Float64,(nMax,length(rois),length(framesPerTrial)))
    ax = range(0,step=metadata["samplingTime"],length=nMax)
    for i in eachindex(rois)
        startPoint = 0
        roiIm = img[rois[i],:]
        for t in 1:length(framesPerTrial)
            nt = framesPerTrial[t]
            for j in 1:nt
                results[j,i,t] = summaryFunc(roiIm[:,startPoint+j])
            end
            startPoint+=nt
        end
    end
    AxisArray(results,Axis{:time}(ax),Axis{:ROI}(1:length(rois)),Axis{:trial}(1:length(framesPerTrial)))
end

function FluorescentSerie(rawImage::AbstractArray,metadata::Dict,roiIm::AbstractArray{Int64},summaryFunc::Function=sum)
    size(rawImage)[1:3] != size(roiIm)[1:3] ? error("ROI image has a different size than the data") :
    rois = Array{Array{CartesianIndex{3},1},1}(undef,maximum(roiIm))
    for i in 1:maximum(roiIm)
        rois[i] = findall(roiIm[:,:,:,1].==i)
    end
    FluorescentSerie(rawImage,metadata,rois,summaryFunc)
end


text:
@time @testset "Foos" begin

  @testset "EnglishFoo" begin

    foo = EnglishFoo()

    @test phrase(foo) == "Hi everybody! I am "

  end

  @testset "GermanFoo" begin

    foo = GermanFoo()

    @test phrase(foo) == "Guten Tag! Ich bin "

  end

end


text:
function build_nabla2_z(LF)

  D2jl = LF.LFz.D2jl

  Nx = LF.Nx
  Ny = LF.Ny
  Nz = LF.Nz

  Npoints = Nx*Ny*Nz

  # initialize with
  rowval = zeros(Int64,Nz*Npoints)
  colptr = zeros(Int64,Npoints+1)
  colptr[1] = 1
  nzval = zeros(Float64,Nz*Npoints)

  idx = 0
  idxCol = 1

  #print("\n")
  for ic = 1:Nx*Ny
    for colLoc = 1:Nz
      for rowLoc = 1:Nz
        idx = idx + 1
        rowGbl = rowLoc + (ic-1)*Nz
        rowval[idx] = rowGbl
        nzval[idx] = D2jl[rowLoc,colLoc]
        #println("colGbl = ", colGbl, " rowLoc = ", rowLoc, " colLoc = ", colLoc)
      end
      idxCol = idxCol + 1
      colptr[idxCol] = colptr[idxCol-1] + Nz
    end
  end

  return SparseMatrixCSC( Npoints, Npoints, colptr, rowval, nzval )

end


text:
function add(a, b)
	a + b
end


text:
const nnp_status = UInt32
const nnp_status_success = (UInt32)(0)
const nnp_status_invalid_batch_size = (UInt32)(2)
const nnp_status_invalid_channels = (UInt32)(3)
const nnp_status_invalid_input_channels = (UInt32)(4)
const nnp_status_invalid_output_channels = (UInt32)(5)
const nnp_status_invalid_input_size = (UInt32)(10)
const nnp_status_invalid_input_stride = (UInt32)(11)
const nnp_status_invalid_input_padding = (UInt32)(12)
const nnp_status_invalid_kernel_size = (UInt32)(13)
const nnp_status_invalid_pooling_size = (UInt32)(14)
const nnp_status_invalid_pooling_stride = (UInt32)(15)
const nnp_status_invalid_algorithm = (UInt32)(16)
const nnp_status_invalid_transform_strategy = (UInt32)(17)
const nnp_status_invalid_output_subsampling = (UInt32)(13)
const nnp_status_invalid_activation = (UInt32)(14)
const nnp_status_invalid_activation_parameters = (UInt32)(15)
const nnp_status_unsupported_input_size = (UInt32)(20)
const nnp_status_unsupported_input_stride = (UInt32)(21)
const nnp_status_unsupported_input_padding = (UInt32)(22)
const nnp_status_unsupported_kernel_size = (UInt32)(23)
const nnp_status_unsupported_pooling_size = (UInt32)(24)
const nnp_status_unsupported_pooling_stride = (UInt32)(25)
const nnp_status_unsupported_algorithm = (UInt32)(26)
const nnp_status_unsupported_transform_strategy = (UInt32)(57)
const nnp_status_unsupported_activation = (UInt32)(28)
const nnp_status_unsupported_activation_parameters = (UInt32)(29)
const nnp_status_uninitialized = (UInt32)(50)
const nnp_status_unsupported_hardware = (UInt32)(51)
const nnp_status_out_of_memory = (UInt32)(52)
const nnp_status_insufficient_buffer = (UInt32)(53)
const nnp_status_misaligned_buffer = (UInt32)(54)

const nnp_activation = UInt32
const nnp_activation_identity = (UInt32)(0)
const nnp_activation_relu = (UInt32)(1)

const nnp_convolution_algorithm = UInt32
const nnp_convolution_algorithm_auto = (UInt32)(0)
const nnp_convolution_algorithm_ft8x8 = (UInt32)(1)
const nnp_convolution_algorithm_ft16x16 = (UInt32)(2)
const nnp_convolution_algorithm_wt8x8 = (UInt32)(3)
const nnp_convolution_algorithm_implicit_gemm = (UInt32)(4)
const nnp_convolution_algorithm_direct = (UInt32)(5)
const nnp_convolution_algorithm_wt8x8_fp16 = (UInt32)(6)

const nnp_convolution_transform_strategy = UInt32
const nnp_convolution_transform_strategy_compute = (UInt32)(1)
const nnp_convolution_transform_strategy_precompute = (UInt32)(2)
const nnp_convolution_transform_strategy_reuse = (UInt32)(3)

const pthreadpool_t = Ptr{Nothing}

mutable struct nnp_size
    width::Csize_t
    height::Csize_t
    nnp_size() = new(Csize_t(0), Csize_t(0))
    nnp_size(w, h) = new(Csize_t(w), Csize_t(h))
end

Base.unsafe_convert(::Type{Ptr{nnp_size}}, a::nnp_size) = Ptr{a}

mutable struct nnp_padding
    top::Csize_t
    right::Csize_t
    bottom::Csize_t
    left::Csize_t
    nnp_padding() = new(Csize_t(0), Csize_t(0), Csize_t(0), Csize_t(0))
    nnp_padding(val) = new(Csize_t(val), Csize_t(val), Csize_t(val), Csize_t(val))
    nnp_padding(t, r, b, l) = new(Csize_t(t), Csize_t(r), Csize_t(b), Csize_t(l))
end

Base.unsafe_convert(::Type{Ptr{nnp_padding}}, a::nnp_padding) = Ptr{a}

mutable struct nnp_profile
    total::Cdouble
    input_transform::Cdouble
    kernel_transform::Cdouble
    output_transform::Cdouble
    block_multiplication::Cdouble
    nnp_profile() = new(Cdouble(0.0), Cdouble(0.0), Cdouble(0.0), Cdouble(0.0), Cdouble(0.0))
    nnp_profile(t, it, kt, ot, bm) = new(Cdouble(t), Cdouble(it), Cdouble(kt), Cdouble(ot), Cdouble(bm))
end

Base.unsafe_convert(::Type{Ptr{nnp_profile}}, a::nnp_profile) = Ptr{a}


text:
using Trapz

export RHS_Spiral, get_g_ft, extend_symmetric, conj

import Base:conj
conj(S::State) = State(gp=S.gp, u=conj.(S.u))

function coupling_kernel(x::Real, y::Real, σ::Real)
    x = if (x > π)  x - 2π else x end
    y = if (y > π)  y - 2π else y end
    if (x^2 + y^2 ≤ (π * σ)^2)  1 / (π^3 * σ^2) else 0.0 end
end

function get_g_ft(gp::GenericParams)
    xs = getXs(gp)
    fft([coupling_kernel(x, y, gp.p[:σ]) for y in xs, x in xs])
end

get_g_ft(S::State) = get_g_ft(S.gp)

function extend_symmetric(A::Matrix; dim::Int=1)
    cat(A, reverse(A, dims=dim), dims=dim)
end

function RHS_Spiral(S::State, S₀_conj::State, g_ft::Matrix)
    @unpack α, γ, Ω, vx, vy = S.gp.p

    Gz = ifft(g_ft .* fft(S.u)) * (exp(-im * α) * (2π / S.gp.N) ^ 2)
    z2_G_conjz = S.u .^ 2 .* ifft(g_ft .* fft(conj.(S.u))) * (exp(im * α) * (2π / S.gp.N) ^ 2)

    V = vx * ∂ₓ(S) .+ vy * ∂ₓ(S, axis='y') .- (γ + im * Ω) * S.u .+ 0.5 * Gz .- 0.5 * z2_G_conjz

    V1 = - (γ + im * Ω) * S.u
    V2 = vx * ∂ₓ(S) .+ vy * ∂ₓ(S, axis='y')
    V3 = 0.5 * Gz .- 0.5 * z2_G_conjz

    xs = range(-π, π, length=S.gp.N)
    g = real(trapz((xs, xs), S.u .* ∂ₓ(S₀_conj)))
    h = real(trapz((xs, xs), S.u .* ∂ₓ(S₀_conj, axis='y')))
    i = real(trapz((xs, xs), im * S.u .* S₀_conj.u))

    V, V1, V2, V3, [g, h, i]
end

function RHS_Spiral(V::Array, V₀_conj::Array, g_ft::Array, gp::GenericParams)
    S = State(gp=gp, u=reshape(V, (gp.N, gp.N)))
    S₀_conj = State(gp=gp, u=reshape(V₀_conj, (gp.N, gp.N)))
    V, pinning = RHS_Spiral(S, S₀_conj, g_ft)
    append!(vec(V), pinning)
end

function RHS_Spiral_Sym()
    # a
end


text:
using VulkanApp
using Test

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


text:
export @class, @-
export makefunc, MethodDef, Signature, Param

# TODO: Remove. Just added to help debug in REPL
export allocclass, createclass, register

# function createmethod end

"""
    allocclass(name::Union{AbstractString, Symbol}, super::Class) -> Class

Allocate memory for structures to represent a new class named `name` which is
a subclass of `super`. This just creates the class objects. You still need to
register the class returned with `register(class)`.
"""
function allocclass(name::Union{AbstractString, Symbol}, super::Class)
  ptr = ccall(:objc_allocateClassPair, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cchar}, Csize_t),
              super, name, 0)
  ptr == C_NULL && error("Couldn't allocate class $name")
  return Class(ptr)
end

function register(class::Class)
  ccall(:objc_registerClassPair, Cvoid, (Ptr{Cvoid},),
        class)
  return class
end

"""
    createclass(name::Union{AbstractString, Symbol}, super::Class) -> Class
Allocates and register a new class named `name` which is sublcassing `super`.
"""
function createclass(name::Union{AbstractString, Symbol}, super::Class)
    allocclass(name, super) |> register
end

getmethod(class::Class, sel::Selector) =
  ccall(:class_getInstanceMethod, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}),
        class, sel)

methodtypeenc(method::Ptr) =
  ccall(:method_getTypeEncoding, Ptr{Cchar}, (Ptr{Cvoid},),
        method) |> unsafe_string

methodtypeenc(class::Class, sel::Selector) = methodtypeenc(getmethod(class, sel))

methodtype(args...) = methodtypeenc(args...) |> parseencoding

replacemethod(class::Class, sel::Selector, imp::Ptr{Cvoid}, types::String) =
  ccall(:class_replaceMethod, Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cchar}),
        class, sel, imp, types)

function setmethod(class::Class, sel::Selector, imp::Ptr{Cvoid}, types::String)
  meth = getmethod(class, sel)
  meth ≠ C_NULL && methodtype(meth) != parseencoding(types) &&
    error("New method $(name(sel)) of $class must match $(methodtype(meth))")
  replacemethod(class, sel, imp, types)
end

### Syntax ###################

"""
    @class(classname, body)

Define a new Objective-C class with method. Instance methods are marked
with a `@-`. Class methods are not yet supported.

# Examples
```julia-repl
julia> @class Foobar begin
          @- (Cdouble) multiply:(Cdouble)x by:(Cdouble)y begin
            x * y
          end
       end

julia> foo = @objc [[Foobar alloc] init]

julia> @objc [foo multiply:4 by:5]
20.0
```
"""
macro class(cname, body)
    name, supername = class_name(cname)
    if !isexpr(body, :block)
        error("You need to provide a block for class methods")
    end

    if classexists(name)
        classobj = Class(name)
    else
        classobj = createclass(name, Class(supername))
    end

    methodexpressions = filter(body.args) do arg
        isexpr(arg, :macrocall)
    end

    for mexpr in methodexpressions
        sign, fnbody = eval(mexpr)
        mdef = makefunc(sign, fnbody)
        setmethod(classobj, mdef.selector, mdef.cfunc, mdef.typestr)
    end
    :(const $(esc(name)) = Class($(Expr(:quote, name))))
end

function class_name(cname::Expr)
    if cname.head == :(<:)
        (cname.args[1], cname.args[2])
    else
        error("You need to write class name or 'classname <: superclass'")
    end
end

function class_name(cname::Symbol)
    (cname, :NSObject)
end

### Methods ###################

"Represent an ObjC parameter such as `multiply:(Ddouble)x`"
struct Param
    name::Symbol
    kind::Symbol
    value::Symbol
end

"ObjC method signature"
struct Signature
    params::Vector{Param}
    returntype::Symbol
end

"Method definition which we can register with a ObjC class"
struct MethodDef
    selector::Selector  # Selector which will cause method to be called
    typestr::String     # Tells ObjC what type each argument each
    cfunc::Ptr{Cvoid}   # C function pointer to implementation of method
end

"Produce a selector from a signature"
function Selector(sig::Signature)
    map(sig.params) do param
      string(param.name, ':')
    end |> join |> Selector
end

"An arbitrary name to use for Julia functions registered as ObjC method"
funcname(sign::Signature) = Symbol(join([p.name for p in sign.params], '_'))

argtypes(sign::Signature) = (param.kind for param in sign.params)

"""
    argtuple(sig::Signature) -> Expr
Returns a tuple expression with all the types of the arguments in the signature `sig`.
"""
argtuple(sign::Signature) = Expr(:tuple, argtypes(sign)...)

macro -(rettype::Symbol, parts...)
    # string(rettype)
    # dump(params[1])
    fnbody = parts[end]
    params = map(parts[1:end-1]) do part
        Param(
            part.args[2],
            part.args[3].args[2],
            part.args[3].args[3]
        )
    end
    sign = Signature(collect(params), rettype)
    (sign, fnbody)
end

"Create a Julia function with signature `sign` and `body`"
function makefunc(sign::Signature, body::Expr)
    atuple = argtuple(sign)
    params = map(sign.params) do param
        :($(param.value)::$(param.kind))
    end

    # Code in method without line number info
    stmts = filter(body.args) do exp
        !isa(exp, LineNumberNode)
    end

    # NOTE: Not sure why we cannot use method name directly for
    # @cfunction. That is why we store function in fn
    methodname = funcname(sign)
    fn = @eval begin
       function $methodname($(params...))
           $(stmts...)
       end
    end

    # Turn Julia function in a C function
    @eval begin
        imp = @cfunction($fn, $(sign.returntype), $atuple)
        typestr = encodetype(
            $(sign.returntype),
            Object,
            Selector,
            $(argtypes(sign)...))
    end

    MethodDef(Selector(sign), typestr, imp)
end


text:
## An example subtype of AbstractRationalFunction
## that might prove useful as a guide
module TransferFunction
# https://github.com/andreasvarga/DescriptorSystems.jl/blob/main/src/types/RationalFunction.jl

using Polynomials
import Polynomials: AbstractPolynomial, AbstractRationalFunction, RationalFunction
import Polynomials: pqs

export RationalTransferFunction
export sampling_time, gain

struct RationalTransferFunction{T,X,P<:AbstractPolynomial{T,X},Ts} <: AbstractRationalFunction{T,X,P}
    num::P
    den::P
    function RationalTransferFunction{T,X,P,Ts}(num::P, den::P) where{T,X,P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        new{T,X,P,Ts}(num,den)
    end
    function RationalTransferFunction{T,X,P,Ts}(num::P, den::P,ts::Union{Real,Nothing}) where{T,X,P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        check_Ts(Ts,ts)
        new{T,X,P,Ts}(num,den)
    end
    # can promote constants to polynomials too
    function  RationalTransferFunction{T,X,P,Ts}(num::S, den::P, ts::Union{Real,Nothing}) where{S, T,X,P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        check_Ts(Ts,ts)
        new{T,X,P,Ts}(convert(P,num),den)
    end
    function  RationalTransferFunction{T,X,P,Ts}(num::P,den::S, ts::Union{Real,Nothing}) where{S, T,X,P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        check_Ts(Ts,ts)
        new{T,X,P,Ts}(num, convert(P,den))
    end
    function RationalTransferFunction{T,X,P}(num::P, den::P, Ts::Union{Real,Nothing}) where{T,X,P<:AbstractPolynomial{T,X}}

        check_den(den)
        Ts′ = standardize_Ts(Ts)
        new{T,X,P,Val(Ts′)}(num,den)
    end
end

(pq::RationalTransferFunction)(x) = Polynomials.eval_rationalfunction(x, pq)

function Base.convert(::Type{PQ}, pq::RationalTransferFunction) where {PQ <:RationalFunction}
    p,q = pq
    p//q
end

# alternate constructor
function RationalTransferFunction(p′::P, q′::Q, Ts::Union{Real,Nothing}) where {T,X,P<:AbstractPolynomial{T,X},
                                                                                S,  Q<:AbstractPolynomial{S,X}}

    p,q = promote(p′, q′)
    R = eltype(p)
    RationalTransferFunction{R,X,typeof(p)}(p,q,Ts)
end

function Polynomials.rational_function(::Type{PQ}, p::P, q::Q) where {PQ <:RationalTransferFunction,
                                                          T,X,   P<:AbstractPolynomial{T,X},
                                                          S,   Q<:AbstractPolynomial{S,X}}
    RationalTransferFunction(promote(p,q)..., sampling_time(PQ))
end

## helpers for constructors
# standardize Ts or throw error
function standardize_Ts(Ts)
    isnothing(Ts) || Ts >= 0 || Ts == -1 ||
        throw(ArgumentError("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)"))
    Ts′ = isnothing(Ts) ? Ts : Float64(Ts)
end
function check_Ts(Ts, ts)
    ValT(Ts) == promote_Ts(ValT(Ts), ts) || throw(ArgumentError("sampling times have mismatch"))
end
function check_den(den)
    iszero(den) && throw(ArgumentError("Cannot create a rational function with zero denominator"))
end

ValT(::Val{T}) where {T} = T
sampling_time(pq::RationalTransferFunction{T,X,P,Ts}) where {T,X,P,Ts} = ValT(Ts)
sampling_time(::Type{𝑷}) where {T,X,P,Ts, 𝑷<:RationalTransferFunction{T,X,P,Ts}} = ValT(Ts)

## ----

function Base.convert(PQ′::Type{PQ}, p::P) where {PQ <: RationalTransferFunction, P<:AbstractPolynomial}
    PQ(p, one(p), sampling_time(PQ))
end
function Base.convert(PQ′::Type{PQ}, p::Number) where {PQ <: RationalTransferFunction}
    PQ(p, one(eltype(PQ)), sampling_time(PQ))
end

function promote_Ts(p,q)
    Ts1,Ts2 = sampling_time.((p,q))
    promote_Ts(Ts1, Ts2)
end

function promote_Ts(Ts1::Union{Float64,Nothing}, Ts2::Union{Float64,Nothing})
    isnothing(Ts1) && (return Ts2)
    isnothing(Ts2) && (return Ts1)
    Ts1 == Ts2 && (return Ts1)
    Ts1 == -1 && (Ts2 > 0 ? (return Ts2) : error("Sampling time mismatch"))
    Ts2 == -1 && (Ts1 > 0 ? (return Ts1) : error("Sampling time mismatch"))
    error("Sampling time mismatch")
end

function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts},
                                                           T′,X′,P′,Ts′,PQ′ <: RationalTransferFunction{T′,X′,P′,Ts′}}

    S = promote_type(T,T′)
    Polynomials.assert_same_variable(X,X′)
    Y = X
    Q = promote_type(P, P′)
    ts = promote_Ts(PQ, PQ′)
    RationalTransferFunction{S,Y,Q,Val(ts)}
end

# zeros (roots) and poles are defined in common. The gain is specific to this  type
"""
    gain(pq::RationalTransferFunction)

The ratio of the leading coefficients
"""
function gain(pq::PQ) where {PQ <: RationalTransferFunction}
    p,q = pqs(pq)
    p[end]/q[end]
end

# need to check here
#
"""
     rt = adjoint(r)
Compute the adjoint `rt(λ)` of the rational transfer function `r(λ)` such that for
`r(λ) = num(λ)/den(λ)` we have:
    (1) `rt(λ) = conj(num(-λ))/conj(num(-λ))`, if `r.Ts = 0`;
    (2) `rt(λ) = conj(num(1/λ))/conj(num(1/λ))`, if `r.Ts = -1` or `r.Ts > 0`.
"""
function Base.adjoint(pq::PQ) where {PQ <: RationalTransferFunction}
    p,q = pqs(pq)
    Ts = sampling_time(pq)
    if Ts != nothing && iszero(Ts)
        # p(-λ)/q(-λ)
        p′ = poly_scale(p, -1)
        q′ = poly_scale(q, -1)
        return RationalTransferFunction(p′, q′, Ts)
    else
       # p(1/λ) / q(1/λ) = poly_inversion(p) / poly_inversion(q)
       # maps oo -> 0
       p′ = poly_inversion(p)
       q′ = poly_inversion(q)
       return RationalTransferFunction(p′, q′, Ts)
    end
end

## XXX confmap ...

end


text:
import Base.convert

convert(::Type{Int}, v::SubString{String}) = parse(Int, v)

max_stars_repo_name: Fatalon/Dashboard

max_issues_repo_name: Fatalon/Dashboard

max_forks_repo_name: Fatalon/Dashboard


text:
# If the number has 3 as a factor, output 'Pling'.
# If the number has 5 as a factor, output 'Plang'.
# If the number has 7 as a factor, output 'Plong'.
# If the number does not have 3, 5, or 7 as a factor, just pass the number's digits straight through.

function raindrops(number::Int)
    drops = []
    number % 3 == 0 && push!(drops, "Pling")
    number % 5 == 0 && push!(drops, "Plang")
    number % 7 == 0 && push!(drops, "Plong")
    if size(drops, 1) == 0
      repr(number)
    else
      join(drops)
    end
end


text:
"""
    hasrowmajororder(::Type{AbstractSparseMatrix})

Check if values are stored in row-major order.
Return false.
"""
function hasrowmajororder(::Type{AbstractSparseMatrix})
  @abstractmethod
end

function hasrowmajororder(S::AbstractSparseMatrix)
  hasrowmajororder(typeof(S))
end

"""
    hascolmajororder(::Type{AbstractSparseMatrix})

Check if values are stored in col-major order.
Return true.
"""
function hascolmajororder(::Type{AbstractSparseMatrix})
  @abstractmethod
end

function hascolmajororder(S::AbstractSparseMatrix)
  hascolmajororder(typeof(S))
end

"""
    getptr(S::AbstractSparseMatrix)

Return columns pointer.
"""
function getptr(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    getindices(S::AbstractSparseMatrix)

Return row indices.
"""
function getindices(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    rowvals(S::AbstractSparseMatrix)

Return row indices or raises an error (Depending on the SparseMatrix type)
"""
function rowvals(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    colvals(S::AbstractSparseMatrix)

Return columns indices or raises an error (Depending on the SparseMatrix type)
"""
function colvals(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    nnz(S::AbstractSparseMatrix)

Returns the number of stored (filled) elements in a sparse array.
"""
function nnz(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    nonzeros(S::AbstractSparseMatrix)

Return a vector of the structural nonzero values in sparse array S.
This includes zeros that are explicitly stored in the sparse array.
The returned vector points directly to the internal nonzero storage of S,
and any modifications to the returned vector will mutate S as well.
"""
function nonzeros(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    count(pred, S::AbstractSparseMatrix) -> Integer

Count the number of elements in nonzeros(S) for which predicate pred returns true.
"""
function count(pred, S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    nzrange(S::AbstractSparseMatrix, index::Integer) where {Bi}

Return the range of indices to the structural nonzero values of a
sparse matrix index (Row or column depending on the compression type).
"""
function nzrange(S::AbstractSparseMatrix, index::Integer)
  @abstractmethod
end

"""
    findnz(S::AbstractSparseMatrix)

Return a tuple (I, J, V) where I and J are the row and column indices
of the stored ("structurally non-zero") values in sparse matrix A,
and V is a vector of the values.
"""
function findnz(S::AbstractSparseMatrix)
  @abstractmethod
end

"""
    convert(::Type{AbstractSparseMatrix}, x::AbstractSparseMatrix)

Convert x to a value of the first type given.
"""
function convert(::Type{AbstractSparseMatrix}, x::AbstractSparseMatrix)
  @abstractmethod
end

function fill_entries!(A::AbstractSparseMatrix{Tv,Ti},v::Number) where {Tv,Ti}
  nonzeros(A) .= convert(Tv,v)
end


text:
using Justly
using Documenter: deploydocs, makedocs

makedocs(
    sitename = "Justly.jl",
    modules = [Justly],
    doctest = false,
    pages = [
        "Public interface" => "index.md"
    ]
)
deploydocs(repo = "github.com/bramtayl/Justly.jl.git")


text:
# METHOD TEMPLATES
# basic method template: "This method has priority over the more general T:<Union{....} template"
#set_wfm(obj::INSTR{:Agilent33500B},func="sin",ch="1") = write(obj,"SOURce$ch:FUNC $func")
# more general method template, superseded by basic method template
#set_wfm(obj::T;func="sin",ch="1") where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC $func")

## TRIGGER
# Trigger source
set_trig_sour(obj::T;ch="1",sour="BUS") where {T<:Union{F332x1,F335x2}} = write(obj,"TRIG$ch:SOUR $sour")
# Trigger delay
set_trig_del(obj::T;ch="1",del=0) where {T<:Union{F332x1,F335x2}} = write(obj,"TRIG$ch:DEL $del")
# Trigger slope
set_trig_slope(obj::T;ch="1",slp="positive") where {T<:Union{F332x1,F335x2}} = write(obj,"TRIG$ch:SLOP $slp")
# Output trigger state
set_out_trig_stat(obj::T;st="normal") where {T<:Union{F332x1,F335x2}} = write(obj,"OUTP:TRIG $st")
# Output trigger slope
set_out_trig_slope(obj::T;slp="positive") where {T<:Union{F332x1,F335x2}} = write(obj,"OUTP:TRIG:SLOP $slp")
# Software trigger
send_soft_trig(obj::T) where {T<:Union{F332x1,F335x2}} = write(obj,"*TRG")

## OUTPUT FUNCTION
# Waveform
set_wfm(obj::T;ch="1",func="sin") where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC $func")
# ARB Waveform
set_arb_wfm(obj::T;ch="1",arb="mywfm") where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNCtion:ARBitrary $arb")
# Amplitude
set_amplit(obj::T;ch="1",unit="vpp",volt=0.01) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:VOLT:UNIT $unit\n;SOURce$ch:VOLT $volt")
# Offset
set_offs(obj::T;ch="1",volt=0.01) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:VOLT:OFFset $volt")
# High Level
set_hilev(obj::T;ch="1",volt=0.01) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:VOLT:HIGH $volt")
# Low Level
set_lolev(obj::T;ch="1",volt=0.01) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:VOLT:LOW $volt")
# Frequency
set_freq(obj::T;ch="1",fr=1000) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FREQ $fr")
# Period
set_per(obj::T;ch="1",func="sin",per=0.001) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC:$func:PER $per")
# Phase
set_phase(obj::T;ch="1",unit="deg",phase=90) where {T<:Union{F332x1,F335x2}} = write(obj,"UNIT:ANGL $unit\n;SOURce$ch:PHAS $phase")
# Symmetry
set_sym(obj::T;ch="1",func="ramp",sym=50) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC:$func:SYMM $sym")
# Both Edge times
set_both_edg(obj::T;ch="1",func="pulse",dt=10e-9) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC:$func:TRAN $dt")
# Lead Edge time
set_lead_edg(obj::T;ch="1",func="pulse",dt=10e-9) where {T<:F335x2} = write(obj,"SOURce$ch:FUNC:$func:TRAN:LEAD $dt")
# Trail Edge time
set_trail_edg(obj::T;ch="1",func="pulse",dt=10e-9) where {T<:F335x2} = write(obj,"SOURce$ch:FUNC:$func:TRAN:TRA $dt")
# Width
set_width(obj::T;ch="1",func="pulse",wd=1e-4) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC:$func:WIDTh $wd")

## OUTPUT
# Output Load
set_load(obj::T;ch="1",load="INF") where {T<:Union{F332x1,F335x2}} = write(obj,"OUTP$ch:LOAD $load")
# Output Polarity
set_pol(obj::T;ch="1",pol="normal") where {T<:Union{F332x1,F335x2}} = write(obj,"OUTP$ch:POL $pol")
# Output on/off
set_outp(obj::T;ch="1",st="off") where {T<:Union{F332x1,F335x2}} = write(obj,"OUTP$ch $st")

## SYNC
# Sync state
set_sync_stat(obj::T;ch="1",st="on") where {T<:Union{F332x1,F335x2}} = write(obj,"OUTP$ch:SYNC $st")

## ARBITRARY WAVEFORMS
# Arb sample rate
set_sample_rate(obj::T;ch="1",srate=1e6) where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC:ARB:SRAT $srate")
# Clear all arbs from volatile memory
clear_arbs(obj::T;ch="1") where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:DATA:VOLatile:CLEar")


text:
using Printf
function rescale_mesh(infile, outfile; scale=1.0f0)
    scale = Float32(scale)
    fid = open(infile)
    out = open(outfile, write=true)
    for line = eachline(fid)
        if length(line) > 0 && line[1] == "v"
            ldata = split(line, " ")
            data = parse.(Float32, ldata[2:4]) * scale
            scaled_line = @sprintf("%s %.6f %.6f %.6f", ldata[1], data...)
            println(out, scaled_line)
            println(scaled_line)
        else
            println(out, line)
        end
    end
    close(out)
    close(fid)
end


text:
# Generic DiscreteNonParametric with symbol support.
struct GenericDiscreteNonParametric
    support::Any
    distribution::DiscreteNonParametric
end

GenericDiscreteNonParametric(vs::T, ps::Ps) where {
        T<:Any,P<:Real,Ps<:AbstractVector{P}} = GenericDiscreteNonParametric([v for v in vs], DiscreteNonParametric([i for i=1:length(vs)], ps))

Distributions.support(d::GenericDiscreteNonParametric) = d.support

Distributions.probs(d::GenericDiscreteNonParametric)  = d.distribution.p

function Base.rand(rng::AbstractRNG, d::GenericDiscreteNonParametric)
    x = support(d)
    p = probs(d)
    n = length(p)
    draw = rand(rng, float(eltype(p)))
    cp = p[1]
    i = 1
    while cp <= draw && i < n
        @inbounds cp += p[i +=1]
    end
    return x[i]
end

function Distributions.pdf(d::GenericDiscreteNonParametric, x::Any)
    s = support(d)
    idx = findfirst(isequal(x), s)
    ps = probs(d)
    if idx <= length(ps) && s[idx] == x
        return ps[idx]
    else
        return zero(eltype(ps))
    end
end
Distributions.logpdf(d::GenericDiscreteNonParametric, x::Any) = log(pdf(d, x))


text:
using Test
using Dorothy

@DorothyAll()

@testset "NDX" begin

	@testset "Read" begin
		index = readf("$(Dorothy.datapath)/PLC.ndx")
		@test index["system"] == 1:1612871
		@test index["protein"] == 1:81923
		@test index["Ahc_mod1"] == 1:5854
		@test index["B2m_mod1"] == 5855:7587
	end

	@testset "Write" begin
		index1 = readf("$(Dorothy.datapath)/PLC.ndx")
		io = IOBuffer()
		write(specify(io, "text/x-ndx"), index1)
		seekstart(io)
		index2 = read(specify(io, "text/x-ndx"))
		@test index2 == index1
	end

end # @testset


text:
# These are the methods that custom LCM types need to overload.
# The expected signatures are:
# encode(::MyMessageType)::Vector{UInt8}
# decode(::Vector{UInt8}, ::Type{MyMessageType})::MyMessageType
function encode end
function decode end

struct SubscriptionOptions{T, F}
    msgtype::Type{T}
    handler::F
    channel::String
end

struct Subscription{T <: SubscriptionOptions}
    options::T
    csubscription::Ptr{Cvoid}
end
unsafe_convert(::Type{Ptr{Cvoid}}, sub::Subscription) = sub.csubscription

# Troubleshooting adapted from:
# https://github.com/RobotLocomotion/drake/blob/f72bb3f465f69f25459a7a65c57b45c795b5e31d/drake/matlab/util/setup_loopback_multicast.sh
# https://github.com/RobotLocomotion/drake/blob/f72bb3f465f69f25459a7a65c57b45c795b5e31d/drake/matlab/util/check_multicast_is_loopback.sh
function troubleshoot()
    lo = loopback_interface()
    if check_loopback_multicast(lo)
        check_multicast_routing(lo)
    end
    error("Failed to create LCM instance.")
end

loopback_interface() = chomp(read(pipeline(`ifconfig`, `grep -m 1 -i loopback`, `cut -d : -f1`), String))

function loopback_multicast_setup_advice(lo::AbstractString)
    if Sys.isapple()
        """Consider running (as root):
        route add -net 224.0.0.0 -netmask 240.0.0.0 -interface $lo"""
    elseif Sys.islinux()
        """Consider running (as root):
        ifconfig $lo multicast
        route add -net 224.0.0.0 netmask 240.0.0.0 dev lo"""
    else
        "OS-specific instructions not available."
    end
end

function check_loopback_multicast(lo::AbstractString)
    pass = if Meta.parse(read(pipeline(`ifconfig $lo`, `grep -c -i multicast`), String)) != 0
        msg = """Loopback interface $lo is not set to multicast.
        The most probable cause for this is that you are not connected to the internet.
        See https://lcm-proj.github.io/multicast_setup.html.
        $(loopback_multicast_setup_advice(lo))"""
        warn(msg)
        false
    else
        true
    end
    pass
end

function check_multicast_routing(lo::AbstractString)
    routing_correct = if Sys.isapple()
        chomp(read(pipeline(`route get 224.0.0.0 -netmask 240.0.0.0`, `grep -m 1 -i interface`, `cut -f2 -d :`, `tr -d ' '`), String)) == lo
    elseif Sys.islinux()
        chomp(read(pipeline(`ip route get 224.0.0.0`, `grep -m 1 -i dev`, `sed 's/.*dev\s*//g'`, `cut -d ' ' -f1`), String)) == lo
    else
        error("Sorry, I only know how to check multicast routing on Linux and macOS")
    end
    pass = if !routing_correct
        msg = """
        Route to multicast channel does not run through the loopback interface.
        The most probable cause for this is that you are not connected to the internet.
        $(loopback_multicast_setup_advice(lo))
        """
        warn(msg)
        false
    else
        true
    end
    pass
end

# must be mutable so that we can attach a finalizer
mutable struct LCM
    pointer::Ptr{Cvoid}
    provider::String
    filedescriptor::RawFD
    subscriptions::Vector{Subscription}

    LCM(provider="") = begin
        pointer = ccall((:lcm_create, liblcm), Ptr{Cvoid}, (Ptr{UInt8},), provider)
        if pointer == C_NULL
            troubleshoot()
        end
        filedescriptor = RawFD(ccall((:lcm_get_fileno, liblcm), Cint, (Ptr{Cvoid},), pointer))
        lcm = new(pointer, provider, filedescriptor, Subscription[])
        finalizer(close, lcm)
        lcm
    end
end
unsafe_convert(::Type{Ptr{Cvoid}}, lcm::LCM) = lcm.pointer

isgood(lcm::LCM) = lcm.pointer != C_NULL

function close(lcm::LCM)
    if isgood(lcm)
        ccall((:lcm_destroy, liblcm), Cvoid, (Ptr{Cvoid},), lcm)
        lcm.pointer = C_NULL
    end
end

function LCM(func::Function)
    lcm = LCM()
    try
        func(lcm)
    finally
        close(lcm)
    end
end

function publish(lcm::LCM, channel::AbstractString, msg::T) where T
    data = encode(msg)
    publish(lcm, convert(String, channel), data)
end

function publish(lcm::LCM, channel::AbstractString, data::Vector{UInt8})
    status = ccall((:lcm_publish, liblcm), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}, Cuint), lcm, channel, data, length(data))
    return status == 0
end

filedescriptor(lcm::LCM) = lcm.filedescriptor

struct RecvBuf
    data::Ptr{UInt8}
    data_size::UInt32
    recv_utime::Int64
    lcm::Ptr{Cvoid}
end

function check_channel_name(channelptr::Ptr{UInt8}, opts::SubscriptionOptions)
    i = 1
    while (byte = unsafe_load(channelptr, i)) != 0x00
        byte == codeunit(opts.channel, i) || error("Mismatch between received channel name and subscription channel name")
        i += 1
    end
end

function onresponse(rbuf::RecvBuf, channelptr::Ptr{UInt8}, opts::SubscriptionOptions{T}) where T
    check_channel_name(channelptr, opts)
    GC.@preserve rbuf begin
        msgdata = UnsafeArray(rbuf.data, (Int(rbuf.data_size), ))
        if T === Nothing
            opts.handler(opts.channel, msgdata)
        else
            msg = decode(msgdata, opts.msgtype)
            opts.handler(opts.channel, msg)
        end
    end
    return nothing
end

function subscribe(lcm::LCM, options::T) where T <: SubscriptionOptions
    csubscription = ccall((:lcm_subscribe, liblcm), Ptr{Cvoid},
        (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Cvoid}, Ptr{Cvoid}),
        lcm,
        options.channel,
        @cfunction(onresponse, Cvoid, (Ref{RecvBuf}, Ptr{UInt8}, Ref{T})),
        Ref(options))
    sub = Subscription(options, csubscription)
    push!(lcm.subscriptions, sub)
    sub
end

function subscribe(lcm::LCM, channel::String, handler, msgtype=Nothing)
    subscribe(lcm, SubscriptionOptions(msgtype, handler, channel))
end

function unsubscribe(lcm::LCM, subscription::Subscription)
    result = ccall((:lcm_unsubscribe, liblcm), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), lcm, subscription)
    result == 0
end

function set_queue_capacity(sub::Subscription, capacity::Integer)
    @assert capacity >= 0
    status = ccall((:lcm_subscription_set_queue_capacity, liblcm), Cint, (Ptr{Cvoid}, Cint), sub, capacity)
    return status == 0
end

function get_queue_size(sub::Subscription)
    ccall((:lcm_subscription_get_queue_size, liblcm), Cint, (Ptr{Cvoid},), sub)
end

lcm_handle(lcm::LCM) = ccall((:lcm_handle, liblcm), Cint, (Ptr{Cvoid},), lcm)

function handle(lcm::LCM)
    fd = filedescriptor(lcm)
    while isgood(lcm)
        event = poll_fd(fd, 10; readable=true)
        if event.readable
            lcm_handle(lcm)
            return true
        end
    end
end

function handle(lcm::LCM, timeout::Period)
    timeout_ms = Dates.value(convert(Millisecond, timeout))
    fd = filedescriptor(lcm)
    event = poll_fd(fd, timeout_ms / 1000; readable=true)
    if event.readable
        lcm_handle(lcm)
        return true
    else
        return false
    end
end


text:
const IntOrTuple = Union{Int, NTuple{N,Int} where N}

gradtest(f, dims::IntOrTuple...; kw...) =
    gradtest(f, randn.(Ref(rng), Float64, dims)...; kw...) # julia v1.3 compat
    # gradtest(f, randn.(rng, Float64, dims)...; kw...)

"""
Compare numerical gradient and automatic gradient
given by Zygote. `f` has to be a scalar valued function.

Applies also `ChainRulesTestUtils.test_rrule` if the rrule for `f` is explicitly defined.
"""
function gradtest(f, xs...; atol=1e-6, rtol=1e-6, fkwargs=NamedTuple(),
                    check_rrule=false,
                    fdm=:central,
                    check_broadcast=false,
                    skip=false, broken=false)

    if check_rrule
        y = f(xs...; fkwargs...)
        simil(x) = x isa Number ? randn(rng, typeof(x)) : randn!(rng, similar(x))
        ȳ =  simil(y)
        xx̄s = [x ⊢ simil(x) for x in xs]
        test_rrule(f, xx̄s...; fkwargs=fkwargs, output_tangent=ȳ)
    end

    if check_broadcast
        length(fkwargs) > 0 && @warn("CHECK_BROADCAST: dropping keywords args")
        h = (xs...) -> sum(sin.(f.(xs...)))
    else
        h = (xs...) -> sum(sin.(f(xs...; fkwargs...)))
    end

    y_true = h(xs...)
    if fdm == :central
        fdm_obj = FiniteDifferences.central_fdm(5, 1)
    elseif fdm == :forward
        fdm_obj = FiniteDifferences.forward_fdm(5, 1)
    elseif fdm == :backward
        fdm_obj = FiniteDifferences.backward_fdm(5, 1)
    end
    # @show fdm fdm_obj

    gs_fd = FiniteDifferences.grad(fdm_obj, h, xs...)

    y_ad, pull = Zygote.pullback(h, xs...)
    gs_ad = pull(one(y_ad))

    @test y_true ≈ y_ad  atol=atol rtol=rtol
    for (g_ad, g_fd) in zip(gs_ad, gs_fd)
        if skip
            @test_skip g_ad ≈ g_fd   atol=atol rtol=rtol
        elseif broken
            @test_broken g_ad ≈ g_fd   atol=atol rtol=rtol
        else
            @test g_ad ≈ g_fd   atol=atol rtol=rtol
        end
    end
    return true
end


text:
include("./imports.jl")

include("mass_hu.jl")


text:
################################################################################
#
# Set of utilities for TWRIdual
#
################################################################################

export judiVector2vec, vec2judiVector, model2vec, vec2model, contr2abs, abs2contr, gradprec_contr2abs, devito_model, gradient_desc

### JUDI data types/array format translations

function judiVector2vec(y::judiVector)
	return vec(cat(y.data..., dims = 2))
end

function vec2judiVector(rcv_geom::GeometryIC, y_vec::Array)
	nsrc = length(rcv_geom.xloc)
	y_arr = Array{Array}(undef, nsrc)
	for i = 1:nsrc
		nrcv = length(rcv_geom.xloc[i])
		nt = rcv_geom.nt[i]
		y_arr[i] = reshape(y_vec[nrcv*nt*(i-1)+1:nrcv*nt*i], nt, nrcv)
	end
	return judiVector(rcv_geom, y_arr)
end

function model2vec(m)
	return vec(m.m)
end

function vec2model(n, d, o, m)
	return Model(n, d, o, reshape(m, n))
end

### Utilities for model input preprocessing and gradient postprocessing

function contr2abs(xvec::Array{R, 1}, mask::BitArray{2}, mb::Array{R, 2})
# Preprocessing steps: - effective to global domain
#                      - contrast to absolute properties

    x = zeros(R, size(mb))
    x[mask] .= xvec
    return mb.*(R(1).+x)

end
function abs2contr(m::Array{R, 2}, mask::BitArray{2}, mb::Array{R, 2})

	return ((m.-mb)./mb)[mask]

end

function gradprec_contr2abs(grad::Array{R, 2}, mask::BitArray{2}, mb::Array{R, 2})
# Postprocessing steps: - preconditioning by diagonal background model

    return mb[mask].*grad[mask]

end
function gradprec_contr2abs(grad::Array{R, 2}, mask::BitArray{2}, mb::Model)

    return gradprec_contr2abs(grad, mask, mb.m)

end

### Extra devito utilities

function devito_model(model::Model, op::String, mode::Int64, opt, dm)

	pm = load_pymodel()
	length(model.n) == 3 ? dims = [3, 2, 1] : dims = [2, 1] # model dimensions for Python are (z,y,x) and (z,x)

	# Set up Python model structure
	if op == 'J' && mode == 1
		modelPy = pm.Model(origin = model.o, spacing = model.d, shape = model.n,
		vp = process_physical_parameter(sqrt.(1f0./model.m), dims), nbpml = model.nb, rho = process_physical_parameter(model.rho, dims),
		dm = process_physical_parameter(reshape(dm, model.n), dims), space_order = SPACE_ORDER)
	else
		modelPy = pm.Model(origin = model.o, spacing = model.d, shape = model.n, vp = process_physical_parameter(sqrt.(1f0./model.m), dims), nbpml = model.nb,
		rho = process_physical_parameter(model.rho, dims), space_order = SPACE_ORDER)
	end

end

### Optimization routines

function gradient_desc(fun!, x, fun_proj, niter; alpha = 1f0, nFunEval_max = 1000, tolfun = 1f-10, tolg = 1f-10)

	x = fun_proj(x)
	f = 0f0
	g = zeros(R, size(x))
	nFunEval = 0
	t0 = time()
	for i = 1:niter
		f = fun!(true, g, x); nFunEval += 1
		normg = norm(g)
		println(string(i), ": f = ", string(f), " --- g = ", string(normg), ", alpha = ", string(alpha), ", nFunEval = ", string(nFunEval))
		println("* t = ", string(Int(floor(time()-t0))), " s")
		if f < tolfun || normg < tolg
			return x, f
		end
		f_ = f
		while f_ >= f
			if nFunEval > nFunEval_max
				return x, f
			end
			f_ = fun!(true, nothing, fun_proj(x.-alpha*g)); nFunEval += 1
			if f_ >= f
				println(string("      f = ", string(f_), ", halving step length..."))
				alpha /= 2
			end
		end
		x = fun_proj(x.-alpha*g)
		f = f_
		if f < tolfun
			return x, f
		end
	end
	return x, f

end


text:
using Distributions, Turing, Plots, StatsPlots, MCMCChains, DataFrames

function BiExponential(theta1::Real, theta2::Real, p::Real)
    return MixtureModel([Exponential(theta1),
                         Exponential(theta2)], [p, 1 - p])
end

function generate_samples(theta1::Real, theta2::Real, p::Real, n::Int64)
    m = BiExponential(theta1, theta2, p)
    rand(m, n)
end

@model bi_exponential(samples::Array{Real, 1}) = begin
    theta1 ~ Gamma(1, 10)
    theta2 ~ truncated(Normal(), theta1, theta1 + 100)
    p ~ Beta(1, 1)
    m = BiExponential(theta1, theta2, p)
    for i in 1:length(samples)
        samples[i] ~ m
    end
end

function svf(x::Array{Real, 1}, digits::Int64)
    x = round.(x, digits=digits)
    unqx = unique(x) |> sort
    counts = map(x_ -> sum(x .== x_), unqx)
    DataFrame(x = unqx, count = counts, svr = 1. .- cumsum(counts ./ sum(counts)))
end

IRTs = Array{Real, 1}(generate_samples(0.1, 5, 0.9, 2000))
IRTs_svf = svf(IRTs, 2)
plot(IRTs_svf.x, log10.(IRTs_svf.svr))
ylims!(-3, 0)

_ = sample(bi_exponential(IRTs), NUTS(), MCMCThreads(), 2, 4)
chains = sample(bi_exponential(IRTs), NUTS(), MCMCThreads(), 500, 4)
plot(chains)


text:
# Reglarized regressions
#

#######################################
#
#   LASSO
#
#######################################

function lasso(X, y, lambda; t = 0.1, K = 100)
    # Minimizes the objective function:
    # 1 /(2*T) *||y - Xb||^2_2 + lambda * ||b||_1
    #
    # Parameters
    # ----------
    # X : array of training data of dimension T x p
    # y : array of dependent variable of dimension T x 1
    # lambda : regularization parameter 1 x 1
    # t : step size for the proximal coordinate descent
    # K : number of iterations
    T, p = size(X)
    bhat = zeros(p, 1)
    fhat = zeros(K-1, 1)

    for k in 1:K-1
        grad = X' * (X * bhat - y)/T
        v = bhat - t*grad
        bhat = max.(v .- t*lambda, 0) - max.(-v .- t*lambda, 0)
        fhat[k] = 0.5*sum((X*bhat - y).^2)/T + lambda*sum(abs.(bhat))
    end

    return bhat, fhat
end

#######################################
#
#   group LASSO
#
#######################################
function glasso(X, y, lambda, G0, G; t = 0.1, K = 100)
    # Minimizes the objective function:
    # 1 /(2*T) *||y - Xb||^2_2 + lambda * ||b||_{2,1}
    #
    # Parameters
    # ----------
    # X : array of training data of dimension T x p
    # y : array of dependent variable of dimension T x 1
    # lambda : regularization parameter 1 x 1
    # alpha : relative wegiht of LASSO and group-LASSO penalties 1x1
    # G0 : size of the first group
    # G : size of the remaining groups
    # t : step size for the proximal coordinate descent
    # K : number of iterations

    T, p = size(X)
    bhat = zeros(p, 1)
    fhat = zeros(K-1, 1)

    for k in 1:K-1
        grad = X' * (X * bhat - y)/T
        v = bhat - t*grad
        bhat0 = max.(1 - t*lambda ./ sqrt(sum(v[1:G0].^2)), 0) .* v[1:G0]
        penalty = sqrt(sum(bhat0.^2))
        bhat[1:G0] = bhat0
        for g = 1:div(p-G0, G)
            j1 = G0 + (g - 1)*G + 1
            j2 = G0 + g*G
            bhatg = max.(1 - t*lambda ./ sqrt(sum(v[j1:j2].^2)), 0) .* v[j1:j2]
            bhat[j1:j2] = bhatg
            penalty = penalty + sqrt(sum(bhatg.^2))
        end
        fhat[k] = 0.5*sum((X*bhat - y).^2)/T + lambda*penalty
    end

    return bhat, fhat
end

#######################################
#
#   sparse-group LASSO
#
#######################################
function sglasso(X, y, lambda, alpha, G0, G; t = 0.1, K = 100)
    # Minimizes the objective function:
    # 1 /(2*T) *||y - Xb||^2_2 + lambda * alpha * ||b||_1
    # + lambda * (1-alpha) * ||b||_{2,1}
    #
    # Parameters
    # ----------
    # X : array of training data of dimension T x p
    # y : array of dependent variable of dimension T x 1
    # lambda : regularization parameter 1 x 1
    # alpha : relative wegiht of LASSO and group-LASSO penalties 1x1
    # G0 : size of the first group
    # G : size of the remaining groups
    # t : step size for the proximal coordinate descent
    # K : number of iterations

    T, p = size(X)
    bhat = zeros(p, 1)
    fhat = zeros(K-1, 1)

    for k in 1:K-1
        grad = X' * (X * bhat - y)/T
        v = bhat - t*grad
        w = max.(v .- t*lambda*alpha, 0) - max.(-v .- t*lambda*alpha, 0)
        bhat0 = max.(1 - t*lambda*(1-alpha) ./ sqrt(sum(w[1:G0].^2)), 0) .* w[1:G0]
        penalty = sqrt(sum(bhat0.^2))
        bhat[1:G0] = bhat0
        for g = 1:div(p-G0, G)
            j1 = G0 + (g - 1)*G + 1
            j2 = G0 + g*G
            bhatg = max.(1 - t*lambda*(1-alpha) ./ sqrt(sum(w[j1:j2].^2)), 0) .* w[j1:j2]
            bhat[j1:j2] = bhatg
            penalty = penalty + sqrt(sum(bhatg.^2))
        end
        fhat[k] = 0.5*sum((X*bhat - y).^2)/T + lambda*(alpha*sum(abs.(bhat)) + (1-alpha)*penalty)
    end

    return bhat, fhat
end


text:
function testuniverse(;N=560, M=5000)::Universe
  local cols::Vector{Symbol} = (n->Symbol("X",n)).(1:M)
  local crosscor::MFloat64 = 0.25
  local stdev::MFloat64 = 0.1
  local allmean = .02

  ρ::Matrix{MFloat64} = ones(M,M) .* crosscor
  ρ[diagind(ρ)] .= 1.0 #fix the diagonals
  σ::Vector{MFloat64} = (n->stdev).(1:M)
  Σ::Matrix{MFloat64} = ρ .* stdev^2
  L = cholesky(Σ).L

  #make the test matrix
  local testmatrix::Matrix{MFloat64} = rand(Normal(), N, M)

  for r ∈ 1:N
    testmatrix[r,:] .= L * testmatrix[r,:] .+ allmean
  end

  #make missing values
  for c ∈ 1:M
    for r ∈ 1:N
      (rand(1:10) == 1) && (testmatrix[r,c] = missing)
    end
  end

  #make the dataframe
  df = DataFrame(testmatrix)
  names!(df, cols)
  df[:date] = Vector{MDate}(missing, N)

  local datebegin::Date = Date(2011,11,11)

  df[:date] .= datebegin .+ (Day).(1:N)

  U = Universe(df)

  return U
end

#testing
#simulates an artificial portfolio and reweights it iter times
function testportfolio(;N=560, M=5000, iter=1,
  U::Universe = testuniverse(N=N, M=M))::Portfolio

  wₜ = ones(M) ./ M

  P = Portfolio(U,wₜ)

  local μs::Vector{Float64} = Vector{Float64}(undef, iter)

  for i ∈ 1:iter
    wₜ .= rand(Normal(2,1), M)
    wₜ ./= sum(wₜ)
    Portfolio!(P, wₜ, checkadjustment=true)
    μs[i] = P.pμ[]
  end

  iter>1 && println("Sum of means: $(sum(μs))")

  return P
end

function testcov(;P::Portfolio=testportfolio(iter=1), iter = 1000)::Nothing
  local P1::Portfolio = deepcopy(P)
  local P2::Portfolio = deepcopy(P)
  local N::Int = P.U.N

  covs::Vector{Float64} = Vector{Float64}(undef, iter)
  @time for i::Int ∈ 1:iter
    P1.R .= rand(N)
    P2.R .= rand(N)

    P1.pμ[] = mean(P1.R)
    P2.pμ[] = mean(P2.R)

    covs[i] = cov(P1,P2)
    var(P1)
  end

  print("Sum of covs: $(sum(covs))")

  return nothing
end

function testtrisectedportfolio(;P::Portfolio=testportfolio(iter=1),
  iter = 10, increment = 2, displayiter = 3)

  TP = TrisectedPortfolio(P)
  if displayiter > 0
    display(TP.assignments)
    display(TP.wₛ)
    display(TP.Σ)
  end

  @time for i ∈ 1:iter
    rotatetrisected!(TP, increment)
    if i ≤ displayiter
      display(TP.assignments)
      display(TP.wₛ)
      display(TP.Σ)
    end
  end

  return nothing
end

function testω(;iter::Int = 1)

  errSG::Vector{Float64} = Vector{Float64}(undef, iter)
  errSGP::Vector{Float64} = Vector{Float64}(undef, iter)
  errSGPratio::Vector{Float64} = Vector{Float64}(undef, iter)
  errsum::Vector{Float64} = Vector{Float64}(undef, iter)

  U::Universe=testuniverse()
  G::Portfolio = testportfolio(iter=1, U=U)
  P::Portfolio = initializetestportfolio(G)
  TG::TrisectedPortfolio = TrisectedPortfolio(G)

  P = initializetestportfolio(G)

  @time for i ∈ 1:iter

    initializetestportfolio(G)
    redrawtestportfolio!(P)
    refreshtrisected!(TG)
    #print("var(TG.P: $(var(TG.P)), 1'S1: $(ones(3)' * TG.S * ones(3))\n\n")

    Θ::Parameters = Parameters(
      σ²G = D_σ²G,
      ζ²G = D_ζ²G,
      ζ²P = D_ζ²P,
      SG = var(G),
      SGP = cov(G,P)
    )

    oldS::Matrix{Float64} =  deepcopy(TG.S) #need thsi for checks later

    SGP::Vector{Float64} = cov(TG, P)
    ω::Vector{Float64} = getω(Θ, TG, SGP)
    rescaletrisected!(TG, ω)

    errSG[i] = abs( ω' * oldS * ω - var(TG.P))
    errSGP[i] = abs( ω' * SGP - cov(TG.P,P))
    errSGPratio[i] = abs( cov(TG.P,P) / (ω' * SGP))
    errsum[i] = abs( 1.0 - sum(TG.P.wₜ))

    (iter==1) && print("SG: $(Θ.SG)\n")
    (iter==1) && print("= ω' S ω: $(ω' * oldS * ω)\n")
    (iter==1) && print("= var(TG.P): $(var(TG.P))\n\n")

    (iter==1) && print("SGP: $(Θ.SGP)\n")
    (iter==1) && print("= ω' SGP: $(ω' * SGP)\n")
    (iter==1) && print("= cov(TG.P,P): $(cov(TG.P,P))\n\n")

    (iter==1) && print("1.0 \n")
    (iter==1) && print("= sum(TG.P.wₜ): $(sum(TG.P.wₜ))\n\n")

  end

  if (iter > 1) && (iter<100)
    print("Error SG - w'Sw: $errSG\n")
    print("Error SGP - w'SGP: $errSGP\n")
    print("Error sum 1'wt: $errsum\n")
  elseif iter ≥ 100
    print("Error SG - w'Sw max: $(maximum(errSG)) mean: $(mean(errSG))\n")
    print("Error SGP - w'SGP max: $(maximum(errSGP)) mean: $(mean(errSGP))\n")
    print("Error sum 1'wt max: $(maximum(errsum)) mean: $(mean(errsum))\n")
  end

end

function testposteriors()
  Π::Priors = Priors()
  Θ1::Parameters = Parameters()
  Θ2::Parameters = Parameters()
  F::Posteriors = Posteriors()

  Random.seed!(11)
  F.σ²G!(Θ1, Π)
  F.ζ²G!(Θ1, Π)
  F.ζ²P!(Θ1, Π)
  F.SG!(Θ1, Π)
  F.SGP!(Θ1, Π)

  F = PosteriorsTest()
  Random.seed!(11)
  F.σ²G!(Θ2, Π)
  F.ζ²G!(Θ2, Π)
  F.ζ²P!(Θ2, Π)
  F.SG!(Θ2, Π)
  F.SGP!(Θ2, Π)

  println("Check Θ1: $(println(Θ1))")
  println("Check Θ2: $(println(Θ2))")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text: _url_from_file(path) = string("file://", path)


text:
# This script create default configuration files that contain the tuning parameters
# for some of the algorithms this package. We assign a default value for each
# compound name key entry in the user-specified
# "compound to parameters filename" JSON file.

# One should manually edit the generated config files according to
# their preferred tuning parameters for the different compounds.

import JSON
import JSON3

name_map_filename = "select_compounds"
load_path = joinpath("/home/roy/Documents/repo/NMRData/input/compound_mapping", "$(name_map_filename).json")
dict_compound_to_filename = JSON.parsefile(load_path)
compound_names = collect( key for (key,value) in dict_compound_to_filename )

##### config file for SH simulation.
function createdefaultSHconfig(compound_names::Vector{String};
    tol_coherence::Float64 = 1e-2,
    α_relative_threshold::Float64 = 0.05,
    Δc_partition_radius::Float64 = 0.3,
    simple_coherence_atol::Float64 = -1.2)

    db_dict = Dict()
    for name in compound_names

        db_dict[name] = Dict("coherence tolerance" => tol_coherence,
            "relative amplitude threshold" => α_relative_threshold,
            "maximum Δc deviation" => Δc_partition_radius,
            "simple coherence absolute tolerance" => simple_coherence_atol)
    end

    return db_dict
end

save_folder = "/home/roy/Documents/repo/NMRData/input/SH_configs"
save_name = "$(name_map_filename)_SH_configs_default.json"
save_path = joinpath(save_folder, save_name)

# write the default SH configurations to json.
dict_out = createdefaultSHconfig(compound_names)

stringdata = JSON.json(dict_out)

open(save_path, "w") do f
    JSON3.pretty(f, stringdata)
    println(f)
end
#@assert 1==2

##### config file for SH parameters.
function createdefaultsurrogateconfig(compound_names::Vector{String};
    Δr = 1.0,
    Δκ_λ = 0.05,
    Δcs_max::Vector{Float64} = Vector{Float64}(undef, 0),
    κ_λ_lb = 0.5,
    κ_λ_ub = 2.5)

    db_dict = Dict()
    for name in compound_names

        db_dict[name] = Dict("radians surrogate sampling step size" => Δr,
            "κ_λ surrogate sampling step size" => Δκ_λ,
            "Δcs_max" => Δcs_max,
            "κ_λ lower bound" => κ_λ_lb,
            "κ_λ upper bound" => κ_λ_ub)
    end

    return db_dict
end

save_folder = "/home/roy/Documents/repo/NMRData/input/surrogate_configs"
save_name = "$(name_map_filename)_surrogate_configs_default.json"
save_path = joinpath(save_folder, save_name)

# write the default SH configurations to json.
dict_out = createdefaultsurrogateconfig(compound_names)

stringdata = JSON.json(dict_out)

open(save_path, "w") do f
    JSON3.pretty(f, stringdata)
    println(f)
end


text:
include("h_file.jl")
const EEPORT = 8000

function http_gatekeeper(req)
    if req.method == "GET"
        handle_file(req)
        # Add a Pragma header (should not be stripped by proxies)?
        # Or better, inject in document?
    else
        Response(501, "Unimplemented method on this server: $(req.method)")
    end
end

WebSockets.addsubproto("relay_frontend") # For testing the template file
WebSockets.addsubproto("pagecontrol") # Use one sub protocol per object from Revise?
function ws_gatekeeper(req, ws)
    @info "ws_gatekeeper", req
    # Check against value sent by this server....
    if WebSockets.subprotocol(req) == "pagecontrol"
        coroutine_pagecontrol(ws)
    end
end

function serveEE(logger= WebSocketLogger())
    server = ServerWS(http_gatekeeper, ws_gatekeeper)
    task = @async   with_logger(logger) do
        WebSockets.serve(server, port = EEPORT)
    end
    @info "http://localhost:$EEPORT"
    return server, task
end

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"

using SymPyCall
using Documenter

makedocs(sitename="My Documentation")

#DocMeta.setdocmeta!(SymPyCall, :DocTestSetup, :(using SymPyCall); recursive=true)

# makedocs(;
#     modules=[SymPyCall],
#     authors="jverzani <jverzani@gmail.com> and contributors",
#     repo="https://github.com/jverzani/SymPyCall.jl/blob/{commit}{path}#{line}",
#     sitename="SymPyCall.jl",
#     format=Documenter.HTML(;
#         prettyurls=get(ENV, "CI", "false") == "true",
#         canonical="https://jverzani.github.io/SymPyCall.jl",
#         assets=String[],
#     ),
#     pages=[
#         "Home" => "index.md",
#         "Examples" => "introduction.md"
#     ],
# )

# # deploydocs(;
#     repo="github.com/jverzani/SymPyCall.jl",
#     devbranch="main",
# )

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
#TODO rename to CarbonAssimilation or so? could be consistently named as CarbonPartition, CarbonAllocation...
@system Photosynthesis begin
    weather ~ hold
    sun ~ hold
    soil ~ hold

    LAI: leaf_area_index ~ hold
    PD: planting_density ~ hold
    water_supply ~ hold
    H2O_weight ~ hold
    CO2_weight ~ hold
    CH2O_weight ~ hold

    radiation(context, sun, leaf_area_index) ~ ::Radiation

    # Calculating transpiration and photosynthesis with stomatal controlled by leaf water potential LeafWP Y
    #TODO: use leaf_nitrogen_content, leaf_width, ET_supply
    sunlit_gasexchange(context, soil, weather, PPFD=Q_sun, LAI=LAI_sunlit) ~ ::GasExchange
    shaded_gasexchange(context, soil, weather, PPFD=Q_sh, LAI=LAI_shaded) ~ ::GasExchange

    leaf_width => begin
        # to be calculated when implemented for individal leaves
        #5.0 # for maize
        1.5 # for Radish
    end ~ preserve(u"cm", parameter)

    #TODO how do we get LeafWP and ET_supply?
    LWP(soil.WP_leaf): leaf_water_potential ~ track(u"MPa")

    evapotranspiration_supply(LAI, PD, ws=water_supply, ww=H2O_weight) => begin
        #TODO common handling logic for zero LAI
        #FIXME check unit conversion (w.r.t water_supply)
        # ? * (1/m^2) / (3600s/hour) / (g/umol) / (cm^2/m^2) = mol/m^2/s H2O
        # ? * (1/m^2) * (hour/3600s) * (umol/g) * (m^2/cm^2) = mol/m^2/s H2O
        # ?(g / hour) * (hour/3600s) * (umol/g) / cm^2
        s = ws * PD / 3600 / ww / LAI
        iszero(LAI) ? zero(s) : s
    end ~ track(u"mol/m^2/s" #= H2O =#)

    LAI_sunlit(radiation.sunlit_leaf_area_index): sunlit_leaf_area_index ~ track
    LAI_shaded(radiation.shaded_leaf_area_index): shaded_leaf_area_index ~ track

    Q_sun(radiation.irradiance_Q_sunlit): sunlit_irradiance ~ track(u"μmol/m^2/s" #= Quanta =#)
    Q_sh(radiation.irradiance_Q_shaded): shaded_irradiance ~ track(u"μmol/m^2/s" #= Quanta =#)

    A_gross(a=sunlit_gasexchange.A_gross_total, b=shaded_gasexchange.A_gross_total): gross_CO2_umol_per_m2_s => begin
        a + b
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    # plantsPerMeterSquare units are umol CO2 m-2 ground s-1
    # in the following we convert to g C plant-1 per hour
    # photosynthesis_gross is umol CO2 m-2 leaf s-1

    A_net(a=sunlit_gasexchange.A_net_total, b=shaded_gasexchange.A_net_total): net_CO2_umol_per_m2_s => begin
        # grams CO2 per plant per hour
        a + b
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    ET(a=sunlit_gasexchange.E_total, b=shaded_gasexchange.E_total): transpiration_H2O_mol_per_m2_s => begin
        #TODO need to save this?
        # when outputting the previous step transpiration is compared to the current step's water uptake
        #self.transpiration_old = self.transpiration
        #FIXME need to check if LAIs are negative?
        #transpiration = sunlit_gasexchange.ET * max(0, sunlit_LAI) + shaded_gasexchange.ET * max(0, shaded_LAI)
        a + b
    end ~ track(u"mmol/m^2/s" #= H2O =#)

    # final values
    assimilation(gross_assimilation) ~ track(u"g/d")

    gross_assimilation(A_gross, PD, w=CH2O_weight) => begin
        # grams carbo per plant per hour
        #FIXME check unit conversion between C/CO2 to CH2O
        A_gross / PD * w
    end ~ track(u"g/d")

    net_assimilation(A_net, PD, w=CH2O_weight) => begin
        # grams carbo per plant per hour
        #FIXME check unit conversion between C/CO2 to CH2O
        A_net / PD * w
    end ~ track(u"g/d")

    transpiration(ET, PD, w=H2O_weight) => begin
        # Units of Transpiration from sunlit->ET are mol m-2 (leaf area) s-1
        # Calculation of transpiration from ET involves the conversion to gr per plant per hour
        ET / PD * w
    end ~ track(u"g/d")

    vapor_pressure_deficit(weather.VPD) ~ track(u"kPa")

    conductance(gs_sun=sunlit_gasexchange.gs, LAI_sunlit, gs_sh=shaded_gasexchange.gs, LAI_shaded, LAI) => begin
        #HACK ensure 0 when one of either LAI is 0, i.e., night
        # average stomatal conductance Yang
        c = ((gs_sun * LAI_sunlit) + (gs_sh * LAI_shaded)) / LAI
        #c = max(zero(c), c)
        iszero(LAI) ? zero(c) : c
    end ~ track(u"mol/m^2/s/bar")
end


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

module NetworkInterfaceClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_NetworkInterfacesApi.jl")

# export models
export convert, EffectiveNetworkSecurityGroup
export convert, EffectiveNetworkSecurityGroupAssociation
export convert, EffectiveNetworkSecurityGroupListResult
export convert, EffectiveNetworkSecurityRule
export convert, EffectiveRoute
export convert, EffectiveRouteListResult
export convert, IPConfigurationPropertiesFormat
export convert, NetworkInterfaceDnsSettings
export convert, NetworkInterfaceIPConfigurationPropertiesFormat
export convert, NetworkInterfaceListResult
export convert, NetworkInterfacePropertiesFormat
export convert, Resource
export convert, SubResource
export convert, ApplicationGatewayBackendAddressPool
export convert, BackendAddressPool
export convert, IPConfiguration
export convert, InboundNatRule
export convert, NetworkInterface
export convert, NetworkInterfaceIPConfiguration
export convert, NetworkSecurityGroup
export convert, PublicIPAddress
export convert, Subnet

# export operations
export convert, NetworkInterfacesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
"""
Abstract coarse quantizer type. The coarse quantizer is ment to
assign a point to a given number of Voronoi cells in which its
neighbors will be searched.
"""
abstract type AbstractCoarseQuantizer{D,T} end

# Naive coarse quantizer
"""
    NaiveQuantizer{D<:Distances.PreMetric, T<:AbstractFloat}

Coarse quantization structure based on brute force search.
The `vector` fields contains the coarse vectors while `distance`
contains the distance that is used to calculate the distance
from a point to the coarse vectors.
"""
struct NaiveQuantizer{D<:Distances.PreMetric, T<:AbstractFloat} <: AbstractCoarseQuantizer{D,T}
    vectors::Matrix{T}
end

Base.show(io::IO, cq::NaiveQuantizer{D,T}) where {D,T} = begin
    nrows, nclusters = size(cq.vectors)
    print(io, "NaiveQuantizer{$D,$T}, $nrows×$nclusters cluster centres")
end

Base.size(cq::NaiveQuantizer{D,T}) where {D,T} = size(cq.vectors)
Base.size(cq::NaiveQuantizer{D,T}, i::Int) where {D,T} = size(cq.vectors, i)

function coarse_search(cq::NaiveQuantizer{D,T}, point::Vector{T}, w::Int) where {D,T}
    coarse_distances = colwise(D(), cq.vectors, point)
    closest_clusters = sortperm(coarse_distances)[1:w]
    return closest_clusters, coarse_distances[closest_clusters]
end

_closest_cluster_residuals(cq::NaiveQuantizer{D,T},
                           point::Vector{T},
                           closest_clusters::Vector{Int}
                          ) where {D,T} = begin
    return point .- cq.vectors[:, closest_clusters]
end

_get_quantizer_vector(cq::NaiveQuantizer, idx::Int) = cq.vectors[:, idx]

# HNSW coarse quantizer
"""
    HNSWQuantizer{U<:Unsigned, V<:Vector{Vector{T}}, D<:Distances.PreMetric, T<:AbstractFloat}

Coarse quantization structure based on HNSW search structure.
The `hnsw` field contains the coarse vectors.
"""
struct HNSWQuantizer{U, V, D<:Distances.PreMetric, T<:AbstractFloat} <: AbstractCoarseQuantizer{D,T}
    hnsw::HierarchicalNSW{U,T,V,D}
end

Base.show(io::IO, cq::HNSWQuantizer{U,V,D,T}) where {U,V,D,T} = begin
    nrows, nclusters = size(cq)
    print(io, "HNSWQuantizer{$U,$V,$D,$T}, $nrows×$nclusters cluster centres")
end

Base.size(cq::HNSWQuantizer{D,T}) where {D,T} = length(cq.hnsw.data[1]), length(cq.hnsw.data)
Base.size(cq::HNSWQuantizer{D,T}, i::Int) where {D,T} = size(cq)[i]

function coarse_search(cq::HNSWQuantizer{U,V,D,T}, point::Vector{T}, w::Int) where {U,V,D,T}
    closest_clusters, coarse_distances = HNSW.knn_search(cq.hnsw, point, w)
    return Int.(closest_clusters), coarse_distances
end

_closest_cluster_residuals(cq::HNSWQuantizer{U,V,D,T},
                           point::Vector{T},
                           closest_clusters::Vector{Int}
                          ) where {U,V,D,T} = begin
    nclusters = length(closest_clusters)
    residuals = Matrix{T}(undef, size(point,1), nclusters)
    @inbounds @simd for i in eachindex(closest_clusters)
        residuals[:, i] = point - cq.hnsw.data[closest_clusters[i]]
    end
    return residuals
end

_get_quantizer_vector(cq::HNSWQuantizer, idx::Int) = cq.hnsw.data[idx]


text:
using JuSwarm, StaticArrays, Test
import Optim

function rastriginfunc(x::AbstractVector)
    d = length(x)

    sum = 0.0
    for i in 1:d
        sum += (x[i]^2 - 10*cos(2*π*x[i]))
    end
    return 10*d + sum
end

N = 2
opts = PSOOptions(N; SwarmSize = 200,
                     LowerBounds = SVector{N}(fill(-5.12, N, 1)),
                     UpperBounds = SVector{N}(fill(5.12, N, 1)),
                     HybridOptimizer = Optim.LBFGS(),
                     Display = "none",
                     UseParallel = true)

sol = psoptimize(rastriginfunc,opts)

@test sol.fbest <= sol.fbestPSO

@test sol.fbestPSO <= 1.0e-6

@test sol.fbest <= 1.0e-10

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm