text:
module AffinityMaps
using ..Types
using ..Domains
using ..Segmentations

export aff2seg, exchangeaffxz!, aff2uniform, gaff2saff, aff2edgelist
export maskaff!, mask_margin!

function downsample(aff::AffinityMap; scale::Union{Vector, Tuple} = (2,2,1,1))
    println("downsampling affinitymap with averaging")
    outSize = map(div, size(aff), scale)
    out = similar(aff, (outSize...,))
    for z in 1:outSize[3]
        for y in 1:outSize[2]
            for x in 1:outSize[1]
                out[x,y,z,1:3] = mean(aff[(x-1)*scale[1]+1:x*scale[1],
                                          (y-1)*scale[2]+1:y*scale[2],
                                          (z-1)*scale[3]+1:z*scale[3],
                                          1:3])
            end
        end
    end
    return out
end

"""
transform google affinity to seung lab affinity
"""
function gaff2saff( gaff::Array{Float32,3} )
    @assert ndims(gaff)==3
    sx,sy,sz = size(gaff)
    saff = reshape(gaff, (sx,sy,Int64(sz/3),Int64(3)));

    # transform the x y and z channel
    ret = zeros(size(saff))
    ret[2:end,:,:, 1] = saff[1:end-1,:,:, 1]
    ret[:,2:end,:, 2] = saff[:,1:end-1,:, 2]
    ret[:,:,2:end, 3] = saff[:,:,1:end-1, 3]

    return ret
end

# exchang X and Z channel of affinity
function exchangeaffxz!(aff::AffinityMap)
    println("exchange x and z of affinity map")
    taffx = deepcopy(aff[:,:,:,1])
    aff[:,:,:,1] = deepcopy(aff[:,:,:,3])
    aff[:,:,:,3] = taffx
    #aff[:,:,:,1], aff[:,:,:,3] = aff[:,:,:,3], aff[:,:,:,1]
    return aff
end

# transform affinity to segmentation
function aff2seg( aff::Array{Float32,4}; dim::Integer = 3, thd::Float32 = Float32(0.5) )
    @assert dim==2 || dim==3
    # note that should be column major affinity map
    # the znn V4 output is row major!!! should exchangeaffxz first!
    xaff = aff[:,:,:,1]
    yaff = aff[:,:,:,2]
    zaff = aff[:,:,:,3]

    # number of voxels in segmentation
    N = UInt32( length(xaff) )

    # initialize and create the disjoint sets
    djsets = Tdjsets( N )

    # union the segments by affinity edges
    X,Y,Z = size( xaff )
    X = UInt32(X);   Y = UInt32(Y);   Z = UInt32(Z);

    # x affinity
    for z in UInt32(1):UInt32(Z)
        for y in UInt32(1):UInt32(Y)
            for x in UInt32(2):UInt32(X)
                if xaff[x,y,z] > thd
                    vid1 = x            + (y-UInt32(1))*X + (z-UInt32(1))*X*Y
                    vid2 = x-UInt32(1)  + (y-UInt32(1))*X + (z-UInt32(1))*X*Y
                    rid1 = find!(djsets, vid1)
                    rid2 = find!(djsets, vid2)
                    union!(djsets, rid1, rid2)
                end
            end
        end
    end

    # y affinity
    for z in UInt32(1):UInt32(Z)
        for y in UInt32(2):UInt32(Y)
            for x in UInt32(1):UInt32(X)
                if yaff[x,y,z] > thd
                    vid1 = x + (y-UInt32(1))*X + (z-UInt32(1))*X*Y
                    vid2 = x + (y-UInt32(2))*X + (z-UInt32(1))*X*Y
                    rid1 = find!(djsets, vid1)
                    rid2 = find!(djsets, vid2)
                    union!(djsets, rid1, rid2)
                end
            end
        end
    end

    # z affinity
    if dim > 2
        # only computed in 3D case
        for z in UInt32(2):UInt32(Z)
            for y in UInt32(1):UInt32(Y)
                for x in UInt32(1):UInt32(X)
                    if zaff[x,y,z] > thd
                        vid1 = x + (y-UInt32(1))*X + (z-UInt32(1))*X*Y
                        vid2 = x + (y-UInt32(1))*X + (z-UInt32(2))*X*Y
                        rid1 = find!(djsets, vid1)
                        rid2 = find!(djsets, vid2)
                        union!(djsets, rid1, rid2)
                    end
                end
            end
        end
    end

    # get current segmentation
    setallroot!( djsets )
    # marking the singletones as boundary
    # copy the segment to avoid overwritting of djsets
    seg = deepcopy( djsets.sets )
    seg = reshape(seg, size(xaff) )
    singleton2boundary!( seg )
    return seg
end

function aff2uniform(aff, alg=QuickSort)
    print("map to uniform distribution...")
    tp = typeof(aff)
    sz = size(aff)
    N = length(aff)

    # get the indices
    print("get the permutation by sorting......")
    @time p = sortperm(aff[:], alg=alg)
    println("done :)")
    q = zeros(eltype(p), size(p))
    q[p[1:N]] = 1:N

    # generating values
    v = range(0, stop=1, length=N)
    # making new array
    v = v[q]
    v = reshape(v, sz)
    v = tp( v )
    println("done!")
    return v
end

#function aff2uniform!(aff::AffinityMap)
 #   println("transfer to uniform distribution...")
 #   for z in 1:size(aff,3)
 #       aff[:,:,z,:] = arr2uniform( aff[:,:,z,:] )
 #   end
#end

"""
transfer affinity map to edge list
"""
function aff2edgelist(aff::Array{T,4}; is_sort::Bool=true) where T
    # initialize the edge list
    elst = Array{Tuple{Float32,UInt32,UInt32},1}([])
    sizehint!(elst, div(length(aff),3))
    # get the sizes
    sx,sy,sz,sc = size(aff)
    @assert sc==3

    for z in UInt32(1):UInt32(sz)
        for y in UInt32(1):UInt32(sy)
            for x in UInt32(1):UInt32(sx)
                vid1 = x + (y-1)*sx + (z-1)*sx*sy
                # x affinity
                if x>1
                    vid2 = x-1 + (y-1)*sx + (z-1)*sx*sy
                    push!(elst, (aff[x,y,z,1], vid1, vid2))
                end
                # y affinity
                if y>1
                    vid2 = x + (y-2)*sx + (z-1)*sx*sy
                    push!(elst, (aff[x,y,z,2], vid1, vid2))
                end
                # z affinity
                if z>1
                    vid2 = x + (y-1)*sx + (z-2)*sx*sy
                    push!(elst, (aff[x,y,z,3], vid1, vid2))
                end
            end
        end
    end
    if is_sort
        @time sort!(elst, rev=true)
    end
    return elst
end

"""
    maskaff!(img::EMImage, aff::AffinityMap)

set the affinity edge value to 0 if the connecting voxel is 0 in image
"""
function maskaff!(img::EMImage, aff::AffinityMap)
    @assert size(img) == size(aff[:,:,:,1])
    # mask the affinity
    for z in 2:size(img, 3)
        for y in 1:size(img, 2)
            for x in 1:size(img, 1)
                if img[x,y,z]==0x00 || img[x,y,z-1]==0x00
                    aff[x,y,z,3] = 0.0f0
                end
            end
        end
    end

    for z in 1:size(img, 3)
        for y in 2:size(img, 2)
            for x in 1:size(img, 1)
                if img[x,y,z]==0x00 || img[x,y-1,z]==0x00
                    aff[x,y,z,2] = 0.0f0
                end
            end
        end
    end

    for z in 1:size(img, 3)
        for y in 1:size(img, 2)
            for x in 2:size(img, 1)
                if img[x,y,z]==0x00 || img[x-1,y,z]==0x00
                    aff[x,y,z,1] = 0.0f0
                end
            end
        end
    end
end

function maskaff!(mask::Array{Bool,3}, aff::AffinityMap)
    @assert size(mask) == size(aff[:,:,:,1])
    # mask the affinity
    for z in 2:size(mask, 3)
        for y in 1:size(mask, 2)
            for x in 1:size(mask, 1)
                if mask[x,y,z] || mask[x,y,z-1]
                    aff[x,y,z,3] = 0.0f0
                end
            end
        end
    end

    for z in 1:size(mask, 3)
        for y in 2:size(mask, 2)
            for x in 1:size(mask, 1)
                if mask[x,y,z] || mask[x,y-1,z]
                    aff[x,y,z,2] = 0.0f0
                end
            end
        end
    end

    for z in 1:size(mask, 3)
        for y in 1:size(mask, 2)
            for x in 2:size(mask, 1)
                if mask[x,y,z] || mask[x-1,y,z]
                    aff[x,y,z,1] = 0.0f0
                end
            end
        end
    end
end

function mask_margin!(aff::AffinityMap, maskSize::Vector)
    sx,sy,sz,sc = size(aff)
    for z in 1:sz
        for y in 1:sy
            for x in 1:sx
                if  z<=maskSize[3] || z>=sz-maskSize[3] ||
                    y<=maskSize[2] || y>=sy-maskSize[2] ||
                    x<=maskSize[1] || x>=sx-maskSize[1]
                    aff[x,y,z,:] = 0f0
                end
            end
        end
    end
end

end # end of module


text:
export SuperpositionMeasure

"""
    struct SuperpositionMeasure{X,NT} <: AbstractMeasure
        components :: NT
    end

Superposition of measures is analogous to mixture distributions, but (because
measures need not be normalized) requires no scaling.

The superposition of two measures μ and ν can be more concisely written as μ + ν.
"""
struct SuperpositionMeasure{X,NT} <: AbstractMeasure
    components :: NT
end

SuperpositionMeasure(ms :: AbstractMeasure...) where {X} = SuperpositionMeasure{X,length(ms)}(ms)

# SuperpositionMeasure(m::NTuple{N, Measure{X}}) where {N,X} = SuperpositionMeasure(m...)

# Base.length(::SuperpositionMeasure{X,N}) where {X,N} = N

function Base.:+(μ::SuperpositionMeasure{X,N1}, ν::SuperpositionMeasure{X,N2}) where {X,N1,N2}
    components = (μ.components..., ν.components...)
    SuperpositionMeasure{X, N1+N2}(components)
end

function Base.:+(μ::AbstractMeasure, ν::SuperpositionMeasure{X,N}) where {X,N}
    components = (μ, ν.components...)
    SuperpositionMeasure{X,N+1}(components)
end

function Base.:+(μ::SuperpositionMeasure{X,N}, ν::AbstractMeasure) where {X,N}
    components = (μ.components..., ν)
    SuperpositionMeasure{X,N+1}(components)
end

function Base.:+(μ::AbstractMeasure, ν::AbstractMeasure) where {X}
    components = (μ, ν)
    SuperpositionMeasure{X,2}(components)
end

# TODO: Fix `rand` method (this one is wrong)
# function Base.rand(μ::SuperpositionMeasure{X,N}) where {X,N}
#     return rand(rand(μ.components))
# end

ext: jl


text:
using Gillespie
import Random: seed!

function F_sir(x,parms)
    (S,I,R) = x
    (beta,gamma) = parms
    infection = beta*S*I
    recovery = gamma*I
    [infection,recovery]
end

x0 = [999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/1000.0,0.01]
tf = 250.0
seed!(1234)

result = jensen(x0,F_sir,nu,parms,tf,35.0)

data = ssa_data(result)


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e16d5b2e-6361-11eb-3766-8b163c393ce7
using PlutoUI

# ╔═╡ 8c5bc2d0-6366-11eb-28bd-c335be197923
using ScikitLearn

# ╔═╡ 6c488620-635e-11eb-31ba-e11945ef7055
include("classifiers.jl")

# ╔═╡ 17d37150-6362-11eb-3516-43ed8dadcd21
names(datawithoutmissing)

# ╔═╡ be4fa400-6361-11eb-0237-c7f346aff057
md"""
# COVID Risk Analyzer App
**Intubed?:** $(@bind intubed Select(["0" => "Yes", "1" => "No"]))

**Have Pneumonia?:** $(@bind pneumonia Select(["0" => "Yes", "1" => "No"]))

**Age:** $(@bind age TextField(default="20"))

**Pregnant?:** $(@bind pregnant Select(["0" => "Yes", "1" => "No"]))

**Diabetic?:** $(@bind diabetic Select(["0" => "Yes", "1" => "No"]))

**Have Chronic Obstructive Pulmonary Disease?:** $(@bind copd Select(["0" => "Yes", "1" => "No"]))

**Have Asthma?:** $(@bind asthma Select(["0" => "Yes", "1" => "No"]))

**Immunosuppressed?:** $(@bind inmsupr Select(["0" => "Yes", "1" => "No"]))

**Hypertension?:** $(@bind hypertension Select(["0" => "Yes", "1" => "No"]))

**Have other diseases?:** $(@bind other_disease Select(["0" => "Yes", "1" => "No"]))

**Have Cardiovascular diseases?:** $(@bind cardiovascular Select(["0" => "Yes", "1" => "No"]))

**Obese?:** $(@bind obesity Select(["0" => "Yes", "1" => "No"]))

**Chronic Renal Disease?:** $(@bind renal_chronic Select(["0" => "Yes", "1" => "No"]))

**Smoker?:** $(@bind tobacco Select(["0" => "Yes", "1" => "No"]))

**Came in contact with COVID patient?:** $(@bind contact_other_covid Select(["0" => "Yes", "1" => "No"]))

"""

# ╔═╡ 98a0e810-6369-11eb-300e-bf80f1c9e8db
@bind go Button("Analyze")

# ╔═╡ 0f549420-6365-11eb-12ee-2d14ffe1de96
let
    go
    input =
        parse.(
            Int,
            [
                intubed,
                pneumonia,
                age,
                pregnant,
                diabetic,
                copd,
                asthma,
                inmsupr,
                hypertension,
                other_disease,
                cardiovascular,
                obesity,
                renal_chronic,
                tobacco,
                contact_other_covid,
            ],
        )
    md"""$(get_prediction(input))"""
end

# ╔═╡ 7ae6bcf2-6369-11eb-317c-772ddd065e0a

# ╔═╡ Cell order:
# ╠═e16d5b2e-6361-11eb-3766-8b163c393ce7
# ╠═8c5bc2d0-6366-11eb-28bd-c335be197923
# ╠═6c488620-635e-11eb-31ba-e11945ef7055
# ╠═17d37150-6362-11eb-3516-43ed8dadcd21
# ╟─be4fa400-6361-11eb-0237-c7f346aff057
# ╟─98a0e810-6369-11eb-300e-bf80f1c9e8db
# ╟─0f549420-6365-11eb-12ee-2d14ffe1de96
# ╠═7ae6bcf2-6369-11eb-317c-772ddd065e0a

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
mutable struct Framebuffer
    mPass::RenderPass
    mViews::Vector{ImageView}
    mViewHandles::Vector{vk.VkImageView}
    mCreateInfo::vk.VkFramebufferCreateInfo
    mOwnedViews

    mHandle::vk.VkFramebuffer

    function Framebuffer(pass::RenderPass, views::Vector{ImageView})
        this = new()
        this.mPass = pass
        this.mViews = views
        this.mViewHandles = Vector{vk.VkImageView}()
        this.mOwnedViews = nothing
        init(this)
        return this
    end

    function Framebuffer(pass::RenderPass, images::Vector{Image})
        this = new()
        this.mPass = pass
        this.mViews = Vector{ImageView}()
        this.mViewHandles = Vector{vk.VkImageView}()
        this.mOwnedViews = this.mViews
        for img in images
            range = vk.VkImageSubresourceRange(
                aspectsOf(img.mCreateInfo.format), # aspectMask::VkImageAspectFlags
                0, # baseMipLevel::UInt32
                1, # levelCount::UInt32
                0, # baseArrayLayer::UInt32
                img.mCreateInfo.arrayLayers # layerCount::UInt32
            )
            push!(this.mViews, img.createView(range))
        end
        init(this)
        return this
    end
end

@class Framebuffer

function destroy(this::Framebuffer)
    vk.vkDestroyFramebuffer(getVkDevice(this.mPass), this.mHandle, C_NULL)
    this.mHandle = C_NULL

    if this.mOwnedViews != nothing
        for v in this.mOwnedViews
            v.destroy()
        end
    end
end

function init(this::Framebuffer)
    empty!(this.mViewHandles)
    for imageView in this.mViews
        push!(this.mViewHandles, imageView.handle())
    end

    width = typemax(Int32)
    height = typemax(Int32)
    layers = typemax(Int32)
    for imageView in this.mViews
        image = imageView.getImage()
        if (image != nothing)
            width = min(width, image.getWidth())
            height = min(height, image.getHeight())
        end
        layers = min(layers, imageView.getLayers())
    end

    this.mCreateInfo = vk.VkFramebufferCreateInfo(
        vk.VK_STRUCTURE_TYPE_FRAMEBUFFER_CREATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        0, #flags::VkFramebufferCreateFlags
        this.mPass.handleRef()[], #renderPass::VkRenderPass
        length(this.mViewHandles), #attachmentCount::UInt32
        pointer(this.mViewHandles), #pAttachments::Ptr{VkImageView}
        width, #width::UInt32
        height, #height::UInt32
        layers, #layers::UInt32
    )

    framebuffer= Ref{vk.VkFramebuffer}()
    infoRef = Ref(this.mCreateInfo)
    GC.@preserve infoRef begin
        if (vk.vkCreateFramebuffer(this.mPass.getVkDevice(), infoRef, C_NULL, framebuffer) != vk.VK_SUCCESS)
            error("Failed to create framebuffer!")
        end
    end
    this.mHandle = framebuffer[]
end

function handle(this::Framebuffer)::vk.VkFramebuffer
    return this.mHandle
end

function pass(this::Framebuffer)::RenderPass
    return this.mPass
end

function getWidth(this::Framebuffer)::UInt32
    return this.mCreateInfo.width
end

function getHeight(this::Framebuffer)::UInt32
    return this.mCreateInfo.height
end


text:
# Own version of gensym:
#const symbol_number = [1]
const symbol_numbers = Dict{Symbol, Int}()

doc"""Return a new, unique symbol like _z10_"""

make_symbol() = make_symbol(:z)

function make_symbol(s::Symbol)

    i = get(symbol_numbers, s, 1)
    symbol_numbers[s] = i + 1

    Symbol("_", s, "_", i, "_")
end

doc"""
`insert_variables` takes a Julia `Expr`ession (i.e. an abstract syntax tree) and
recursively replaces operations like `a+b` by assignments
of the form `_z10_ = a+b`, where `_z10_` is a distinct symbol,
created using `make_symbol` (which is like `gensym`, but more readable).

Returns:

1. generated variable at head (top) of tree, which will contain the result of the whole tree;
2. sorted vector of (external) variables contained in the tree;
3. vector of intermediate variables (both introduced by the user, and generated);
4. code generated.

Usage: `IntervalConstraintProgramming.insert_variables(:(x^2 + y^2))`
"""
function insert_variables(ex)  # numbers are leaves
    ex, Symbol[], Symbol[], quote end
end

function insert_variables(ex::Symbol)  # symbols are leaves
    ex, [ex], Symbol[], quote end
end

function insert_variables(ex::Expr)

    if ex.head == :$   # process constants of form $a
        return :(esc($(ex.args[1]))), Symbol[], Symbol[], quote end

    elseif ex.head == :call
        process_call(ex)

    elseif ex.head == :(=)
        process_assignment(ex)

    elseif ex.head == :block
        process_block(ex)

    elseif ex.head == :tuple
        process_tuple(ex)

    elseif ex.head == :return
        process_return(ex)
    end
end

function process_tuple(ex)
    @show ex, ex.args

    new_code = quote end
    current_args = []  # the arguments in the current expression that will be added
    all_vars = Set{Symbol}()  # all variables contained in the sub-expressions
    top_vars = Symbol[]
    generated_variables = Symbol[]

    top_vars = Symbol[]

    for arg in ex.args
        top, contained_vars, generated, code = insert_variables(arg)

        union!(all_vars, contained_vars)
        append!(new_code.args, code.args)  # add previously-generated code
        append!(generated_variables, generated)

        push!(top_vars, top)
    end

    @show top_vars
    @show new_code
    #exit(1)

    return top_vars, sort(collect(all_vars)), generated_variables, new_code

end

function process_return(ex)
    top, contained_vars, generated, code = insert_variables(ex.args[1])
    append!(code.args, :(return $top))

    top, contained_vars, generated, code
end

function process_block(ex)

    new_code = quote end
    current_args = []  # the arguments in the current expression that will be added
    all_vars = Set{Symbol}()  # all variables contained in the sub-expressions
    generated_variables = Symbol[]

    local top

    for arg in ex.args[1:end]

        isa(arg, LineNumberNode) && continue
        (isa(arg, Expr) && arg.head == :line) && continue

        top, contained_vars, generated, code = insert_variables(arg)

        push!(current_args, top)
        union!(all_vars, contained_vars)
        append!(new_code.args, code.args)  # add previously-generated code
        append!(generated_variables, generated)
    end

    return top, sort(collect(all_vars)), generated_variables, new_code
end

function process_assignment(ex)
    # assumes ex is an assignment
    process_call(ex.args[2], ex.args[1])

end

function process_call(ex, new_var=nothing)
    # new_var is an optional variable name to assign the result of the call to
    # if none is given, then a new, unique variable name is created

    op = ex.args[1]

    if isa(op, Expr) && op.head == :line
         return quote end, Symbol[], Symbol[], quote end
     end

    # rewrite +(a,b,c) as +(a,+(b,c)):
    # TODO: Use @match here!

    if op in (:+, :*) && length(ex.args) > 3
        return insert_variables( :( ($op)($(ex.args[2]), ($op)($(ex.args[3:end]...) )) ))
    end

    new_code = quote end
    current_args = []  # the arguments in the current expression that will be added
    all_vars = Set{Symbol}()  # all variables contained in the sub-expressions
    generated_variables = Symbol[]

    for arg in ex.args[2:end]

        isa(arg, LineNumberNode) && continue

        top, contained_vars, generated, code = insert_variables(arg)

        push!(current_args, top)
        union!(all_vars, contained_vars)
        append!(new_code.args, code.args)  # add previously-generated code
        append!(generated_variables, generated)
    end

    top_level_code = quote end

    #@show op

    if op ∈ keys(rev_ops)  # standard operator
        if new_var == nothing
            new_var = make_symbol()
        end

        push!(generated_variables, new_var)

        top_level_code = :($(new_var) = ($op)($(current_args...)))  # new top-level code

    else  # assume user-defined function

        function_name = :($(op).forward)  # need esc?

        func_args = registered_functions[op]
        @show func_args

        new_generated_vars = Symbol[]
        for i in func_args.generated
            push!(new_generated_vars, make_symbol())
        end

        append!(generated_variables, new_generated_vars)
        new_var = new_generated_vars[end]

        new_generated_vars = Expr(:tuple, new_generated_vars...)

        top_level_code = :($(new_generated_vars) = $(function_name)($(current_args...)))  # new top-level code

    end

    push!(new_code.args, top_level_code)

    return new_var, sort(collect(all_vars)), generated_variables, new_code

end

function constraint(root_var, constraint_interval)
    # if constraint == Interval(-∞, ∞)
    #     constraint_code = :($(root_var) = $(root_var) ∩ _A_)
    #     # push!(all_vars, :_A_)
    #
    # else
        constraint_code = :( $(root_var) = $(root_var) ∩ $(constraint_interval) )
    # end

    return constraint_code

    # new_code = quote end
    # push!(new_code.args, constraint_code)
end

function forward_pass(ex::Expr)
    root, all_vars, generated, code = insert_variables(ex)
    forward_pass(root, all_vars, generated, code)
end

function forward_pass(root, all_vars, generated, code)
    make_function(all_vars, generated, code)
end

function backward_pass(ex::Expr) #, constraint::Interval)
    root, all_vars, generated, code = insert_variables(ex)
    backward_pass(root, all_vars, generated, code)
end

doc"""`backward_pass` replaces e.g. `z = a + b` with
the corresponding reverse-mode function, `(z, a, b) = plusRev(z, a, b)`
"""

function backward_pass(root_var, all_vars, generated, code)

    new_code = quote end

    for line in reverse(code.args)  # run backwards

        line.head == :line  && continue  # ignore line number nodes

        (var, op, args) = @match line begin
            (var_ = op_(args__))  => (var, op, args)
        end

        #@show "hello", op

        if (@capture(op, f_.forward))  # user-defined forward
            rev_op = :($f.backward)
            #println("Hello there")
            @show var, op, args

            return_args = [args...; var.args...]
            rev_code = :($(rev_op)($(return_args...)))

            #return_args =

        else
            return_args = [var, args...]
            rev_op = rev_ops[op]  # find reverse operation

            rev_code = :($(rev_op)($(return_args...)))

            # delete non-symbols in return args:
            for (i, arg) in enumerate(return_args)
                if !(isa(arg, Symbol))
                    return_args[i] = :_
                end
            end

         end

        return_tuple = Expr(:tuple, return_args...)  # make tuple out of array
        # or: :($(return_args...),)

        new_line = :($(return_tuple) = $(rev_code))
        push!(new_code.args, new_line)

    end

    sort!(all_vars)

    make_function(vcat(all_vars, generated), all_vars, new_code)
end

doc"""
`forward_backward` takes in an expression like `x^2 + y^2` and outputs
code for the forward-backward contractor

TODO: Add intersections in forward direction
"""
function forward_backward(ex::Expr, constraint::Interval=entireinterval())

    new_ex = copy(ex)

    # Step 1: Forward pass using insert_variables

    root_var, all_vars, generated, code = insert_variables(new_ex)

    # Step 2: Add constraint code:

    local constraint_code

    if constraint == Interval(-∞, ∞)
        constraint_code = :($(root_var) = $(root_var) ∩ _A_)
        push!(all_vars, :_A_)

    else
        constraint_code = :($(root_var) = $(root_var) ∩ $constraint)
    end

    new_code = copy(code)
    push!(new_code.args, constraint_code)

    # Step 3: Backwards pass
    # replace e.g. z = a + b with reverse mode function plusRev(z, a, b)

    for line in reverse(code.args)  # run backwards

        if line.head == :line  # line number node
            continue
        end

        (var, op, args) = @match line begin
            (var_ = op_(args__))  => (var, op, args)
        end

        new_args = []
        push!(new_args, var)
        append!(new_args, args)

        rev_op = rev_ops[op]  # find the reverse operation

        rev_code = :($(rev_op)($(new_args...)))

        return_args = copy(new_args)

        # delete non-symbols in return args:
        for (i, arg) in enumerate(return_args)
            if !(isa(arg, Symbol))
                return_args[i] = :_
            end
        end

        return_tuple = Expr(:tuple, return_args...)  # make tuple out of array
        # or: :($(return_args...),)

        new_line = :($(return_tuple) = $(rev_code))
        push!(new_code.args, new_line)
    end

    sort(all_vars), new_code
end

function make_function(all_vars, code)

    vars = Expr(:tuple, all_vars...)  # make a tuple of the variables

    if all_vars[1] == :_A_
        vars2 = Expr(:tuple, (all_vars[2:end])...)  # miss out _A_
        push!(code.args, :(return $(vars2)))
    else
        push!(code.args, :(return $(vars)))
    end

    # @show code

    function_code = :( $(vars) -> $(code) )

    function_code
end

doc"""
Generate code for an anonymous function with given
input arguments, output arguments, and code block.
"""
function make_function(input_args, output_args, code)

    input = Expr(:tuple, input_args...)  # make a tuple of the variables
    output = Expr(:tuple, output_args...)  # make a tuple of the variables

    new_code = copy(code)
    push!(new_code.args, :(return $output))

    return :( $input -> $new_code )
end

doc"""`parse_comparison` parses comparisons like `x >= 10`
into the corresponding interval, expressed as `x ∈ [10,∞]`

Returns the expression and the constraint interval

TODO: Allow something like [3,4]' for the complement of [3,4]'"""

function parse_comparison(ex)
    expr, limits =
    @match ex begin
       ((a_ <= b_) | (a_ < b_))   => (a, (-∞, b))
       ((a_ >= b_) | (a_ > b_))   => (a, (b, ∞))

       ((a_ == b_) | (a_ = b_))   => (a, (b, b))

       ((a_ <= b_ <= c_)
        | (a_ < b_ < c_)
        | (a_ <= b_ < c)
        | (a_ < b_ <= c))         => (b, (a, c))

       ((a_ >= b_ >= c_)
       | (a_ > b_ > c_)
       | (a_ >= b_ > c_)
       | (a_ > b_ >= c))          => (b, (c, a))

       ((a_ ∈ [b_, c_])
       | (a_ in [b_, c_])
       | (a_ ∈ b_ .. c_)
       | (a_ in b_ .. c_))        => (a, (b, c))

       _                          => (ex, (-∞, ∞))

   end

   a, b = limits

   return (expr, a..b)   # expr ∈ [a,b]

end

type Contractor
    variables::Vector{Symbol}
    constraint_expression::Expr
    contractor::Function
    code::Expr
end

function Contractor(ex::Expr)
    expr, constraint_interval = parse_comparison(ex)

    vars, code = forward_backward(expr, constraint_interval)

    fn = eval(make_function(vars, code))

    Contractor(vars, expr, fn, code)
end

# new call syntax to define a "functor" (object that behaves like a function)
@compat (C::Contractor)(x...) = C.contractor(x...)

function Base.show(io::IO, C::Contractor)
    println(io, "Contractor:")
    println(io, "  - variables: $(C.variables)")
    print(io, "  - constraint: $(C.constraint_expression)")
end

doc"""Usage:
```
C = @contractor(x^2 + y^2 <= 1)
x = y = @interval(0.5, 1.5)
C(x, y)

`@contractor` makes a function that takes as arguments the variables contained in the expression, in lexicographic order
```

TODO: Hygiene for global variables, or pass in parameters
"""

function contractor(ex)
    expr, constraint = parse_comparison(ex)
    @show expr, constraint

    all_vars, code = forward_backward(expr, constraint)
    @show all_vars, code

    make_function(all_vars, code)
end

macro contractor(ex)
    ex = Meta.quot(ex)
    :(Contractor($ex))
end

show_code(c::Contractor) = c.code


text:
include("src/quick.jl")
include("src/merge.jl")
include("src/selection.jl")
include("src/insertion.jl")
include("src/bubble.jl")
include("src/counting.jl")

function main(args)
    if length(args) ≠ 6
        error("usage: julia driver.jl <algo> <nrep> <start> <stop> <inc> <filename>\nfound $(args)")
    end

    algo::String     =            args[1]
    nrep::Int        = parse(Int, args[2])
    start::Int       = parse(Int, args[3])
    stop::Int        = parse(Int, args[4])
    inc::Int         = parse(Int, args[5])
    filename::String =            args[6]

    dtime::Float64      = 0.0
    dtime_best::Float64 = 1000.0

    @show algo
    @show nrep
    @show start
    @show stop
    @show inc
    @show filename

    buffer = "size, time\n"

    if algo == "quick"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = randn(len)
                dtime = time()
                quicksort!(a, 1, len);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end
    if algo == "merge"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = randn(len)
                dtime = time()
                mergesort!(a, 1, len);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end
    if algo == "insertion"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = randn(len)
                dtime = time()
                insertionsort!(a);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end
    if algo == "selection"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = randn(len)
                dtime = time()
                selectionsort!(a);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end
    if algo == "bubble"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = randn(len)
                dtime = time()
                bubblesort!(a);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end
    if algo == "counting"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = rand(0:1_000_000, len)
                dtime = time()
                countingsort(a);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end
    if algo == "inbuilt"
        for len ∈ start:inc:stop
            for rep ∈ 1:nrep
                a = randn(len)
                dtime = time()
                sort!(a);
                dtime = time() - dtime
                if rep == 1 dtime_best = dtime
                else dtime_best = ( dtime < dtime_best ? dtime : dtime_best )
                end
            end
            buffer *= "$len, $(dtime_best*1.0e09)\n"
        end
    end

    open(filename, "w") do io
        println("Writing Results...")
        write(io, buffer)
    end
end

main(ARGS)


text:
using ShortStrings, SortingAlgorithms, BenchmarkTools, SortingLab, DataFrames, CSV
N = 100_000_000; K=100;
using Plots
using StatPlots, RCall

function string_sort_perf(N, K)
    tic()
    ss ="id".*dec.(1:N÷K,10);
    srand(1);
    ss1 = rand(ss, N);
    svec = ShortString.(ss1);

    rs = @benchmark sort($svec, by = x->x.size_content, alg = RadixSort) samples = 5

    rs1 = @benchmark SortingLab.radixsort($ss1) samples = 5

    @rput ss1
    r_times = R"""
    memory.limit(2^31-1)
    replicate(max(5, $(length(rs.times))), system.time(sort(ss1, method="radix"))[3])
    """

    df = DataFrame(xlabel = ["Julia - Strings", "Julia - ShortStrings", "R"], timing = mean.([rs1.times/1e9, rs.times./1e9, r_times]))
    CSV.write("String sort perf $(N÷1_000_000)m.csv", df)
    gb = bar(
        df[:xlabel],
        df[:timing],
        title = "String sort perf $(N÷1_000_000)m",
        label="seconds"
        )
    savefig(gb, "String sort perf $(N÷1_000_000)m.png")
    toc()
end

string_sort_perf(100, 100) # warm up

@time string_sort_perf.((1:10).*10^8, 100)


text:
using Documenter
using CKANQuery

makedocs(
    sitename = "CKANQuery",
    format = Documenter.HTML(),
    modules = [CKANQuery]
)

deploydocs(
    repo = "github.com/VLucet/CKANQuery.jl"
)


text:
function raytrace_image(scene::Scene, cam::Camera{T}, image_height::Integer, samples_per_pixel::Integer, max_bounces_per_ray::Integer, max_steps_per_bounce::Integer, distance_tolerance::T, rng::MersenneTwister) where {N<:Integer, T<:Real}
    aspect_ratio_cam = aspect_ratio(cam)
    image_width = Int64(floor(aspect_ratio_cam * image_height))
    out_image = zeros(RGB{T}, image_height, image_width)
    for i in 1:image_height, j in 1:image_width
        pixel_color = zero(Vec{3, T})
        for rands in 1:samples_per_pixel
            s = T((j - 1 + rand(rng))/(image_width - 1))
            t = T((i - 1 + rand(rng))/(image_height - 1))
            ray = generate_ray(cam, s, t)
            ray_color = raymarch(scene, ray, max_bounces_per_ray, max_steps_per_bounce, distance_tolerance, rng::MersenneTwister)
            pixel_color += ray_color ./ samples_per_pixel
        end
        out_image[i, j] = RGB{T}(pixel_color...)
    end
    out_image
end

function raytrace_image_cuda(scene::Scene, cam::Camera{T}, image_height::Integer, samples_per_pixel::Integer, max_bounces_per_ray::Integer, max_steps_per_bounce::Integer, distance_tolerance::T, rng::MersenneTwister) where {N<:Integer, T<:Real}
    aspect_ratio_cam = aspect_ratio(cam)
    image_width = Int64(floor(aspect_ratio_cam * image_height))
    out_image = zeros(RGB{T}, image_height, image_width)
    for i in 1:image_height, j in 1:image_width
        pixel_color = zero(Vec{3, T})
        for rands in 1:samples_per_pixel
            s = T((j - 1 + rand(rng))/(image_width - 1))
            t = T((i - 1 + rand(rng))/(image_height - 1))
            ray = generate_ray(cam, s, t)
            ray_color = raymarch(scene, ray, max_bounces_per_ray, max_steps_per_bounce, distance_tolerance, rng::MersenneTwister)
            pixel_color += ray_color ./ samples_per_pixel
        end
        out_image[i, j] = RGB{T}(pixel_color...)
    end
    out_image
end

#=

CUDA.@sync begin
        #@cuda threads=numthreads blocks=numblocks mandelbrotandregiongpu!(cudaimage,
            #centerx, xstart, xrangeextent, centery, ystart, yrangeextent,
            #numiters, height, width, cosrot, sinrot, a, b, c, d)
        @cuda threads=numthreads blocks=numblocks mandelbrotandregiongpurational!(cudaimage,
            centerx, xstart, xrangeextent, centery, ystart, yrangeextent,
            numiters, height, width, cosrot, sinrot, conformal)
    end

    CUDA.copyto!(outimage[], cudaimage

    function mandelbrotandregiongpurational!(escape_color, centerx::T, xstart::T, xrangeextent::T, centery::T, ystart::T, yrangeextent::T,
        numiters, height, width, cosrot::T, sinrot::T, rational::ComplexRational{T, N}) where {T <: Real, N}
    indexx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    indexy = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stridex = blockDim().x * gridDim().x
    stridey = blockDim().y * gridDim().y
    #@cuprintln("thread $indexx, $indexy, stride $stridex, $stridey, weight $width $height")
    for i in indexx:stridex:width, j in indexy:stridey:height
        prex_ij = (i - 1) / (width - 1) * xrangeextent + xstart
        prey_ij = (j - 1) / (height - 1) * yrangeextent + ystart
        z_ij = Complex{T}(cosrot * prex_ij - sinrot * prey_ij + centerx, sinrot * prex_ij + cosrot * prey_ij + centery)
        z_ij_rational = transform(rational, z_ij)
        escape_ij = mandelbrot(z_ij_rational, numiters) / numiters
        escape_ij_f32 = Float32(escape_ij)
        @inbounds escape_color[i, j] = ColorTypes.RGBA{Float32}(escape_ij_f32, escape_ij_f32, escape_ij_f32)
    end
    return nothing
end
=#

function raymarch(scene::Scene, ray::Ray{T}, max_bounces_left::Integer, max_steps_per_bounce::Integer, distance_tolerance::T, rng::MersenneTwister) where {T<:Real}
    if max_bounces_left > 0
        t = T(0)
        for i in 1:max_steps_per_bounce
            current_point = origin(ray) + t*direction(ray)
            current_min_dist, current_min_index = min_dist_index(scene, current_point)
            if current_min_dist < t * distance_tolerance
                min_shape = shapes(scene)[current_min_index]
                outward_normal = normal(sdf(min_shape), current_point)
                hitrecord = HitRecord(ray, outward_normal, t)
                color, ray_out, absorbed = scatter(mat(min_shape), ray, hitrecord, rng)
                if absorbed
                    return color
                else
                    return color .* raymarch(scene, ray_out, max_bounces_left - 1, max_steps_per_bounce, distance_tolerance, rng)
                end
            else
                # Assumes the ray direction is normalized
                t += current_min_dist
            end
        end
    end
    # max bounces <= 0 or ran out of march steps
    zero(Vec{3, T})
end

saveimage(savename::AbstractString, image) = save("./images/$(savename)", map(clamp01nan, image))


text:
using DashBootstrapComponents, DashHtmlComponents

progress = html_div([
    dbc_progress(value=80, id="animated-progress", animated=false, striped=true),
    dbc_button(
        "Toggle animation",
        id="animation-toggle",
        className="mt-3",
        n_clicks=0,
    ),
]);

callback!(
    app,
    Output("animated-progress", "animated"),
    Input("animation-toggle", "n_clicks"),
    State("animated-progress", "animated"),
) do n, animated
    return n > 0 ? animated == 0 : animated
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using EllipticalSliceSampling

using Distances
using Distributions
using FillArrays

using Distributed
using LinearAlgebra
using Random
using Statistics
using Test

Random.seed!(0)

# add additional processes
addprocs(2)
@everywhere begin
    using EllipticalSliceSampling
    using Distributions

    using LinearAlgebra
end

@testset "EllipticalSliceSampling" begin
    println("Simple tests")
    @time @testset "Simple tests" begin
        include("simple.jl")
    end

    println("GP regression tests")
    @time @testset "GP regression tests" begin
        include("regression.jl")
    end
end


text:
module SMR

include("./SMRTypes.jl")

using .SMRTypes

export read_wavemark_channel, read_continuous_channel, read_event_channel,
       read_marker_channel, read_channel_info, get_channel_type, channel_string,
       get_read_function

include("./helpers.jl")

# ============================================================================ #
# global constant path to libsmr shared object file
const LIBSMR = @libpath("libsmr", @__DIR__, "..", "lib")

# ============================================================================ #
"""
`idx = get_channel_index(ifile, label)`
In:
    ifile - path to .smr file
    label - channel label

Out:
    idx - integer index of channel, -1 indicates failure
"""
function get_channel_index(ifile::String, label::String)::Int64
    return ccall((:channel_label_path_to_index, LIBSMR), Cint, (Cstring,
        Cstring), ifile, label)
end
# ============================================================================ #
"""
`wmrk = read_wavemark_channel(ifile::String, idx::Integer)` *OR*\n
`wmrk = read_wavemark_channel(ifile::String, label::String)`
### Input:
* ifile - the path to a .smr file
* idx - an integer channel index
* label - the channel label as a String

### Output:
* wmrk - a SMRWMrkChannel type with fields:\n
            timestamps: Nx1 Vector{Float64} of spike timestamps in seconds
            wavemarks: MxN Array{UInt8,2} of spike waveforms
            markers: 4xN Array{UInt8,2} of marker codes
"""
function read_wavemark_channel(ifile::String, idx::Integer)
    @calllib(ifile, idx, "wavemark", SMRWMrkChannel)
end
function read_wavemark_channel(ifile::String, label::String)
    return read_wavemark_channel(ifile, get_channel_index(ifile, label))
end
# ============================================================================ #
"""
`cont = read_continuous_channel(ifile::String, idx::Integer)` *OR*\n
`cont = read_continuous_channel(ifile::String, label::String)`
### Input:
 * see `read_wavemark_channel`

### Output:
* cont - a SMRContChannel type with fields:\n
            data: Nx1 Vector{Int16} of voltage values
            sampling_rate: channel sampling rate in Hz
"""
function read_continuous_channel(ifile::String, idx::Integer)
    @calllib(ifile, idx, "continuous", SMRContChannel)
end
function read_continuous_channel(ifile::String, label::String)
    return read_continuous_channel(ifile, get_channel_index(ifile, label))
end
# ============================================================================ #
"""
`evt = read_event_channel(ifile::String, idx::Integer)` *OR*\n
`evt = read_event_channel(ifile::String, label::String)`
### Input:
 * see `read_wavemark_channel`

### Output:
* evt - a SMREventChannel type with fields:\n
            data: Nx1 Vector{Float64} of event timestamps
"""
function read_event_channel(ifile::String, idx::Integer)
    @calllib(ifile, idx, "event", SMREventChannel)
end
function read_event_channel(ifile::String, label::String)
    return read_event_channel(ifile, get_channel_index(ifile, label))
end
# ============================================================================ #
"""
`mrk = read_marker_channel(ifile::String, idx::Integer)` *OR*\n
`mrk = read_marker_channel(ifile::String, label::String)`
### Input:
 * see `read_wavemark_channel`

### Output:
* mrk - a SMRMarkerChannel type with fields:\n
            timestamps: Nx1 Vector{Float64} of marker timestamps in seconds
            markers: 4xN Array{UInt8,2} of marker codes
            text: an Nx1 Vector{String} of text info for each marker
"""
function read_marker_channel(ifile::String, idx::Integer)
    @calllib(ifile, idx, "marker", SMRMarkerChannel)
end
function read_marker_channel(ifile::String, label::String)
    return read_marker_channel(ifile, get_channel_index(ifile, label))
end
# ============================================================================ #
"""
`ifo = read_channel_info(ifile)`
### Input:
* ifile - the path to a .smr file

### Output:
* ifo - a Vector{SMRChannelInfo} with one element per channel, try `display(ifo)` to see channel infomration
"""
function read_channel_info(ifile::String)

    ptr = ccall((:read_channel_array, LIBSMR), Ptr{cSMRChannelInfoArray}, (Cstring,), ifile)

    if ptr != C_NULL
        err = false
        ifo_array = unsafe_load(ptr)
        if ifo_array.ifo != C_NULL
            ifo = Array{SMRChannelInfo}(undef, ifo_array.length)

            for k = 1:ifo_array.length
                tmp = unsafe_load(ifo_array.ifo ,k)
                if tmp != C_NULL
                    ifo[k] = SMRChannelInfo(unsafe_load(tmp))
                else
                    print("WARNING: Failed to read array element " * string(k))
                end
            end
        else
            err = true
        end

        ccall((:free_channel_info_array, LIBSMR), Cvoid, (Ptr{cSMRChannelInfoArray},), ptr)

        if err
            error("failed to load array of channel info structs")
        end
    else
        error("call to read_channel_array failed")
    end

    return ifo
end
# ============================================================================ #
function get_channel_type(ifile::String, label::String)
    return get_channel_type(read_channel_info(ifile), label)
end
# ============================================================================ #
function get_channel_type(ifo::Array{T}, label::String) where T<:SMRChannelInfo
    tp = 0
    label = lowercase(label)

    for x in ifo
        if label == lowercase(x.title)
            tp = x.kind
            break
        end
    end
    return tp
end
# ============================================================================ #
function get_read_function(ifile::String, label::String)
    typ = get_channel_type(ifile, label)
    if typ == 1
        f = read_continuous_channel
    elseif typ == 2
        f = read_event_channel
    elseif typ == 3
        f = read_event_channel
    elseif typ == 4
        f = read_event_channel
    elseif typ == 5
        f = read_marker_channel
    elseif typ == 6
        f = read_wavemark_channel
    elseif typ == 7
        f = read_marker_channel
    elseif typ == 8
        f = read_marker_channel
    elseif typ == 9
        error("real wave channels are not yet supported!")
    else
        error("Channel $(label) cannot be located")
    end
end
# ============================================================================ #
function test_all(ifile::String="../b1_con_006.smr")

    ifo = read_channel_info(ifile)

    for chan in ifo
        str = channel_string(chan.kind)
        f = get_read_function(ifile, chan.title)
        print("Reading ", chan, "... ")
        x = f(ifile, chan.title)
        println("success")
    end

    return ifo
end
# ============================================================================ #
end #END MODULE

max_stars_repo_name: scottiealexander/libsmr

max_issues_repo_name: scottiealexander/libsmr

max_forks_repo_name: scottiealexander/libsmr


text:
using Revise
using PyPlot, DelimitedFiles, Random, Statistics, transD_GP

Random.seed!(23)

zfixed = [-1e5]
ρfixed = [1e12]
nmax = 100

zstart = 0.0
extendfrac, dz = 1.03, 1.5
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=65, showplot=true)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)

## geometry parameters for tempest
zTx = -120
zRx = -82
x_rx = -113.0
y_rx = 0.
rx_roll = 0.
rx_pitch = 2.
rx_yaw = 0.
tx_roll = 0.
tx_pitch = 0.
tx_yaw = 0.
# electronics and stuff
include("electronics_halt.jl")
## fill in detail in ohm-m
## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<50)] .= 20.
ρ[(z.>=50) .&(z.<80)] .= 1
ρ[(z.>=80) .&(z.<100)] .= 20
ρ[(z.>=100) .&(z.<200)] .= 50
ρ[(z.>=200) .&(z.<250)] .= 80
ρ[(z.>=250)]            .= 150
# add jitter to model in log10 domain
Random.seed!(11)
ρ = 10 .^(0.1*randn(length(ρ)) + log10.(ρ))
## create total field operator (required for nuisance inversion)
tempest = transD_GP.TEMPEST1DInversion.Bfield(
    zTx = zTx, zRx = zRx, x_rx = x_rx, y_rx = y_rx,
    rx_roll = rx_roll, rx_pitch = rx_pitch, rx_yaw = rx_yaw,
    tx_roll = tx_roll, tx_pitch = tx_pitch, tx_yaw = tx_yaw,
	ramp = ramp, times = times,
	z=z,
	ρ=ρ,
	addprimary = true #this ensures that the geometry update actually changes everything that needs to be
)
# plot before adding noise
transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,z,ρ)
## compute noisy data to invert
# remember noise in electronics_halt.jl are in fT!!
transD_GP.TEMPEST1DInversion.set_noisy_data!(tempest, z, ρ,
                        noisefracx = 0.02, noisefracz = 0.02,
                        halt_X = Hx_add_noise*1e-15, halt_Z = Hz_add_noise*1e-15)
# but model with a coarser grid
extendfrac, dz = 1.06, 1.15
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=50, showplot=true)
zgrid, ρgrid, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
tempest.z, tempest.ρ = zgrid, copy(ρgrid)
# only primary field stuff if you want for GA-AEM
# Hxp, Hyp, Hzp = transD_GP.TEMPEST1DInversion.returnprimary!(tempest)
# Xnoisy, Znoisy = tempest.dataHx - Hxp, tempest.dataHz - Hzp # raw SI units!!


text:
struct EnumerableTake{T,S} <: Enumerable
    source::S
    n::Int
end

function take(source::Enumerable, n::Integer)
    T = eltype(source)
    S = typeof(source)
    return EnumerableTake{T,S}(source, Int(n))
end

Base.IteratorSize(::Type{EnumerableTake{T,S}}) where {T,S} = haslength(S)

Base.eltype(::Type{EnumerableTake{T,S}}) where {T,S} = T

Base.length(iter::EnumerableTake{T,S}) where {T,S} = min(length(iter.source),iter.n)

function Base.iterate(iter::EnumerableTake{T,S}) where {T,S}
    ret = iterate(iter.source)

    if ret===nothing
        return nothing
    elseif iter.n==0
        return nothing
    else
        return ret[1], (ret[2],1)
    end
end

function Base.iterate(iter::EnumerableTake{T,S}, state) where {T,S}
    if state[2]==iter.n
        return nothing
    else
        ret = iterate(iter.source, state[1])

        if ret===nothing
            return nothing
        else
            return ret[1], (ret[2], state[2]+1)
        end
    end
end


text:
"""
    compute_scc(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta=1.45)

Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiDICE2013 model.
If no model is provided, the default model from MimiDICE2013.get_model() is used.
Constant discounting is used from the specified pure rate of time preference `prtp`.
"""
function compute_scc(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta=1.45)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing
    !(year in model_years[1]:5:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):5:$last_year.") : nothing

    mm = get_marginal_model(m; year=year)

    return _compute_scc(mm, year=year, last_year=last_year, prtp=prtp, eta=eta)
end

"""
    compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta=1.45)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiDICE2013 model.
If no model is provided, the default model from MimiDICE2013.get_model() is used.
Constant discounting is used from the specified pure rate of time preference `prtp`.
"""
function compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta=1.45)
    year === nothing ? error("Must specify an emission year. Try `compute_scc_mm(m, year=2020)`.") : nothing
    !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing
    !(year in model_years[1]:5:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):5:$last_year.") : nothing

    mm = get_marginal_model(m; year=year)
    scc = _compute_scc(mm; year=year, last_year=last_year, prtp=prtp, eta=eta)

    return (scc = scc, mm = mm)
end

# helper function for computing SCC from a MarginalModel, not to be exported or advertised to users
function _compute_scc(mm::MarginalModel; year::Int, last_year::Int, prtp::Float64, eta::Float64)
    ntimesteps = findfirst(isequal(last_year), model_years)     # Will run through the timestep of the specified last_year
    run(mm, ntimesteps=ntimesteps)

    # Go from trillion$ to $; multiply by -1 so that damages are positive; pulse
    # was already converted to CO2 so we don't need to multiply by 12/44
    marginal_damages = -1 * mm[:neteconomy, :C][1:ntimesteps] * 1e12

    cpc = mm.base[:neteconomy, :CPC]

    year_index = findfirst(isequal(year), model_years)

    df = [zeros(year_index - 1)..., ((cpc[year_index] / cpc[i])^eta * 1 / (1 + prtp)^(t - year) for (i, t) in enumerate(model_years) if year <= t <= last_year)...]
    scc = sum(df .* marginal_damages * 5)  # currently implemented as a 5year step function; so each timestep of discounted marginal damages is multiplied by 5
    return scc
end

"""
    get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing)

Creates a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
If no Model m is provided, the default model from MimiDICE2013.get_model() is used as the base model.
"""
function get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2015)`.") : nothing
    !(year in model_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    # note here that the pulse size will be used as the `delta` parameter for
    # the `MarginalModel` and thus allow computation of the SCC to return units of
    # dollars per ton, as long as `pulse_size` is in tons
    mm = create_marginal_model(m, 5 * 1e9)    # 1 GtCO2 per year for 5 years, so 5 * 10^9
    add_marginal_emissions!(mm.modified, year)

    return mm
end

"""
    add_marginal_emissions!(m::Model, year::Int)

Adds a marginal emission component to year m which adds 1Gt of additional CO2 emissions per year for five years starting in the specified `year`.
"""
function add_marginal_emissions!(m::Model, year::Int)
    add_comp!(m, Mimi.adder, :marginalemission, before=:co2cycle)

    time = Mimi.dimension(m, :time)
    addem = zeros(length(time))
    addem[time[year]] = 1.0     # 1 GtCO2 per year for five years

    set_param!(m, :marginalemission, :add, addem)
    connect_param!(m, :marginalemission, :input, :emissions, :E)
    connect_param!(m, :co2cycle, :E, :marginalemission, :output)
end

# Old available marginal model function
function getmarginal_dice_models(;emissionyear=2010)

    DICE = constructdice()
    run(DICE)

    mm = MarginalModel(DICE)
    m1 = mm.base
    m2 = mm.modified

    add_comp!(m2, adder, :marginalemission, before=:co2cycle)

    time = dimension(m1, :time)
    addem = zeros(length(time))
    addem[time[emissionyear]] = 1.0

    set_param!(m2, :marginalemission, :add, addem)
    connect_param!(m2, :marginalemission, :input, :emissions, :E)
    connect_param!(m2, :co2cycle, :E, :marginalemission, :output)

    run(m1)
    run(m2)

    return m1, m2
end


text:
# Wrapper types and functions for XML interaction.

const Node = EzXML.Node

"""
    $(DocStringExtensions.SIGNATURES)

Get the root node of an XML document.
"""
xml_hasroot(xml) = EzXML.hasroot(xml)

"""
    $(DocStringExtensions.SIGNATURES)

Parse a string into an XML document.
"""
parse_xml(xmlstring) = EzXML.parsexml(xmlstring)


text:
using DynamicEnergyBudgets, Unitful, Test

using DynamicEnergyBudgets: rate_formula, half_saturation

@testset "Half saturation curve" begin
    mx = 10.0
    half = 5.0
    x = 5.0
    @test half_saturation(mx, half, x) == 5.0

    # Maximum
    x = 100000000000000000000000000.0
    @test half_saturation(mx, half, x) == 10.0

    # Minimum
    x = 0.0
    @test half_saturation(mx, half, x) == 0.0
end

@testset "Rate CNE" begin
    j_E_mai    = 0.003
    y_E_CH_NO  = 1.0/1.5
    y_E_EN     = 1.0/0.7
    y_V_E      = 0.7
    kap_soma   = 0.6
    r          = 0.001
    reserve    = (1.0, 1.0, 1.0)
    turnover   = ((0.2, 0.2, 0.2) .* 0.05)
    su = ParallelComplementarySU()

    rf = rate_formula(r, su, reserve, turnover, j_E_mai, y_E_CH_NO, y_E_EN, y_V_E, kap_soma)

    @test rf == 0.002874195250659631

    @testset "test rate with units" begin
        j_E_mai = 0.003u"mol/mol*d^-1"
        turnover = ((0.2u"mol/mol*d^-1", 0.2u"mol/mol*d^-1", 0.2u"mol/mol*d^-1") .* 0.05)
        mass = (1.0u"mol/mol", 1.0u"mol/mol", 1.0u"mol/mol")
        r = 0.001u"mol/mol*d^-1"
        rf = rate_formula(r, su, reserve, turnover, j_E_mai, y_E_CH_NO, y_E_EN, y_V_E, kap_soma)
        @test rf == 0.002874195250659631u"mol/mol*d^-1"
    end
end

@testset "Rate CN" begin
    j_E_mai    = 0.003
    y_E_CH_NO  = 1.0/1.5
    y_E_EN     = 1.0/0.7
    y_V_E      = 0.7
    kap_soma   = 0.6
    r          = 0.001
    reserve    = (10.0, 1.0)
    turnover   = ((0.2, 0.2) .* 0.05)
    su = ParallelComplementarySU()

    rf = rate_formula(r, su, reserve, turnover, j_E_mai, y_E_CH_NO, y_E_EN, y_V_E, kap_soma)

    # @test rf == 0.002874195250659631

    @testset "test rate with units" begin
        j_E_mai = 0.003u"mol/mol*d^-1"
        turnover = ((0.2u"mol/mol*d^-1", 0.2u"mol/mol*d^-1") .* 0.05)
        mass = (1.0u"mol/mol", 1.0u"mol/mol")
        r = 0.001u"mol/mol*d^-1"
        rf = rate_formula(r, su, reserve, turnover, j_E_mai, y_E_CH_NO, y_E_EN, y_V_E, kap_soma)
        #@test rf == 0.002874195250659631u"mol/mol*d^-1"
    end
end

nothing


text:
using LinearAlgebra
using Einsum
using Test
using NLsolve
using Plots
plotly()

# Tensors
I_ = Matrix(1.0I,3,3) # Second order identity tensor
II = zeros(3,3,3,3)
@einsum II[i,j,k,l] = 0.5*(I_[i,k]*I_[j,l] + I_[i,l]*I_[j,k]) # Fourth order symmetric identity tensor
IxI = zeros(3,3,3,3)
@einsum IxI[i,j,k,l] = I_[i,j]*I_[k,l] # "Trace" tensor
P = II - 1/3*IxI # Deviatoric projection tensor

# Functions
function double_contraction(x::AbstractArray{<:Number,2},y::AbstractArray{<:Number,2})
    return sum(x.*y)
end

function double_contraction(x::AbstractArray{<:Number,4},y::AbstractArray{<:Number,2})
    retval = zeros(3,3)
    @einsum retval[i,j] = x[i,j,k,l]*y[k,l]
    return retval
end

A = rand(3,3) # Create Random second order tensor
A += A' # Symmetrize it

@test isapprox(double_contraction(II, A), A)
@test isapprox(double_contraction(IxI, A), I_*tr(A))

function deviator(x::AbstractArray{<:Number,2})
    s = zeros(3,3)
    @einsum s[i,j] = P[i,j,k,l]*x[k,l]
    return s
end

function von_mises_stress(stress::AbstractArray{<:Number,2})
    s = deviator(stress)
    return sqrt(3/2*double_contraction(s,s))
end

S = [100 0 0; 0 0 0; 0 0 0]
@test isapprox(von_mises_stress(S), 100)

### Material parameters ###
# Isotropic elasticity: \dot{sigma} = \mathcal{C}:(\dot{\varepsilon}_{tot} - \dot{\varepsilon}_{pl})
E = 210000.0 # Young's modulus
nu = 0.3 # Poisson's ratio
K = E/(3*(1-2*nu)) # Bulk modulus
G = E/(2*(1+nu)) # Shear modulus
C = K*IxI + 2*G*P # Elasticity Tensor \mathcal{C}

@test isapprox(double_contraction(C, [0.001 0 0; 0 -nu*0.001 0; 0 0 -nu*0.001]), [E*0.001 0 0; 0 0 0; 0 0 0])

# Non-linear isotropic hardening: \dot{R} = b(Q-R)\dot{p}
# where \dot{p} = \sqrt{2/3 \dot{\varepsilon}_{pl}:\dot{\varepsilon}_{pl}}} - equivalent plastic strain rate
R0 = 100.0 # Initial proportionality limit
Q = 50.0 # Hardening magnitude
b = 0.1 # Hardening rate

# Non-linear kinematic hardening: \dot{X}_i = 2/3C_i\dot{p}(n - \frac{3D_i}{2C_i}X_i)
# where n = \frac{\partial f}{\partial \sigma} - plastic strain direction
# and X = \sum_{i=1}^N X_i
C_1 = 10000.0 # Slope parameter 1
D_1 = 100.0 # Rate parameter 1

C_2 = 50000.0 # Slope parameter 2
D_2 = 1000.0 # Rate parameter 2

# Viscoplasticity: Norton viscoplastic potential \phi = \frac{K_n}{n_n+1}\left( \frac{f}{K_n} \right)^{n_n+1}
# \dot{\varepsilon}_{pl} = \frac{\partial \phi}{\partial \sigma} = \frac{\partial \phi}{\partial f}\frac{\partial f}{\partial \sigma}
# => \dot{p} = \frac{\partial \phi}{\partial f} = \left( \frac{f}{K_n} \right)^n_n
# => n = \frac{\partial f}{\partial \sigma}
# => \dot{\varepsilon}_{pl} = \dot{p} n
K_n = 100.0 # Drag stress
n_n = 3.0 # Viscosity exponent

# Initialize variables
sigma = zeros(3,3)
R = R0
X_1 = zeros(3,3)
X_2 = zeros(3,3)
varepsilon_pl = zeros(3,3)
varepsilon_el = zeros(3,3)
t = 0.0

# Determine loading sequence
varepsilon_a = 0.01 # Strain amplitude
#varepsilon_tot(t) = sin(t)*[varepsilon_a 0 0; 0 -nu*varepsilon_a 0; 0 0 -nu*varepsilon_a]
dt = 0.01 # Time step
T0 = 1.0
T = 10.0 # Time span
# varepsilon_tot(t) = t/T*[varepsilon_a 0 0; 0 -nu*varepsilon_a 0; 0 0 -nu*varepsilon_a]
function varepsilon_tot(t)
    if t<T0
        return t/T0*[varepsilon_a 0 0; 0 -nu*varepsilon_a 0; 0 0 -nu*varepsilon_a]
    else
        return [varepsilon_a 0 0; 0 -nu*varepsilon_a 0; 0 0 -nu*varepsilon_a]
    end
end

# Initialize result storage
ts = [t]
sigmas = [sigma]
Rs = [R]
X_1s = [X_1]
X_2s = [X_2]
varepsilon_pls = [varepsilon_pl]
varepsilon_els = [varepsilon_el]

# Time integration
while t < T
    global t, sigma, R, X_1, X_2, varepsilon_pl, varepsilon_el, ts, sigmas, Rs, X_1s, X2_s, varepsilon_pls, varepsilon_els
    global C, K_n, n_n, C_1, D_1, C_2, D_2, Q, b
    # Store initial state
    sigma_n = sigma
    R_n = R
    X_1n = X_1
    X_2n = X_2
    varepsilon_pln = varepsilon_pl
    varepsilon_eln = varepsilon_el
    t_n = t

    # Increments
    t = t + dt
    dvarepsilon_tot = varepsilon_tot(t) - varepsilon_tot(t_n)
    # Elastic trial
    sigma_tr = sigma_n + double_contraction(C, dvarepsilon_tot)
    # Check for yield
    f_tr = von_mises_stress(sigma_tr - X_1 - X_2) - R
    println("***************************************")
    if f_tr <= 0 # Elastic step
        # Update variables
        println("Elastic step!")
        sigma = sigma_tr
        varepsilon_el += dvarepsilon_tot
    else # Viscoplastic step
        println("Viscoplastic step!")
        function g!(F, x) # System of non-linear equations
            sigma = reshape(x[1:9], 3,3)
            R = x[10]
            X_1 = reshape(x[11:19], 3,3)
            X_2 = reshape(x[20:28], 3,3)
            dotp = ((von_mises_stress(sigma - X_1 - X_2) - R)/K_n)^n_n
            dp = dotp*dt
            s = deviator(sigma - X_1 - X_2)
            n = 3/2*s/von_mises_stress(sigma - X_1 - X_2)
            dvarepsilon_pl = dp*n
            f1 = vec(sigma_n - sigma + double_contraction(C, dvarepsilon_tot - dvarepsilon_pl))
            f2 = R_n - R + b*(Q-R)*dp
            f3 = vec(X_1n - X_1 + 2/3*C_1*dp*(n - 3*D_1/(2*C_1)*X_1))
            f4 = vec(X_2n - X_2 + 2/3*C_2*dp*(n - 3*D_2/(2*C_2)*X_2))
            F[:] = vec([f1; f2; f3; f4])
        end
        x0 = vec([vec(sigma_tr); R; vec(X_1); vec(X_2)])
        F = similar(x0)
        res = nlsolve(g!, x0)
        x = res.zero
        sigma = reshape(x[1:9],3,3)
        R = x[10]
        X_1 = reshape(x[11:19], 3,3)
        X_2 = reshape(x[20:28], 3,3)
        dotp = ((von_mises_stress(sigma - X_1 - X_2) - R)/K_n)^n_n
        dp = dotp*dt
        s = deviator(sigma - X_1 - X_2)
        n = 3/2*s/von_mises_stress(sigma - X_1 - X_2)
        dvarepsilon_pl = dp*n
        varepsilon_pl += dvarepsilon_pl
        varepsilon_el += dvarepsilon_tot - dvarepsilon_pl
    end

    # Store variables
    push!(ts, t)
    push!(sigmas, sigma)
    push!(Rs, R)
    push!(X_1s, X_1)
    push!(X_2s, X_2)
    push!(varepsilon_pls, varepsilon_pl)
    push!(varepsilon_els, varepsilon_el)
end

qs = [von_mises_stress(sigma_i) for sigma_i in sigmas]
ps = [tr(sigma_i)/3 for sigma_i in sigmas]
xs = [von_mises_stress(X_1s[i] + X_2s[i]) for i in 1:length(ts)]
plot(ps, qs, label="Stress")
plot!(ps, xs+Rs, label="Static yield surface")
xlabel!("Hydrostatic stress")
ylabel!("Von Mises stress")


text:
"""
Functionality for handling the defautl conent of the various Genie files (migrations, models, controllers, etc).
"""
module FileTemplates

import Genie.Inflector

"""
    newtask(module_name::String) :: String

Default content for a new Genie Toolbox task.
"""
function newtask(module_name::String) :: String
  """
  module $module_name

  \"\"\"
  Description of the task here
  \"\"\"
  function runtask()
    # Build something great
  end

  end
  """
end

"""
    newcontroller(controller_name::String) :: String

Default content for a new Genie controller.
"""
function newcontroller(controller_name::String) :: String
  """
  module $(controller_name)Controller
    # Build something great
  end
  """
end

"""
    newtest(plural_name::String, singular_name::String) :: String

Default content for a new test file.
"""
function newtest(plural_name::String, singular_name::String) :: String
  """
  using Genie, App.$(plural_name)

  ### Your tests here
  @test 1 == 1
  """
end

"""
    appmodule(path::String)

Generates a custom app module when a new app is bootstrapped.
"""
function appmodule(path::String)
  path = replace(path, '-'=>'_') |> strip
  appname = split(path, '/', keepempty = false)[end] |> String |> Inflector.from_underscores

  content = """
  module $appname

  using Genie, Logging, LoggingExtras

  function main()
    Base.eval(Main, :(const UserApp = $appname))

    Genie.genie(; context = @__MODULE__)

    Base.eval(Main, :(const Genie = $appname.Genie))
    Base.eval(Main, :(using Genie))
  end

  end
  """

  (appname, content)
end

"""
    dockerfile(; user::String = "genie", supervisor::Bool = false, nginx::Bool = false, env::String = "dev",
                      filename::String = "Dockerfile", port::Int = 8000, dockerport::Int = 80, host::String = "0.0.0.0",
                      websockets_port::Int = port, websockets_dockerport::Int = dockerport)

Generates dockerfile for the Genie app.
"""
function dockerfile(; user::String = "genie", supervisor::Bool = false, nginx::Bool = false, env::String = "dev",
                      filename::String = "Dockerfile", port::Int = 8000, dockerport::Int = 80, host::String = "0.0.0.0",
                      websockets_port::Int = port, websockets_dockerport::Int = dockerport)
  appdir = "/home/$user/app"

  """
  FROM julia:latest

  # user
  RUN useradd --create-home --shell /bin/bash $user

  # app
  RUN mkdir $appdir
  COPY . $appdir
  WORKDIR $appdir

  RUN chown $user:$user -R *

  RUN chmod +x bin/repl
  RUN chmod +x bin/server
  RUN chmod +x bin/runtask

  USER $user

  RUN julia -e "using Pkg; Pkg.activate(\\".\\"); Pkg.instantiate(); Pkg.precompile(); "

  # ports
  EXPOSE $port
  EXPOSE $dockerport

  # websockets ports
  EXPOSE $websockets_port
  EXPOSE $websockets_dockerport

  ENV JULIA_DEPOT_PATH "/home/$user/.julia"
  ENV GENIE_ENV "$env"
  ENV HOST "$host"
  ENV PORT "$port"

  CMD ["bin/server"]
  """
end

end


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.Lorenz96
using StaticArrays

include("test.jl")

setRNGs(0)

model, theta, ys = Lorenz96.defaultLorenzModel(8, 100)

## just run the algorithm a few times

testSMC(model, 1024, 2, false)
testSMCParallel(model, 1024, 2, false)


text:
@testset "webserver" begin
    @async start_webserver()
    sleep(2)

    valid_params = JSON.json(TEST_PARAMETERS)

    @testset "valid plain input" begin
        response = HTTP.request("POST", SERVER_URI, ["Content-Type" => "application/json"],
                                valid_params)
        @test response.status == 200
    end

    @testset "invalid json plain input" begin
        try
            HTTP.request("POST", SERVER_URI, ["Content-Type" => "application/json"],
                                valid_params * "foo;baa")
        catch e
            @test e isa HTTP.ExceptionRequest.StatusError
            @test e.status == 400
        end
    end

    @testset "invalid content-type" begin
        try
            HTTP.request("POST", SERVER_URI, ["Content-Type" => "text/plain"],
                                valid_params)
        catch e
            @test e isa HTTP.ExceptionRequest.StatusError
            @test e.status == 400
        end
    end

    @testset "valid gzip input" begin
        response = HTTP.request("POST", SERVER_URI, ["Content-Type" => "application/json", "Content-Encoding" => "gzip"],
                                transcode(GzipCompressor, valid_params))
        @test response.status == 200
    end

    @testset "invalid content-encoding" begin
        try
            HTTP.request("POST", SERVER_URI, ["Content-Type" => "application/json", "Content-Encoding" => "compress"],
                                transcode(GzipCompressor, valid_params))
        catch e
            @test e isa HTTP.ExceptionRequest.StatusError
            @test e.status == 400
        end
    end
end


text:
"Return text with all redundant whitespace removed"
function remove_redundant_whitespace(text)
    cleaned_text = string(text[1])
    previous_c = text[1]
    for i in 2:length(text)
        c = text[i]
        if !(previous_c == ' ' && c == ' ')
            cleaned_text *= string(c)
        end
        previous_c = c
    end
    # Finally remove whitespaces from start and end
    return strip(cleaned_text)
end

"Utility function which turns vector of vectors into vector by concatenating them together"
function flatten(arr::Vector{Vector{T}}) where {T}
    flattened_arr = arr[1]
    for i = 2:length(arr)
        append!(flattened_arr, arr[i])
    end
    return flattened_arr
end

# If it is not Vector of Vectors, return it unchanged
flatten(arr::Vector) = arr

"Split text at stopwords to produce a vector of strings"
function split_at_stopwords(text::Union{Vector{SubString{String}},Vector{String}}, stopwords::Vector{String})::Vector{Vector{SubString{String}}}

    # A phrase is a vector of words
    phrases = Vector{Vector{SubString{String}}}(0)
    current_phrase = Vector{SubString{String}}(0)

    for sentence in text

        words_in_sentence = split(sentence, ' ')
        for word in words_in_sentence
            stopword_found = false
            for stopword in stopwords
                if lowercase(word) == stopword
                    stopword_found = true
                    break
                end
            end

            if stopword_found
                if length(current_phrase) > 0
                    push!(phrases, current_phrase)
                    current_phrase = Vector{SubString{String}}(0)
                end
            else
                push!(current_phrase, word)
            end
        end

        if length(current_phrase) > 0
            push!(phrases, current_phrase)
        end

    end

    return phrases
end

"Take vector of strings and split each string into all possible keyphrases"
function find_all_possible_keyphrases(phrases, max_phrase_length::Int)
    all_possible_keyphrases = Array{String}(0)
    for phrase in phrases
        for j = 0: max_phrase_length-1
            for i = 1:(length(phrase)-j)
                push!(all_possible_keyphrases, join(phrase[i:i+j], " "))
            end
        end
    end
    return all_possible_keyphrases
end

"Score all keywords using a variety of scoring techniques"
function calculate_keyphrase_scores(all_possible_keyphrases)
    keyphrases = unique(all_possible_keyphrases)
    d = Dict([(i, count(x->x==i,all_possible_keyphrases)) for i in keyphrases])
    return d
end

"Utility function to test if string is a float"
function isa_numberstring(s)
    try
        parse(Float64, s)
        return True
    catch
        return False
    end
end

max_stars_repo_name: EtymoIO/RAKE

max_issues_repo_name: EtymoIO/RAKE


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

#Generic benchmark driver
for (testfunc, testname, longtestname, problem_sizes) in testdata
    for (n, t, size) in problem_sizes
        @timeit testfunc(n, t) string(testname,"_",size) string(uppercase(size[1]),size[2:end]," ",longtestname," test")
    end
end

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
function SymmetricGaussIntegrator(arg0::Pair)
    return SymmetricGaussIntegrator((Pair,), arg0)
end

function SymmetricGaussIntegrator(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return SymmetricGaussIntegrator((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function integrate(obj::SymmetricGaussIntegrator, arg0::UnivariateFunction)
    return jcall(obj, "integrate", jdouble, (UnivariateFunction,), arg0)
end


text:
#=
Proof of concept of 1D localization against a scalar field.
In this example we use a sequence of 1D measurements and a known map with a
correlator-like approach to localize.

2021-03-10_16
=#

using Images, FileIO           # load elevation from PNG
using ImageFiltering
using DSP
using DataInterpolations       # handles the terrain
using Distributions
using Statistics
using Caesar, RoME
using Gadfly, Colors
using Cairo
using RoMEPlotting
using Plots

## Preamble: creating the measurement sequence factor

# 1. This is the (constant) scalar field the elevation measurements refer to
#
# struct ScalarField1D
#     # scalar field over one variable f:R->R (i.e. h= f(x))
#     # [a special case of ScalarField, h = f(·)]
#     # ...
# end

# 2. This is the actual scalar measurement sequence factor
# It takes two arguments:
#  - the elevation measurement sequence, and
#  - the odometry estimate at each elevation measurement
#
# addFactor!(fg, [:x1, :topography], ScalarFieldSequenceFactor1D(x,y))
struct ScalarFieldSequenceFactor1D <: IIF.AbstractPrior
  # FIXME: objective is to <: IIF.AbstractRelativeConstant (final name TBD)
  gridSequence::Vector{Float64} # sample location (odometry)
  scalarSequence::Vector{Float64} #sample value (elevation meas)
end

# 3. Objective/fitness function evaluating match between a map and a template for different displacement values
# This is used within the factor to generate the likelihood over the terrain
#
# NOTE: this function assumes map and template are on an equivalent grid
# this function does not care for where/what the grid is
# NOTE: result will assume base pose is first in sequence
# for each i
#   energy[i] = sum( (map[i:(i+length(template))] - template).^2 )
# density = exp.(-energy)     # Woodward
function _mySSDCorr(map, template)

    maplen = length(map)
    len = length(template)

    mnm = Statistics.mean(map)
    map_ = [map; mnm*ones(len-1)]
    density = zeros(maplen)
    for i in 1:maplen
        density[i] = -sum( (map_[i:(i+len-1)] .- template).^2 )
    end
    adaptive = abs(maximum(density))

    density ./= adaptive
    density .= exp.(density)
    density ./= sum(density)
    return density
end

function IIF.getSample(s::CalcFactor{<:ScalarFieldSequenceFactor1D}, N::Int=1)
    # locality if desired from current variable estimate,
    # X0 = getBelief(s.metadata.fullvariables[1]) # to consider only local map info

    # assuming user addConstant!(fg, :terrain, ScalarField1D(dem))

    global terrain # should come from s.constants[:terrain] ?????

    # ensure equal spacing as terrain
    # assumptions:
    # - terrain.t is equally spaced
    # - s.factor.gridSequence increases monotonically
    # s.factor.gridSequence is the sequence of odometry estimates at which
    #  the elevation measurements (s.factor.scalarSequence) were taken
    measurement = LinearInterpolation(s.factor.scalarSequence,s.factor.gridSequence)

    upsample = 1 # avoid quantization effects from gridded sequence/terrain data
    dx = diff(terrain.t)[1]/upsample
    template = measurement.(s.factor.gridSequence[1]:dx:s.factor.gridSequence[end])
    map_ = [terrain.u;]
    map2_grid = range(terrain.t[1],terrain.t[end],length=upsample*length(map_))
    map2 = terrain.(map2_grid)

    intensity_ = _mySSDCorr(map2, template)

        # intensity_ = xcorr([terrain.u;], template, padmode=:none)
        # chop_ = round(Int, length(template)/2)
        # intensity = intensity_[chop_:(end-chop_)]
        # # force equal length
        # while length(intensity) < length(terrain.t)
        #     push!(intensity, 0)
        # end
        # intensity = intensity[1:length(terrain.t)]

    # AliasingSS only works for 1D at this time
    # NOTE, if you stray out of region of support, things blow up!
    bss = AliasingScalarSampler([map2_grid;], intensity_) # [terrain.t;]

    # buffer with kde for continuous sampling (not grid spaced)
    stagedsmpls = rand(bss,N)
    smpls = rand(kde!(stagedsmpls, [3*dx;]),N)

    return (reshape(smpls,1,N),)
end

# # keep track of IIF #1051
# # parameters are both in map/trajectory space X
# function (s::CalcFactor{<:ScalarFieldSequenceFactor1D})(likelihoodSample, pose)
#     return likelihoodSample - pose
# end

## Terrain

# load dem (18x18km span)
img = load(joinpath(dirname(dirname(pathof(Caesar))), "examples/dev/scalar/dem.png")) .|> Gray
h = 1e3*Float64.( @view img[512,:])

## Size of the map
# x = range(-9000,9000,length = length(h)) # actual span is 18km on each axis
x = range(-100,100,length = length(h))
global terrain
terrain = LinearInterpolation(h, x)
Plots.plot(terrain)

## Simulation parameters
seq_length = 10 # measurement sequence length (will get upsampled to the grid size!)
xt = -70        # true pose (initial value)
delta = 1       # pose increment
sigma_x = 0.1   # odometry measurement uncertainty
sigma_y = 1.0   # elevation measurement uncertainty

fg = initfg() # create empty factor graph
# getSolverParams(fg).inflation = 3.0

addVariable!(fg, :x0, ContinuousEuclid{1})  # add first pose

# add terrain to factor graph (missing for now)
# addConstant!(fg, :em, ScalarField1D(dem))

# simulate for a few poses
xt_seq = zeros(0)
yt_seq = zeros(0)
x_seq = zeros(0)
y_seq = zeros(0)
for i=1:10
    @show i
    xt_seq = zeros(0)
    yt_seq = zeros(0)
    x_seq = zeros(0)
    y_seq = zeros(0)

    # build up sequence
    for j=1:seq_length
        xt = xt + delta         # true position
        yt = terrain(xt)  # true elevation

        append!(xt_seq,xt)
        append!(yt_seq,yt)

        # append to sequence
        append!(x_seq, xt .+ sigma_x*randn(1))
        append!(y_seq, yt .+ sigma_y*randn(1))
    end

    addVariable!(fg, Symbol("x$(i)"), ContinuousEuclid{1}) # BoundedEuclid{1}

    # add propagation
    labels = [Symbol("x$j") for j = (i-1):i]
    addFactor!(fg, labels, LinearRelative(Normal(delta*seq_length, 5)),graphinit=false)

    # grid, scalar
    em = ScalarFieldSequenceFactor1D( x_seq , y_seq )  # .- x_seq[1]

    # addFactor!(fg, [:terrain, Symbol("y$(i-1)")], em)        # The Right Way (TM)
    addFactor!(fg, [Symbol("x$(i)");], em, graphinit=false) # the hacky version
    @info "i=$i, ends"
end

## eval the sampler at first prior
PL = []
for lb in (ls(fg) |> sortDFG)[2:end]
    pts = approxConv(fg, Symbol(lb, "f1"), lb)
    prop = kde!(pts)
    pl = plotKDE(prop)

    pl2 = Gadfly.plot(
        Gadfly.layer(x=terrain.t, y=terrain.u, Geom.line, Theme(default_color=colorant"red")),
        Gadfly.layer(x=xt_seq,y=yt_seq, Geom.line, Theme(default_color=colorant"green")),
        Gadfly.layer(x=x_seq, y=y_seq, Geom.line, Theme(default_color=colorant"blue"))
    )

    push!(PL, vstack(pl, pl2))
end

PL[1]

## Solve
getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true   # still experimental, treeinit

tree, _, = solveTree!(fg);

## Show results

smplOnMap = approxConv(fg, :x10f1, :x10)
histogram(smplOnMap[:], bins=100)

plot(terrain)
plot!([-100; 100],[436; 436])

# all poses
Gadfly.set_default_plot_size(25cm,20cm)
pl1 = plotKDE(fg, sortDFG(ls(fg)))             # plot posteriors for all variables
pl1.coord = Coord.cartesian(xmin=-100,xmax=85)
pl2 = Gadfly.plot(x=terrain.t,y=terrain.u, Geom.line)
vstack(pl1, pl2)

##


text:
#unitary test GenericStatemod
x0 = ones(6)
state0 = GenericState(x0)

@test isnan(state0.start_time) #Default value of start_time is 0
x1 = [1.0]
update!(state0, x = x1) #Check the update of state0
@test state0.x == x1
@test isnan(state0.start_time) #start_time must be unchanged

update!(state0, tmps = 1.0)
@test state0.x == x1 #must be unchanged
@test state0.start_time == 1.0


text:
#   Purpose:   Demonstrate how to solve the following
#              optimization problem
#
#              minimize    f(x) + c^T x
#              subject to  l^c <= a x <= u^c
#                          l^x <=  x  <= u^x
#              where f is a nonlinear (convex) function.
#
#              The program solves the example:
#
#              minimize    - x_0 - ln(x_1 + x_2)
#              subject to  x_0 + x_1 + x_2 = 1
#                          x_0, x_1, x_2 >= 0.
#

using Mosek

maketask() do t
    putstreamfunc(t,MSK_STREAM_LOG,msg -> print(msg))

    appendvars(t,3)
    appendcons(t,1)

    putvarbound(t,1,MSK_BK_LO,0.0, +Inf)  # x0
    putvarbound(t,2,MSK_BK_LO,0.0, +Inf)  # x1
    putvarbound(t,3,MSK_BK_LO,0.0, +Inf)  # x2

    putarow(t,1,[1,2,3],[1.0,1.0,1.0])   # x0 + x1 + x2
    putconbound(t,1,MSK_BK_FX,1.0, 1.0)  # = 1.0

    putclist(t,[1], [-1.0])

    #define the evaluation functions
    function indexof(v,a)
        for i in 1:length(a)
            if a[i] == v return i
            end
        end
        return 0
    end

    evalobj(x :: Array{Float64,1}) = log(x[2]+x[3])

    evalconi(x:: Array{Float64,1},i:: Int32) = 0.0

    function grdobj(x :: Array{Float64,1},sub:: Array{Int32,1}, val:: Array{Float64,1})
        for j in 1:length(sub)
            if sub[j] == 2 || sub[j] == 3
                val[j] = -1.0 / (x[2]+x[3])
            end
        end
    end

    function grdconi(x  :: Array{Float64,1},
                     i  :: Int32,
                     sub:: Array{Int32,1},
                     val:: Array{Float64,1})
        none
    end

    function grdlag(x ::   Array{Float64,1},
                    yo::   Float64,
                    yc::   Array{Float64,1},
                    subi:: Array{Int32,1},
                    val::  Array{Float64,1})
        val[2] = - yo * 1.0 / (x[2]+x[3])
        val[3] = - yo * 1.0 / (x[2]+x[3])
    end

    # NOTE: in this case we have no off-diagonal nonzeros in the hessian.
    # If we do, then _either_ element (i,j) _or_ element (j,i) should be
    # specified.
    function heslag(x ::      Array{Float64,1},
                    yo::      Float64,
                    yc::      Array{Float64,1},
                    subi::    Array{Int32,1},
                    hessubi:: Array{Int32,1},
                    hessubj:: Array{Int32,1},
                    hesval::  Array{Float64,1})

        hessubi[1] = 2; hessubj[1] = 2; hesval[1] = (x[2]+x[3])^(-2)
        hessubi[2] = 3; hessubj[2] = 2; hesval[2] = (x[2]+x[3])^(-2)
        hessubi[3] = 3; hessubj[3] = 3; hesval[3] = (x[2]+x[3])^(-2)
    end

    putnlcallbacks(t,
                   Int[2,3], # subscripts of non-zeros in the gradient of the objective
                   Int[], # subscripts of non-zeros in the gradient of the constraints
                   Int[1,1], # rowptr for subscripts of non-zeros in the gradient of the constraints
                   Int[2,3,3], # hessubi
                   Int[2,2,3], # hessubj
                   evalobj,
                   evalconi,
                   grdlag,
                   grdobj,
                   grdconi,
                   heslag)

    optimize(t)

    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(t,MSK_STREAM_MSG)

    # Get status information about the solution
    solsta = getsolsta(t,MSK_SOL_ITR)

    if solsta in     [ MSK_SOL_STA_OPTIMAL,
                      MSK_SOL_STA_NEAR_OPTIMAL ]
        xx = getxx(t,MSK_SOL_ITR)
        println("Optimal solution:\n")
        show(xx)
        println()
    elseif solsta in [ MSK_SOL_STA_DUAL_INFEAS_CER,
                      MSK_SOL_STA_PRIM_INFEAS_CER,
                      MSK_SOL_STA_NEAR_DUAL_INFEAS_CER,
                      MSK_SOL_STA_NEAR_PRIM_INFEAS_CER ]
        println("Primal or dual infeasibility certificate found.\n")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status\n")
    else
        @printf("Other solution status (%d)\n",solsta)
    end

end


text:
import LanguageServer: LanguageServerInstance, Document
server = LanguageServerInstance(IOBuffer(), IOBuffer(), false)
LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse(init_request)), server)

function getresult(server)
    str = String(take!(server.pipe_out))
    JSON.parse(str[search(str, '{'):end])["result"]["contents"]
end

testtext = """
module testmodule
struct testtype
    a
    b::Float64
    c::Vector{Float64}
end

function testfunction(a, b::Float64, c::testtype)
    return c
end
end
testmodule
"""

server.documents[LanguageServer.URI2("testdoc")] = Document("testdoc", testtext, true)
doc = server.documents[LanguageServer.URI2("testdoc")]
LanguageServer.parse_all(doc, server)

# clear init output
take!(server.pipe_out)

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse("""{"jsonrpc":"2.0","id":1,"method":"textDocument/hover","params":{"textDocument":{"uri":"testdoc"},"position":{"line":11,"character":5}}}""")), server)

res = getresult(server)

@test res[1]["value"] == "module"

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse("""{"jsonrpc":"2.0","id":1,"method":"textDocument/hover","params":{"textDocument":{"uri":"testdoc"},"position":{"line":1,"character":9}}}""")), server)

res = getresult(server)

@test res[1]["value"] == "struct"

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse("""{"jsonrpc":"2.0","id":1,"method":"textDocument/hover","params":{"textDocument":{"uri":"testdoc"},"position":{"line":7,"character":19}}}""")), server)

res = getresult(server)

@test res[1]["value"] == "function"


text:
"""
Format the shapes from the opf and return a [`Dict`] with the id as key and the shape data as values.

# Exemples
> opf= read_opf("simple_OPF_shapes.opf")
> extract_opf_shapes(opf)
"""
function extract_opf_shapes(opf)
    shapes_Dict= Dict{Int32,Any}()
    shapes= opf["shapeBDD"]["shape"]
    for i in shapes
        push!(shapes_Dict, i[:Id] => i)
    end
    return shapes_Dict
end


text:
module NeuroAnalysis

using FileIO, JLD2
using Plots
gr(fmt="png", size=(600, 400))

include("IO/io.jl")
include("Spike/Spike.jl")

end  # module NeuroAnalysis


text:
module MeasureTheory

using Random

using MeasureBase
using ConcreteStructs
using MLStyle
using NestedTuples
using TransformVariables
const TV = TransformVariables

import Base
import Distributions
const Dists = Distributions

# export TV
export ≪
export sampletype
export For

export AbstractMeasure
export Dirac
export Lebesgue
export ℝ, ℝ₊, 𝕀
export ⊙
export SpikeMixture
export CountingMeasure
export TrivialMeasure
export Likelihood
export testvalue

using InfiniteArrays
using ConcreteStructs
using DynamicIterators
using KeywordCalls
using ConstructionBase
using Accessors
using StatsFuns
using SpecialFunctions
using LogExpFunctions

import MeasureBase: testvalue, logdensity, density, basemeasure, kernel, params, ∫

using Reexport
@reexport using MeasureBase

using Tricks: static_hasmethod
const ∞ = InfiniteArrays.∞

export ∞

export as
export Affine
export AffineTransform

if VERSION < v"1.7.0-beta1.0"
    @eval begin
        struct Returns{T}
            value::T
        end

        (f::Returns)(x) = f.value
    end
end

sampletype(μ::AbstractMeasure) = typeof(testvalue(μ))

# sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))

import Distributions: pdf, logpdf

export pdf, logpdf

function logpdf(d::AbstractMeasure, x)
    _logpdf(d, x, logdensity(d,x))
end

function _logpdf(d::AbstractMeasure, x, acc::Real)
    β = basemeasure(d)
    d === β && return acc

    _logpdf(β, x, acc + logdensity(β, x))
end

pdf(d::AbstractMeasure, x) = exp(logpdf(d, x))

"""
    logdensity(μ::AbstractMeasure [, ν::AbstractMeasure], x::X)

Compute the logdensity of the measure μ at the point x. This is the standard way
to define `logdensity` for a new measure. the base measure is implicit here, and
is understood to be `basemeasure(μ)`.
"""
function logdensity end

include("const.jl")
# include("traits.jl")
include("parameterized.jl")
# include("resettablerng.jl")

include("combinators/weighted.jl")
include("combinators/product.jl")
include("combinators/transforms.jl")
include("combinators/chain.jl")
# include("combinators/basemeasure.jl")

include("distributions.jl")

include("parameterized/normal.jl")
include("parameterized/studentt.jl")
include("parameterized/cauchy.jl")
include("parameterized/categorical.jl")
include("parameterized/laplace.jl")
include("parameterized/uniform.jl")
include("parameterized/beta.jl")
include("parameterized/dirichlet.jl")
include("parameterized/gumbel.jl")
include("parameterized/exponential.jl")
include("parameterized/mvnormal.jl")
# include("parameterized/inverse-gamma.jl")
include("parameterized/bernoulli.jl")
include("parameterized/poisson.jl")
include("parameterized/binomial.jl")
include("parameterized/multinomial.jl")
include("parameterized/lkj-cholesky.jl")
include("parameterized/negativebinomial.jl")

include("transforms/corrcholesky.jl")
include("transforms/ordered.jl")

include("distproxy.jl")
end # module


text:
# Generate some test data
t = collect(0:0.1:5.0)
f_(t) = [sin(t); cos(t)]
df_(t) = [cos(t); -sin(t)]
X = hcat(map(i -> f_(i), t)...)
Y = randn(2,2)*X
DX = hcat(map(i -> df_(i), t)...)
p = [2.7]
u_(u, p, t) = [p[1] * u[1] - 3.0 * t]
U = hcat(map(i -> u_(X[:, i], p, t[i]), 1:length(t))...)

@testset "DiscreteProblem" begin
    p1 = DiscreteDataDrivenProblem(X)
    p2 = DiscreteDataDrivenProblem(X, t)
    p3 = DiscreteDataDrivenProblem(X, t, U)
    p4 = DiscreteDataDrivenProblem(X, t, u_, p = p)

    for p_ in [p1;p2;p3;p4]
        @test is_valid(p_)
        @test is_discrete(p_)
        @test_nowarn @is_applicable p_
    end
    @test isequal(p3.U, p4.U)
    @test isequal(p4.p, p)
end

@testset "Collocation" begin
    t_ = 0.0:0.001:10.0
    x_ = permutedims(sin.(t_))
    y_ = permutedims(cos.(t_))
    for m in [EpanechnikovKernel, UniformKernel,TriangularKernel,QuarticKernel,
        TriweightKernel, TricubeKernel, GaussianKernel, CosineKernel,
        LogisticKernel, SigmoidKernel, SilvermanKernel]
        p_ = ContinuousDataDrivenProblem(x_, t_)
        p_2 = ContinuousDataDrivenProblem(x_, t_, crop = true)
        @test norm(y_ .- p_.DX) < 1e-1
        @test norm(y_[:, 2:end-1] - p_2.DX) < 1e-1
    end
end

@testset "ContinuousProblem" begin
    p2 = ContinuousDataDrivenProblem(X, t, GaussianKernel())
    p3 = ContinuousDataDrivenProblem(X, DX)
    p4 = ContinuousDataDrivenProblem(X, t, U)
    p5 = ContinuousDataDrivenProblem(X, t, DX, U)
    p6 = ContinuousDataDrivenProblem(X, t, DX, u_, p = p)

    for p_ in [p2;p3;p4;p5;p6]
        @test is_valid(p_)
        @test is_continuous(p_)
    end

    @test isapprox(p3.X, p2.X, atol = 0.3)
    @test isapprox(p3.DX, p2.DX, atol = 0.3)
    @test isequal(p5.U, p6.U)
    @test is_parametrized(p6)
    @test isequal(p6.p, p)
    @test has_timepoints(p4)
end

@testset "DirectProblem" begin
    p1 = DirectDataDrivenProblem(X, Y)
    p2 = DirectDataDrivenProblem(X, t, Y)
    p3 = DirectDataDrivenProblem(X, t, Y, U)
    p4 = DirectDataDrivenProblem(X, t, Y, u_, p = p)

    for p_ in [p1;p2;p3;p4]
        @test is_valid(p_)
        @test is_direct(p_)
    end

    @test isequal(p3.U, p4.U)
    @test isequal(p4.p, p)
    @test !is_autonomous(p3)
end

@testset "Problem Basis Interaction" begin
    @variables x y z t α β u
    b1 = Basis([α*x; β*y; z*t^2+u], [x; y; z], iv = t, parameters = [α; β], controls = [u])
    b2 = Basis([α*x; β*y; z*t; α], [x; y; z], iv = t, parameters = [α; β])
    sample_size = 100
    X1 = randn(3, sample_size)
    DX1 = randn(3, sample_size)
    DX2 = randn(3, sample_size-1)
    X2 = randn(4, sample_size)
    Y1 = randn(3, sample_size)
    Y2 = randn(4, sample_size)
    t = randn(sample_size)
    U = randn(1, sample_size)
    ps = randn(3)

    p1 = DirectDataDrivenProblem(X1, Y1, t = t, p = ps, U = U)
    p2 = DirectDataDrivenProblem(X1, Y2, t = t, p = ps, U = U)
    p3 = ContinuousDataDrivenProblem(X1, t, DX1,  Y = Y1, p = ps, U = U)
    p4 = DiscreteDataDrivenProblem(X1, Y = Y2, t = t, p = ps, U = U)

    @testset "Check validity" begin
        @test_throws AssertionError @is_applicable p2 b2
        @test_throws AssertionError @is_applicable p1 b2
        @test_throws AssertionError @is_applicable p3 b2
        @test_throws AssertionError @is_applicable p4 b2
        @test_throws AssertionError @is_applicable p4 b1 DX1

        @test_nowarn @is_applicable p1 b1
        @test_nowarn @is_applicable p3 b1
        @test_nowarn @is_applicable p4 b1
        @test_nowarn @is_applicable p1 b1 DX1
        @test_nowarn @is_applicable p3 b1 DX1
    end

    @testset "Apply problem" begin
        DX1 .= 0.0
        @test_nowarn b1(DX1, p1)
        @test iszero(norm(DX1 - b1(p1)))

        DX1 .= 0.0
        @test_nowarn b1(DX1, p3)
        @test iszero(norm(DX1 - b1(p3)))

        DX2 .= 0.0
        @test_nowarn b1(DX2, p4)
        @test iszero(norm(DX2 - b1(p4)))

        DX1 .= 0.0
        @test_nowarn @views b1(DX1[:,1:end-1], p4)
        @test iszero(norm(DX1[:, 1:end-1] - b1(p4)))
    end
end

@testset "DataDrivenDataset" begin
    p1 = ContinuousDataDrivenProblem(X, t)
    p2 = ContinuousDataDrivenProblem(X, t, DX = DX)

    data = (
        prob1 = (X = X, t = t, Y = Y),
        prob2 = (X = X, t = t, Y = Y, DX = DX)
    )

    s1 = DataDrivenDataset(p1, p2)
    s2 = ContinuousDataset(data)
    s3 = DirectDataset(data)
    s4 = DiscreteDataset(data)

    sets = [s1, s2, s3, s4]

    # Information
    @test is_discrete(s4)
    @test is_continuous(s2)
    @test is_direct(s3)

    # Sizes
    for s in sets
        @test size(s) == (first(size(p1)), is_discrete(s) ? 2*size(X, 2)-2 : 2*size(X, 2))
        @test DataDrivenDiffEq.is_valid(s)
    end

    # Targets
    @test DataDrivenDiffEq.get_target(s1) == DataDrivenDiffEq.get_target(s2)
    @test DataDrivenDiffEq.get_target(s1) == hcat(DataDrivenDiffEq.get_target(p1), DataDrivenDiffEq.get_target(p2))
    @test DataDrivenDiffEq.get_target(s3) == hcat(Y, Y)
    @test DataDrivenDiffEq.get_target(s4) == hcat(X[:, 2:end], X[:, 2:end])

    # Args
    @test DataDrivenDiffEq.get_oop_args(s1) == DataDrivenDiffEq.get_oop_args(s2)
    @test DataDrivenDiffEq.get_implicit_oop_args(s1) == DataDrivenDiffEq.get_implicit_oop_args(s2)

    # Basis handling
    @variables x[1:size(X, 1)]
    b = Basis(x, x)
    @test b(s1) == hcat(b(p1), b(p2))
    @test b(s2) == b(s1)
    @test hcat(X, X) == b(s3)
    @test hcat(X[:, 1:end-1], X[:, 1:end-1]) == b(s4)

    # Check if misspecified data is detected
    wrong_data = (
        prob1 = (X = X, Y = Y),
        prob2 = (X = X, t = t, Y = Y)
        )
    @test_throws ArgumentError ContinuousDataset(wrong_data)
end