text:
# See **[`ImagePreprocessing`](#).**

"""
    ImagePreprocessing(means, stds[; C = RGB{N0f8}, T = Float32])

Converts an image to a color `C`, then to a 3D-array of type `T` and
finally normalizes the values using `means` and `stds`.

If no `means` or `stds` are given, uses ImageNet statistics.
"""
struct ImagePreprocessing
    # hold one copy of the transform for every context
    # in case the input sizes differ based on `context`
    # since that would create problems with the buffers
    traintfm
    validtfm
    inferencetfm
end

function ImagePreprocessing(
        means::SVector{N} = IMAGENET_MEANS,
        stds::SVector{N} = IMAGENET_STDS;
        C = RGB{N0f8},
        T = Float32,
        buffered = true) where N
    # TODO: tensor of type T
    tfms = (
        ToEltype(C) |> ImageToTensor() |> Normalize(means, stds),
        ToEltype(C) |> ImageToTensor() |> Normalize(means, stds),
        ToEltype(C) |> ImageToTensor() |> Normalize(means, stds),
    )
    if buffered
        tfms = (
            BufferedThreadsafe(tfms[1]),
            BufferedThreadsafe(tfms[2]),
            # Inference transform is not buffered since it can have
            # varying sizes
            tfms[3]
        )
    end

    return ImagePreprocessing(tfms...)
end

function ImagePreprocessing(means::NTuple{N}, stds::NTuple{N}; kwargs...) where N
    return ImagePreprocessing(SVector{N}(means), SVector{N}(stds); kwargs...)
end

Base.show(io::IO, ::ImagePreprocessing) = print(io, "ImagePreprocessing()")

function run(ip::ImagePreprocessing, context::Context, image)
    tfm = _gettfm(ip, context)
    return DataAugmentation.apply(tfm, DataAugmentation.Image(image)) |> itemdata
end

function run!(x, ip::ImagePreprocessing, context::Context, image)
    tfm = _gettfm(ip, context)
    DataAugmentation.apply!(ArrayItem(x), tfm, DataAugmentation.Image(image)) |> itemdata
    return x
end

_gettfm(ip::ImagePreprocessing, context::Training) = ip.traintfm
_gettfm(ip::ImagePreprocessing, context::Validation) = ip.validtfm
_gettfm(ip::ImagePreprocessing, context::Inference) = ip.inferencetfm

function invert(ip::ImagePreprocessing, x::AbstractArray{T, N}) where {N, T}
    tfm = _parenttfm(ip.traintfm).transforms[end]
    return DataAugmentation.tensortoimage(DataAugmentation.denormalize(x, tfm.means, tfm.stds))
end

const IMAGENET_MEANS = SVector{3, Float32}(.485, 0.456, 0.406)
const IMAGENET_STDS = SVector{3, Float32}(0.229, 0.224, 0.225)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: braket
using AWS.Compat
using AWS.UUIDs

"""
    cancel_quantum_task(client_token, quantum_task_arn)
    cancel_quantum_task(client_token, quantum_task_arn, params::Dict{String,<:Any})

Cancels the specified task.

# Arguments
- `client_token`: The client token associated with the request.
- `quantum_task_arn`: The ARN of the task to cancel.

"""
cancel_quantum_task(clientToken, quantumTaskArn; aws_config::AbstractAWSConfig=global_aws_config()) = braket("PUT", "/quantum-task/$(quantumTaskArn)/cancel", Dict{String, Any}("clientToken"=>clientToken); aws_config=aws_config)
cancel_quantum_task(clientToken, quantumTaskArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("PUT", "/quantum-task/$(quantumTaskArn)/cancel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>clientToken), params)); aws_config=aws_config)

"""
    create_quantum_task(action, client_token, device_arn, output_s3_bucket, output_s3_key_prefix, shots)
    create_quantum_task(action, client_token, device_arn, output_s3_bucket, output_s3_key_prefix, shots, params::Dict{String,<:Any})

Creates a quantum task.

# Arguments
- `action`: The action associated with the task.
- `client_token`: The client token associated with the request.
- `device_arn`: The ARN of the device to run the task on.
- `output_s3_bucket`: The S3 bucket to store task result files in.
- `output_s3_key_prefix`: The key prefix for the location in the S3 bucket to store task
  results in.
- `shots`: The number of shots to use for the task.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"deviceParameters"`: The parameters for the device to run the task on.
- `"tags"`: Tags to be added to the quantum task you're creating.
"""
create_quantum_task(action, clientToken, deviceArn, outputS3Bucket, outputS3KeyPrefix, shots; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/quantum-task", Dict{String, Any}("action"=>action, "clientToken"=>clientToken, "deviceArn"=>deviceArn, "outputS3Bucket"=>outputS3Bucket, "outputS3KeyPrefix"=>outputS3KeyPrefix, "shots"=>shots); aws_config=aws_config)
create_quantum_task(action, clientToken, deviceArn, outputS3Bucket, outputS3KeyPrefix, shots, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/quantum-task", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("action"=>action, "clientToken"=>clientToken, "deviceArn"=>deviceArn, "outputS3Bucket"=>outputS3Bucket, "outputS3KeyPrefix"=>outputS3KeyPrefix, "shots"=>shots), params)); aws_config=aws_config)

"""
    get_device(device_arn)
    get_device(device_arn, params::Dict{String,<:Any})

Retrieves the devices available in Amazon Braket.

# Arguments
- `device_arn`: The ARN of the device to retrieve.

"""
get_device(deviceArn; aws_config::AbstractAWSConfig=global_aws_config()) = braket("GET", "/device/$(deviceArn)"; aws_config=aws_config)
get_device(deviceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("GET", "/device/$(deviceArn)", params; aws_config=aws_config)

"""
    get_quantum_task(quantum_task_arn)
    get_quantum_task(quantum_task_arn, params::Dict{String,<:Any})

Retrieves the specified quantum task.

# Arguments
- `quantum_task_arn`: the ARN of the task to retrieve.

"""
get_quantum_task(quantumTaskArn; aws_config::AbstractAWSConfig=global_aws_config()) = braket("GET", "/quantum-task/$(quantumTaskArn)"; aws_config=aws_config)
get_quantum_task(quantumTaskArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("GET", "/quantum-task/$(quantumTaskArn)", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Shows the tags associated with this resource.

# Arguments
- `resource_arn`: Specify the resourceArn for the resource whose tags to display.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = braket("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    search_devices(filters)
    search_devices(filters, params::Dict{String,<:Any})

Searches for devices using the specified filters.

# Arguments
- `filters`: The filter values to use to search for a device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned in the response. Use the
  token returned from the previous request continue results where the previous request ended.
"""
search_devices(filters; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/devices", Dict{String, Any}("filters"=>filters); aws_config=aws_config)
search_devices(filters, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/devices", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filters"=>filters), params)); aws_config=aws_config)

"""
    search_quantum_tasks(filters)
    search_quantum_tasks(filters, params::Dict{String,<:Any})

Searches for tasks that match the specified filter values.

# Arguments
- `filters`: Array of SearchQuantumTasksFilter objects.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned in the response. Use the
  token returned from the previous request continue results where the previous request ended.
"""
search_quantum_tasks(filters; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/quantum-tasks", Dict{String, Any}("filters"=>filters); aws_config=aws_config)
search_quantum_tasks(filters, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/quantum-tasks", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filters"=>filters), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Add a tag to the specified resource.

# Arguments
- `resource_arn`: Specify the resourceArn of the resource to which a tag will be added.
- `tags`: Specify the tags to add to the resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Remove tags from a resource.

# Arguments
- `resource_arn`: Specify the resourceArn for the resource from which to remove the tags.
- `tag_keys`: pecify the keys for the tags to remove from the resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = braket("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = braket("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

ext: jl


text:
export preprocess_one_hot_encoding

using mlpack._Internal.io

import mlpack_jll
const preprocess_one_hot_encodingLibrary = mlpack_jll.libmlpack_julia_preprocess_one_hot_encoding

# Call the C binding of the mlpack preprocess_one_hot_encoding binding.
function preprocess_one_hot_encoding_mlpackMain()
  success = ccall((:preprocess_one_hot_encoding, preprocess_one_hot_encodingLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module preprocess_one_hot_encoding_internal
  import ..preprocess_one_hot_encodingLibrary

end # module

"""
    preprocess_one_hot_encoding(dimensions, input; [verbose])

This utility takes a dataset and a vector of indices and does one-hot encoding
of the respective features at those indices. Indices represent the IDs of the
dimensions to be one-hot encoded.

The output matrix with encoded features may be saved with the `output`
parameters.

So, a simple example where we want to encode 1st and 3rd feature from dataset
`X` into `X_output` would be

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> X_ouput = preprocess_one_hot_encoding(1, X)
```

# Arguments

 - `dimensions::Vector{Int}`: Index of dimensions thatneed to be one-hot
      encoded.
 - `input::Array{Float64, 2}`: Matrix containing data.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Float64, 2}`: Matrix to save one-hot encoded features
      data to.

"""
function preprocess_one_hot_encoding(dimensions::Vector{Int},
                                     input;
                                     verbose::Union{Bool, Missing} = missing,
                                     points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, preprocess_one_hot_encodingLibrary), Nothing, ());

  IORestoreSettings("One Hot Encoding")

  # Process each input argument before calling mlpackMain().
  IOSetParam("dimensions", dimensions)
  IOSetParamMat("input", input, points_are_rows)
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  # Call the program.
  preprocess_one_hot_encoding_mlpackMain()

  return IOGetParamMat("output", points_are_rows)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libsigcpp_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libsigcpp")
JLLWrappers.@generate_main_file("libsigcpp", UUID("8855df87-3ca1-5a3d-a68b-4f0ddc198ba8"))
end  # module libsigcpp_jll


text:
module PenaltyFunctions

using LinearAlgebra
using InteractiveUtils
using RecipesBase

# trait functions
include("traits.jl")

# penalty functions
include("common.jl")

# IO and plot recipes
include("printing.jl")
include("plotrecipes.jl")

export
    Penalty,
        ElementPenalty,
            ProxableElementPenalty,
                NoPenalty,
                L1Penalty,
                L2Penalty,
                ElasticNetPenalty,
            SCADPenalty,
            MCPPenalty,
            LogPenalty,
        ArrayPenalty,
            NuclearNormPenalty,
            GroupLassoPenalty,
            MahalanobisPenalty,
    addgrad, prox, prox!, deriv, value, grad, grad!, addgrad!, scaled

end


text:
""" given two generators `g1`, `g2`, `vertices` of their common boundary
(in generator representation) and the list of all generators, compute the boundary volume.
It works by constructing a polytope from halfspaces between the generators of the boundary
vertices. We use an affine transformation to get rid of the dimension in direction g1->g2
"""
function boundary_area(g1::Int, g2::Int, vertices::AbstractVector{<:Sigma}, generators)
    A = generators[g1]
    B = generators[g2]
    transform = transformation(A, B)
    A = transform(A)
    gen_inds = unique!(sort!(reduce(vcat,vertices)))

    halfspaces = []
    for gen_ind in gen_inds
        gen_ind in (g1,g2) && continue
        B = transform(generators[gen_ind])
        u = normalize(B-A)
        b = dot(u, (A+B)/2)
        hf = HalfSpace(u[2:end], b)
        push!(halfspaces, hf)
    end
    halfspaces = [h for h in halfspaces]

    # Performance with different libraries on random 6x100 data:
    # nolib: 30 mins
    # CDDLib: 40 mins
    # QHull: 5h (with warnings about not affine polyhedron using a solver)
    poly = polyhedron(hrep(halfspaces))
    local vol
    try
        vol = hasrays(poly) ? Inf : volume(poly)
    catch e
        if isa(e, AssertionError)
            @warn e
            vol = NaN
        else
            rethrow(e)
        end
    end

    return vol
end

""" affine transformation rotatinig and translating such that the boundary is aligned with
the first dimension. A->B will be mapped to const*[1,0,0,...] and (A+B)/2 to [0,0,...] """
function transformation(A, B)
    R = diagm(ones(length(A)))
    R[:,1] = B-A
    R = collect(inv(qr(R).Q))  # Shortcut for Gram-Schmidt orthogonalization.
    t = R*(A+B)/2
    transform(x) = R*x - t
    return transform
end

function transformation(A::SVector{N, T}, B::SVector{N, T}) where {N, T}
    R = MMatrix{N,N,T}(I)
    R[:,1] = B-A
    R = inv(qr(R).Q)
    t = R * (A+B)/2
    transform(x) = R*x - t
end

""" build the connectivity matrix for the SQRA from adjacency and boundary information """
function volumes(vertices, P::AbstractVector)
    dim = length(P[1])
    conns = adjacency(vertices)

    I = Int[]
    J = Int[]
    As = Float64[]
    Vs = zeros(length(P))

    #=@showprogress 1 "Voronoi volumes "=# for ((g1,g2), sigs) in conns
        a = boundary_area(g1, g2, sigs, P)
        h = norm(P[g1] - P[g2]) / 2
        v = a * h / dim  # Volume computation
        push!(I, g1)
        push!(J, g2)
        push!(As, a)
        Vs[g1] += v
        Vs[g2] += v
    end
    A = sparse(I, J, As, length(P), length(P))
    A = A + A'

    return A, Vs
end

""" given vertices in generator-coordinates,
collect the verts belonging to generator pairs, i.e. boundary vertices """
function adjacency(v::Vertices)
    conns = Dict{Tuple{Int,Int}, Vector{Vector{Int}}}()

    for (sig, r) in v
        for a in sig
            for b in sig
                a >= b && continue
                v = get!(conns, (a,b), [])
                push!(v, sig)
            end
        end
    end

    return conns
end

""" similar to `boundary_area`, however uses a vector representation and is slower """
function boundary_area_vrep(g1::Int, g2::Int, inds::AbstractVector{<:Sigma}, vertices::Vertices, generators)
    A = generators[g1]
    B = generators[g2]
    dim = length(A)
    vertex_coords = map(i->vertices[i], inds)
    push!(vertex_coords, A)  # Append one Voronoi center for full-dimensional volume.

    poly = polyhedron(vrep(vertex_coords), QHull.Library())
    vol = hasrays(poly) ? Inf : volume(poly)

    h = norm(B - A) / 2
    A = dim * vol / h
    return A, h, vol
end


text:
#   This file is part of Similitude.jl
#   It is licensed under the MIT license
#   Similitude Copyright (C) 2022 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

abstract type AbstractModule <: Number end

abstract type AbelianGroup <: AbstractModule end

struct Group{T,N,S} <: AbelianGroup
    v::Values{N,T}
    c::S
    @pure function Group{T,N,S}(v,c) where {N,T,S}
        r = promoteint(UnitSystems.cache(c))
        new{T,N,typeof(r)}(v,r)
    end
    @pure function Group{T,N,S}(v,c) where {N,T<:Rational,S}
        r,p = promoteints(v),promoteint(UnitSystems.cache(c))
        new{eltype(r),N,typeof(p)}(r,p)
    end
    @pure function Group{T,N,S}(v::Values{T,N},c) where {N,T<:Rational,S}
        r,p = promoteints(v),promoteint(UnitSystems.cache(c))
        new{eltype(r),N,typeof(p)}(r,p)
    end
end

@pure Group{T,N}(v,c::S=1) where {T,N,S} = Group{T,N,S}(v,c)
@pure Group{T}(v::Values{N,T},c::S=1) where {N,T,S} = Group{T,N,S}(v,c)
@pure Group(v::Values{N,T},c::S=1) where {N,T,S} = Group{T,N,S}(v,c)

@pure function Group(v::Values{N,<:Rational},c::S=1) where {N,S}
    r = promoteints(v)
    Group{eltype(r),N,S}(r,c)
end

import AbstractTensors: value
value(g::Group) = g.v
coef(g::Group) = UnitSystems.measure(g.c)

Base.:(==)(a::Group,b::Group) = a.v == b.v && a.c == b.c

Base.abs(a::Group{T,N,S}) where {T,N,S} = Group{T,N,S}(a.v,abs(a.c))

@pure isonezero(x) = isone(x) || iszero(x)
@pure checkint(v::T) where T<:Integer = v
@pure checkint(v::T) where T<:Rational = isone(denominator(v))
@pure checkint(v) = isonezero(v)
@pure checkints(v::Values{N} where N) = prod(isonezero.(v.v))
@pure checkints(v::Values{N,<:Rational} where N) = prod(isone.(denominator.(v.v)))
@pure checkints(v::Values{N,<:Integer} where N) = v
@pure promoteint(v::T) where T<:Integer = v
@pure promoteint(v) = checkint(v) ? Int(v) : v
@pure promoteint(v::UnitSystems.Constant) = isone(v) ? 1 : v
@pure promoteints(v::Values{N,<:Integer} where N) = v
@pure promoteints(v) = checkints(v) ? Int.(v) : v

const expos = Values('⁰','¹','²','³','⁴','⁵','⁶','⁷','⁸','⁹')
const chars = Dict([[string(i-1)[1]=>expos[i] for i ∈ 1:length(expos)];['.'=>'⋅','-'=>'⁻','e'=>'ᵉ','v'=>'ᵛ','₀'=>'⁰','₁'=>'¹','₂'=>'²','₃'=>'³','₄'=>'⁴','₅'=>'⁵','₆'=>'⁶','₇'=>'⁷','₈'=>'⁸','₉'=>'⁹','*'=>'*','⋅'=>'⋅']])

makeint(x) = x
@pure makint(x::Int) = x
@pure function makeint(x::AbstractFloat)
    ax,rem,ne = abs(x),abs(x%1),sqrt(eps()*norm(x))
    if ne < 1
        if log10(ax)-log(1.7,rem) > 20
            return Int(x÷1)
        elseif log10(ax)-log10(1-rem) > 17
            return Int(x÷1)+1
        else
            return x
        end
    else
        return x
    end
end

findpower(x) = 0
findpower(x::Int,i=Int(round(log10(x),RoundDown))) = i<0 ? 0 : (d = (x÷(10^i))%10; !iszero(d) ? findpower(x,i-1) : i+1)

printexpo(io::IO, d, x) = !iszero(x) && (print(io, d); printexpo(io,x))
function printexpo(io::IO, d, x::AbstractFloat)
    if !iszero(x)
        ix = makeint(x)
        if abs(x) < 1
            mix = makeint(inv(x))
            if typeof(mix) == Int
                printexpo(io, d, 1//mix)
            else
                if (d == 10 || d == "10") && length(string(abs(x)))>5
                    x < 0 && print(io, '/')
                    print(io, makeint(10^abs(x)))
                    !(x<0) && print(io, '⋅')
                else
                    print(io, d)
                    printexpo(io, x)
                end
            end
        elseif d == 10 || d == "10"
            mx = makeint(abs(x))
            x < 0 && print(io, '/')
            if typeof(mx)== Int
                printexpo(io, d, mx)
            else
                ten = makeint(10^abs(x))
                pow = findpower(ten)
                if !iszero(pow)
                    net = ten÷10^pow
                    if !isone(net)
                        print(io, net)
                        print(io, x < 0 ? '/' : '⋅')
                    end
                    print(io, d)
                    printexpo(io, pow)
                elseif length(string(abs(x)))>5
                    print(io, ten)
                    !(x<0) && print(io, '⋅')
                else
                    printexpo(io, d, rationalize(x))
                end
                if !iszero(pow)

                end
            end
        elseif typeof(ix) == Int
            printexpo(io,d,ix)
        else
            print(io, d)
            printexpo(io,x)
        end
    end
end

function printexpo(io::IO, x::Integer)
    if !isone(x)
        print(io, (x<0 ? ('⁻',) : ())..., expos[reverse(digits(abs(x)).+1)]...)
    end
end

function printexpo(io::IO, x::AbstractFloat)
    if !isone(x)
        print(io, (x<0 ? ('⁻',) : ())..., [chars[i] for i ∈ string(abs(x))]...)
    end
end

function printexpo(io::IO, x::Rational)
    if !isone(x)
        print(io, (x.num<0 ? ('⁻',) : ())..., expos[reverse(digits(abs(x.num)).+1)]...)
        !isone(x.den) && print(io, 'ᐟ', expos[reverse(digits(x.den).+1)]...)
    end
end

function printexpo(io::IO, x::Complex)
    if !isone(x)
        if !iszero(x.re)
            isone(x.re) ? print(io, '¹') : printexpo(io, x.re)
            !iszero(x.im) && print(io, x.im<0 ? '⁻' : '⁺')
        else
            x.im<0 && print(io, '⁻')
        end
        if !iszero(x.im)
            printexpo(io, abs(x.im))
            print(io, "ⁱᵐ")
        end
    end
end

function printdims(io::IO,x::Group{T,N},name) where {T,N}
    M = 0
    if haskey(ENV,"GROUPAREN")
        for i ∈ 1:N-M
            printnum(io, name[i], x.v[i])
        end
        n = sum(first(x.v,N-M).<0)
        n>0 && print(io, '/')
        n>1 && print(io, '(')
        for i ∈ 1:N-M
            printnum(io, name[i], -x.v[i])
            typeof(name[i])==String && isone(x.v[i]) && !iszero(norm(last(first(x.v,N-M),N-i-1))) && print(io,'*')
        end
        n>1 && print(io, ')')
    else
        for i ∈ 1:N-M
            printexpo(io, name[i], makeint(x.v[i]))
            typeof(name[i])==String && isone(x.v[i]) && !iszero(norm(last(first(x.v,N-M),N-i-M))) && print(io,'⋅')
        end
    end
end

function printnum(io, b, e)
    me = makeint(e)
    !(me<0) && printexpo(io, b, me)
end

Base.show(io::IO,x::Group) = showgroup(io,x)
function showgroup(io::IO,x::Group{T,N},u=Natural,c='𝟙') where {T,N}
    #back = T<:AbstractFloat && x.v[N]<0
    #!back && printexpo(io, 10, x.v[N])
    printdims(io,x,u)
    iz = iszero(norm(x.v))
    xc = coef(x)
    iz && (isone(xc)||abs(xc)<1) && print(io, c)
    #back && printexpo(io, 10, last(x.v))
    if !isone(xc)
        if abs(xc)<1 && !isconstant(xc)
            print(io,'/',makeint(inv(xc)))
        else
            !iz && print(io, '⋅')
            if isconstant(xc)
                print(io, '(', makeint(xc), ')')
            else
                print(io, makeint(xc))
            end
        end
    end
end

function showgroup(io::IO, x::AbelianGroup, u=Metric, c='𝟙')
    showfun(io,x)
    showgroup(io,dimensions(x),u,c)
    print(io,')')
end

# log

struct LogGroup{B,T<:AbelianGroup} <: AbelianGroup
    v::T
    @pure LogGroup{B,T}(v) where {B,T} = new{B,T}(v)
end

@pure LogGroup{B}(v::T) where {B,T} = LogGroup{B,T}(v)
@pure LogGroup(d::AbelianGroup) = LogGroup{ℯ}(d)

@pure islog(::LogGroup) = true
@pure islog(x) = false

@pure base(x::LogGroup{B}) where B = B
@pure dimensions(x::LogGroup) = x.v

value(g::LogGroup) = g.v

Base.show(io::IO, x::LogGroup) = (showfun(io,x); print(io,dimensions(x),')'))

showfun(io::IO, x::LogGroup{B}) where B = print(io,"log(",B,',')
showfun(io::IO, x::LogGroup{ℯ}) = print(io,"log(")
showfun(io::IO, x::LogGroup{2}) = print(io,"log2(")
showfun(io::IO, x::LogGroup{10}) = print(io,"log10(")
showfun(io::IO, x::LogGroup{exp10(0.1)}) = print(io,"dB(")

Base.log(x::AbelianGroup) = LogGroup(x)
Base.log2(x::AbelianGroup) = LogGroup{2}(x)
Base.log10(x::AbelianGroup) = LogGroup{10}(x)
Base.log(b::Number,x::AbelianGroup) = LogGroup{b}(x)
Base.exp(x::LogGroup{ℯ}) = dimensions(x)
Base.exp2(x::LogGroup{2}) = dimensions(x)
Base.exp10(x::LogGroup{10}) = dimensions(x)
Base.exp(x::LogGroup{B}) where B = dimensions(x)^inv(log(B))
Base.exp2(x::LogGroup) = exp2(x)
Base.exp10(x::LogGroup) = exp10(x)
Base.:^(b::T,x::LogGroup) where T<:Number = exp(x*log(b))

Base.:+(x::LogGroup{B},y::LogGroup{B}) where B = LogGroup{B}(x.v*y.v)
Base.:-(x::LogGroup{B},y::LogGroup{B}) where B = LogGroup{B}(x.v/y.v)
Base.:/(x::LogGroup{B},y::T) where {B,T<:Number} = LogGroup{B^y}(x.v)
Base.:*(x::LogGroup,y::T) where T<:Number = x/inv(y)
Base.:*(x::T,y::LogGroup) where T<:Number = y*x

# exp

struct ExpGroup{B,T<:AbelianGroup} <: AbelianGroup
    v::T
end

@pure ExpGroup{B}(v::T) where {B,T} = ExpGroup{B,T}(v)
@pure ExpGroup(d::AbelianGroup) = ExpGroup{ℯ}(d)

@pure base(x::ExpGroup{B}) where B = B
@pure dimensions(x::ExpGroup) = x.v

value(g::ExpGroup) = g.v

Base.show(io::IO, x::ExpGroup) = (showfun(io,x); print(io,dimensions(x),')'))

showfun(io::IO, x::ExpGroup{B}) where B = print(io,B,"^(")
showfun(io::IO, x::ExpGroup{ℯ}) = print(io,"exp(")
showfun(io::IO, x::ExpGroup{2}) = print(io,"exp2(")
showfun(io::IO, x::ExpGroup{10}) = print(io,"exp10(")

Base.exp(x::AbelianGroup) = ExpGroup(x)
Base.exp2(x::AbelianGroup) = ExpGroup{2}(x)
Base.exp10(x::AbelianGroup) = ExpGroup{10}(x)
Base.:^(b::T,x::AbelianGroup) where T<:Number = ExpGroup{b}(x)
Base.log(x::ExpGroup{ℯ}) = dimensions(x)
Base.log2(x::ExpGroup{2}) = dimensions(x)
Base.log10(x::ExpGroup{10}) = dimensions(x)
Base.log(b::Number,x::ExpGroup{B}) where B = dimensions(x)/log(B,b)
Base.log(x::ExpGroup) = log(ℯ,x)
Base.log2(x::ExpGroup) = log(2,x)
Base.log10(x::ExpGroup) = log(10,x)

Base.:^(x::ExpGroup{B},y::T) where {B,T<:Number} = iszero(y) ? one(x) : x
Base.:^(x::ExpGroup{B},y::T) where {B,T<:Integer} = iszero(y) ? one(x) : x
Base.:^(x::ExpGroup,y::ExpGroup) = ExpGroup{x}(y)
Base.:*(x::ExpGroup{B},y::ExpGroup{B}) where B = ExpGroup{B}(x.v+y.v)
Base.:/(x::ExpGroup{B},y::ExpGroup{B}) where B = ExpGroup{B}(x.v-y.v)

Base.:*(x::ExpGroup{X},y::ExpGroup{Y}) where {X,Y} = ExpGroup{X*Y}(x.v+y.v)

# other

logdb(x::AbelianGroup) = log(exp10(0.1),x)

Base.one(::AbelianGroup) = Group(zeros(Values{dims,Int}))
Base.isone(x::Group) = iszero(norm(x.v)) && isone(x.c)
Base.isone(x::AbelianGroup) = false
Base.zero(x::AbelianGroup) = log(one(x))
Base.iszero(x::AbelianGroup) = false
Base.iszero(x::LogGroup) = isone(dimensions(x))

coefprod(a,b) = a*b
coefprod(a::UnitSystems.Constant,b) = a*Constant(b)
coefprod(a,b::UnitSystems.Constant) = Constant(a)*b
coefprod(a::UnitSystems.Constant,b::UnitSystems.Constant) = a*b
Base.:*(a::Group,b::Group) = Group(a.v+b.v,coefprod(coef(a),coef(b)))
Base.:/(a::Group,b::Group) = Group(a.v-b.v,coefprod(coef(a),inv(coef(b))))
Base.:^(a::Group,b::Number) = Group(b*a.v,coef(a)^b)
Base.:^(a::Group,b::Integer) = Group(b*a.v,coef(a)^b)
Base.:^(a::Group,b::Rational) = Group(b*a.v,coef(a)^b)
Base.sqrt(a::Group{T}) where T = Group{T}(a.v/2,sqrt(coef(a)))
Base.sqrt(a::Group{Int}) = Group{Rational{Int}}(a.v//2,sqrt(coef(a)))
Base.cbrt(a::Group{T}) where T = Group{T}(a.v/3,cbrt(coef(a)))
Base.cbrt(a::Group{Int}) = Group{Rational{Int}}(a.v//3,cbrt(coef(a)))
Base.inv(a::Group{T}) where T = Group{T}(-a.v,inv(coef(a)))

@pure valueat(::Val{j},::Val{k}) where {j,k} = valueat(j,k)
valueat(j,k,z::T=1) where T = Group(Values{k,T}([i==j ? z : 0 for i ∈ 1:k]))
Base.:*(a::Number,b::Group{T,N}) where {T,N} = Group{T,N}(b.v,coefprod(a,coef(b)))
Base.:*(a::Group{T,N},b::Number) where {T,N} = Group{T,N}(a.v,coefprod(coef(a),b))
Base.:/(a::Number,b::Group) = a*inv(b)#
Base.:/(a::Group,b::Number) = a*inv(b)#


text:
# Julia translation of http://nbviewer.jupyter.org/github/barbagroup/AeroPython/blob/master/lessons/01_Lesson01_sourceSink.ipynb
# Lession 2 Source and Sink in a Freestream

using PyPlot
using Distributions

close("all")
meshgrid(x,y) = (repmat(x',length(y),1),repmat(y,1,length(x)))

N = 200                                # number of points in each direction
x_start, x_end = -4.0, 4.0            # boundaries in the x-direction
y_start, y_end = -2.0, 2.0            # boundaries in the y-direction
x = linspace(x_start, x_end, N)    # creates a 1D-array with the x-coordinates
y = linspace(y_start, y_end, N)    # creates a 1D-array with the y-coordinates

X,Y=meshgrid(x,y)

u_inf = 1.0

# computes the freestream velocity field
u_freestream = u_inf * ones((N, N))
v_freestream = zeros((N, N))

# computes the stream-function
psi_freestream = u_inf * Y

function get_velocity(strength, xs, ys, X, Y):
    """Returns the velocity field generated by a source/sink.

    Arguments
    ---------
    strength -- strength of the source/sink.
    xs, ys -- coordinates of the source/sink.
    X, Y -- mesh grid.
    """
    u = strength/(2*pi).*(X-xs)./((X-xs).^2+(Y-ys).^2)
    v = strength/(2*pi).*(Y-ys)./((X-xs).^2+(Y-ys).^2)

    return u, v
  end

function get_stream_function(strength, xs, ys, X, Y):
      """Returns the stream-function generated by a source/sink.

      Arguments
      ---------
      strength -- strength of the source/sink.
      xs, ys -- coordinates of the source/sink.
      X, Y -- mesh grid.
      """
      psi = strength/(2*pi).*atan2((Y-ys), (X-xs))

      return psi
    end

strength_source = 5.0            # strength of the source
x_source, y_source = -1.0, 0.0   # location of the source

# computes the velocity field
u_source, v_source = get_velocity(strength_source, x_source, y_source, X, Y)

# computes the stream-function
psi_source = get_stream_function(strength_source, x_source, y_source, X, Y)

  # superposition of the source on the freestream
  u = u_freestream + u_source
  v = v_freestream + v_source
  psi = psi_freestream + psi_source

  # plots the streamlines

  size = 10
  PyPlot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
  PyPlot.grid(true)
  PyPlot.xlabel("x", fontsize=16)
  PyPlot.ylabel("y", fontsize=16)
  PyPlot.xlim(x_start, x_end)
  PyPlot.ylim(y_start, y_end)
  PyPlot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle="->")
  PyPlot.scatter(x_source, y_source, color="#CD2305", s=80, marker="o")

  # calculates the stagnation point
  x_stagnation = x_source - strength_source/(2*pi*u_inf)
  y_stagnation = y_source

  # adding the stagnation point to the figure
  PyPlot.scatter(x_stagnation, y_stagnation, color="g", s=80, marker="o")

  # adds the dividing line to the figure
  PyPlot.contour(X, Y, psi,
              levels=[-strength_source/2, +strength_source/2],
              colors="#CD2305", linewidths=2, linestyles="solid");

strength_sink = -5.0        # strength of the sink
x_sink, y_sink = 1.0, 0.0   # location of the sink

# computes the velocity field on the mesh grid
u_sink, v_sink = get_velocity(strength_sink, x_sink, y_sink, X, Y)

# computes the stream-function on the grid mesh
psi_sink = get_stream_function(strength_sink, x_sink, y_sink, X, Y)

# superposition of a source and a sink on the freestream
u = u_freestream + u_source + u_sink
v = v_freestream + v_source + v_sink
psi = psi_freestream + psi_source + psi_sink

# plots the streamlines
size = 10
PyPlot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
PyPlot.xlabel("x", fontsize=16)
PyPlot.ylabel("y", fontsize=16)
PyPlot.xlim(x_start, x_end)
PyPlot.ylim(y_start, y_end)
PyPlot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle="->")
PyPlot.scatter([x_source, x_sink], [y_source, y_sink], color="#CD2305", s=80, marker="o")
PyPlot.contour(X, Y, psi, levels=[0.], colors="#CD2305", linewidths=2, linestyles="solid");

# computes the pressure coefficient field
cp = 1.0 - (u.^2+v.^2)./u_inf^2

# plots the pressure coefficient field
size = 10
PyPlot.figure(figsize=(1.1*size, (y_end-y_start)/(x_end-x_start)*size))
PyPlot.xlabel("x", fontsize=16)
PyPlot.ylabel("y", fontsize=16)
PyPlot.xlim(x_start, x_end)
PyPlot.ylim(y_start, y_end)
contf = PyPlot.contourf(X, Y, cp, levels=linspace(-2.0, 1.0, 100), extend="both")
cbar = PyPlot.colorbar(contf)
cbar[:set_label](L"$C_p$", fontsize=16)
cbar[:set_ticks]([-2.0, -1.0, 0.0, 1.0])
PyPlot.scatter([x_source, x_sink], [y_source, y_sink], color="#CD2305", s=80, marker="o")
PyPlot.contour(X, Y, psi, levels=[0.], colors="#CD2305", linewidths=2, linestyles="solid");


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using PkgBenchmark
using Oceananigans
using Benchmarks

baseline = BenchmarkConfig(id="main")
script = joinpath(@__DIR__, "benchmarkable_nonhydrostatic_model.jl")
resultfile = joinpath(@__DIR__, "regression_benchmarks.json")

print_system_info()

judgement = judge(Oceananigans, baseline, script=script, resultfile=resultfile, verbose=true)
results = PkgBenchmark.benchmarkgroup(judgement)

for (case, trial) in results
    println("Results for $case")
    display(trial)
end


text:
@json struct Release
    tag_name::String
    description::String
end

@json struct Tag
    commit::Commit
    release::Release
    name::String
    target::String
    message::String
end

endpoint(::GitLabAPI, ::typeof(get_tags), project::Integer) =
    Endpoint(:GET, "/projects/$project/repository/tags")
endpoint(::GitLabAPI, ::typeof(get_tags), owner::AStr, repo::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner, repo))/repository/tags")
into(::GitLabAPI, ::typeof(get_tags)) = Vector{Tag}


text:
# basic one
Y = [
     0  3;
     1  0;
     1 -3;
    -2  3;
     0  0;
]

D, X = ksvd(Y, 8, max_iter_mp = 600, sparsity_allowance = 1.0)
@test norm(Y-D*X) < 1e-4

# sparsity_allowance must be > 0
@test_throws ArgumentError ksvd(Y, 20, sparsity_allowance = -0.1)
@test_throws ArgumentError ksvd(Y, 20, sparsity_allowance = 1.1)

# should work normally
ksvd(Y, 20, max_iter = 1)
ksvd(Y, 20, sparsity_allowance = 0.0)
ksvd(Y, 20, sparsity_allowance = 1.0)

# But should work well when size(Y, 1) == n_atoms
Y = [
    -1 1 2;
     1 0 1
]
D, X = ksvd(Y, 2, max_iter_mp = 4000)
@test norm(Y-D*X) < 0.001  # relax the constraint since the dictionary is small

# Return only if X is sparse enough
Y = [
	0  2  3 -1  1;
	1 -3  1  3  0
]

# More than 20% of elements in X must be zeros
sparsity_allowance = 0.2
D, X = ksvd(Y, 5, max_iter = Int(1e10), sparsity_allowance = sparsity_allowance)
@test sum(iszero, X) / length(X) > sparsity_allowance


text:
module RobustAmpliconDenoising

	include("include_all.jl")

	export

	denoise,
	denoise_pipeline,
	alignment_consensus,
	get_centroid,
	consensus_seq,
	refine_ref,
	consensus_viz,
	disagreements,
	diff_in_homopolymer_region,
	get_coarse_centroid,
	trim_to_refs

end


text:
# bring in the previous definitions
include("rocketparts.jl")
include("rigidbody.jl")
# this is slightly different implementation than rocketcoll.jl
# so we start from just the parts

# copied code that follows from rocketcoll.jl
# since we need it and it is not in rocketparts.jl
struct NoPayload <: Payload
end

mass(payload::NoPayload) = 0.0

Rocket(tank::Tank, engine::Engine) = Rocket(NoPayload(), tank, engine)

const nopayload = NoPayload()

# new definition of SpaceVehicle with rigid body and boolean to toggle gravity
mutable struct SpaceVehicle
    active_stage::Payload
    body::RigidBody
    gravity::Bool
end

# new constructor
function SpaceVehicle(rockets::Array{Rocket})
    first_stage = nopayload
    for r in reverse(rockets)
        if first_stage != nopayload
            r.payload = first_stage
        end
        first_stage = r
    end
    body = RigidBody(mass(first_stage), 0.0)
    SpaceVehicle(first_stage, body, true)
end

SpaceVehicle(rockets::Rocket...) = SpaceVehicle(collect(rockets))

# also copied the following from rocketcoll.jl
function stage_separate!(ship::SpaceVehicle)
    stage = ship.active_stage
    if stage isa Rocket
       ship.active_stage = stage.payload
       stage.payload
    else
       nothing
    end
 end
# thrust functions build upon each other
thrust(engine::SingleEngine)   = engine.thrust
thrust(cluster::EngineCluster) = thrust(cluster.engine) * cluster.count
thrust(ship::SpaceVehicle) = thrust(ship.active_stage)
thrust(payload::Payload) = 0.0
thrust(r::Rocket) = thrust(r.engine)

# mass functions also built upon each other. We get mass for
# higher level objects by getting the mass of lower level objects:
mass(tank::Tank) = tank.dry_mass + tank.propellant
mass(payload::NoPayload) = 0.0
mass(probe::SpaceProbe) = probe.mass
mass(capsule::Capsule) = capsule.mass
mass(satellite::Satellite) = satellite.mass

function mass(r::Rocket)
    mass(r.payload) + mass(r.tank) + mass(r.engine)
end

mass(ship::SpaceVehicle) = mass(ship.active_stage)

# The launch! function relies on knowing how much propellant is left
# to decide whether stage separation needs to be performed.
propellant(non_rocket) = 0.0
propellant(r::Rocket)          = r.tank.propellant
propellant(r::SpaceVehicle)    = propellant(r.active_stage)
# code used for the actual simulation

# the rotation function we defined in linearalgebra.jl
# copied it here so as not to include all the other baggage
function rotation(deg::Real)
    rad = deg2rad(deg)
    cosθ = cos(rad)
    sinθ = sin(rad)

    [cosθ -sinθ;
     sinθ  cosθ]
  end

  const earth_radius = 6.38e6
  const earth_mass   = 5.98e24
  const gravitation_constant = 6.673e-11

# Engheim has the annoying habit of using new things without explainging
# No problem that he wants to encourage you to look stuff up, but at least
# provide the link!
# Here is documentation on doc strings in Julia: https://docs.julialang.org/en/v1/manual/documentation/
"""
   earth_gravity(d)
The acceleration of gravity distance `d`
from earth surface. Should be `g₀` when `d == 0`.
"""
function earth_gravity(d::Real)
    (gravitation_constant * earth_mass)/(earth_radius + d)^2
end

function gravity_force(body::RigidBody)
  _, elevation = body.position
  [-earth_gravity(elevation) * body.mass, 0]
end

# note Engheim's signature which differentiates this update from the one called inside
# I used a different function name to make it less confusing#
#function update!(ship::SpaceVehicle, Δt::Number)
function simupdate!(ship::SpaceVehicle, Δt::Number)
    stage = ship.active_stage
    body  = ship.body

    # call the update function defined in rocketparts.jl
    update!(stage, Δt)

    # Engheim correctly points out it is better to copy the ship data to the rigid body
    # instead of including the ship in hte rigidbody for separation of concerns
    # rigidbody.jl is a separate module so can be reused elsewhere
    body.mass   = mass(ship)
    # Engheim: "We multiply a vector representing the force produced by the rocket engine
    # with a rotation matrix. This matrix is setup to rotate the force to the same
    # direction as the one the ship is pointing.”
    body.force  = rotation(body.orientation) * [thrust(ship), 0.0]
    # add the force of gravity if it is toggled as true
    if ship.gravity
        body.force += gravity_force(body)
    end

    integrate!(body, Δt)
end

# simulation code
# max duration is to ensure no infinite loops
function launch!(ship::SpaceVehicle, Δt::Number; max_duration::Number = 5000)
    t = 0 # start time
    while ship.active_stage isa Rocket
        while propellant(ship) > 0 && t <= max_duration
            simupdate!(ship, Δt)
            t += Δt
        end
        stage_separate!(ship)

        # Current stage cannot get us higher
        if thrust(ship) <= 0
            return ship
        end
    end
    ship
end

# helper functions
pos(ship::SpaceVehicle) = ship.body.position
velocity(ship::SpaceVehicle) = ship.body.velocity

# Example run
# merlin  = Engine(845e3, 282, 470)
# cluster = EngineCluster(merlin, 9)

# firststage  = Rocket(Tank(23.1e3, 418.8e3), cluster)
# secondstage = Rocket(SpaceProbe(22e3),
#                      Tank(4e3, 111.5e3),
#                      Engine(934e3, 348, 470))

# ship = SpaceVehicle(firststage, secondstage)
# Δt = 0.1
# launch!(ship, Δt)
# yay! it works
# julia> velocity(ship)
# 2-element Vector{Float64}:
#  3703.572857010847
#     0.0

# julia> pos(ship)
# 2-element Vector{Float64}:
#  937346.6102144092
#       0.0

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
###################################################################
#  Common definitions
###################################################################

type VegaLiteVis
    vis::Dict{}
end

+(a::VegaLiteVis, b::VegaLiteVis) = VegaLiteVis(softmerge(a.vis, b.vis))
*(a::VegaLiteVis, b::VegaLiteVis) = VegaLiteVis(merge(a.vis, b.vis))

function softmerge(a::Dict, b::Dict)
    ck = intersect(keys(a), keys(b))
    nd = merge(a, b)
    for k in ck
        nd[k] = softmerge(a[k], b[k])
    end
    nd
end
softmerge(a::Dict, b::Any) = a
softmerge(a::Any, b::Dict) = b
softmerge(a, b) = b

#  Switch for plotting in SVGs or canvas
SVG = true
svg() = SVG
svg(b::Bool) = (global SVG ; SVG = b)

#  Switch for showing or not the "save as PNG buttons"
SAVE_BUTTONS = true
buttons() = SAVE_BUTTONS
buttons(b::Bool) = (global SAVE_BUTTONS ; SAVE_BUTTONS = b)

# build Dict from kwargs, checks against signature
function _mkdict(signature, properties)
  d = Dict()
  for (p,v) in properties
    haskey(signature, p) || error("unknown property $p")
    isa(v, signature[p]) || error("property $p : expected $(signature[p]), got $(typeof(v))")
    d[p] = v
  end
  d
end

function _mkvis(pos::Tuple, signature, properties)
  d = _mkdict(signature, properties)
  for k in reverse(pos)
    d = Dict(k => d)
  end
  VegaLiteVis(d)
end


text:
# Abdessalem, A., Dervilis, N., Wagg, D. & Worden, K. Model selection and parameter estimation in structural dynamics using approximate Bayesian computation. Mechanical Systems and Signal Processing 99, 306–325 (2018).
@everywhere using DifferentialEquations

@everywhere function model1(ydot, y, p, t, sitp)
    m,c,k,k3 = p

    ydot[1] = y[2]
    ydot[2] = -(c/m)*y[2] - (k/m)*y[1] - (k3/m)*y[1]^3 + (1/m)*sitp[t]
end

@everywhere function model2(ydot, y, p, t, sitp)
    m,c,k,k3,k5 = p

    ydot[1] = y[2]
    ydot[2] = -(c/m)*y[2] - (k/m)*y[1] - (k3/m)*y[1]^3 - (k5/m)*y[1]^5 + (1/m)*sitp[t]
end

@everywhere function model_generator(identifier, scaled_interpolation)
    id = identifier + 1
    models = [model1, model2]
    m(ydot, y, p, t) = models[id](ydot, y, p, t, scaled_interpolation)
    return m
end

@everywhere function obj_func(dataset, model, params::Array{Float64, 1}, lbs, ubs)
    """
    :L: number of dataset
    :timepoints: timepoints of experimental data.
    :params: model parameters
    :eval_idx: indecies of components used for evaluating objective function value.
    """
    y0 = [0.0,0.0]
    params = (lbs.*(ubs./lbs).^(params./10)) # rescaling parameters
    alg = CVODE_BDF(max_hnil_warns = 0, max_error_test_failures = 0)

    L = 1
    rss_per_data = zeros(1)
    timepoints = linspace(0,5,500)
    tspan = (timepoints[1], timepoints[end])
    eval_idx = [1]
    # res = zeros(length(timepoints), length(eval_idx))
    for i in 1:L   # doseの数だけ回す。
        prob = ODEProblem(model, y0, tspan, params)   # rescaling parameters
        sol = solve(prob, alg, saveat=timepoints, maxiters=1e7,
                    abstol = 1e-8, reltol = 1e-8, verbose=false)
        if sol.retcode == :Success
            res = hcat(sol.u...)'
            rss_per_data[i] = mean((dataset[i][:, eval_idx] - res[:, eval_idx]).^2)
        else
            rss_per_data[i] = -1
        end
    end

    if -1 in rss_per_data
        return -1
    else
        rss = mean(rss_per_data)
        return sqrt(rss)
    end
end

@everywhere function test_modelgenerator(identifier)
    n_digits = 1
    identifier = bin(parse(Int,identifier), n_digits)
    id = Array{Int, 1}(zeros(n_digits))
    for i in 1:n_digits
        id[i] = parse(Int, identifier[i])
    end
    id = id .+ 1
    return id
end

@everywhere function test_all()
    ret = test_modelgenerator(ARGS[1])
    println(ret)
end


text:
"""
    synchronize!(x::NamedTuple; root_rank::Integer = 0)
    synchronize!(x::Tuple; root_rank::Integer = 0)
    synchronize!(x::AbstractArray; root_rank::Integer = 0)
    synchronize!(x::ComponentArray; root_rank::Integer = 0)
    synchronize!(x::Nothing; root_rank::Integer = 0)
    synchronize!(x::Symbol; root_rank::Integer = 0)
    synchronize!(x::Leaf; root_rank::Integer = 0)
    synchronize!(x::Number; root_rank::Integer = 0)

Synchronize `x` across all processes.
"""
function synchronize!(ps::Union{NamedTuple,Tuple}; root_rank::Integer=0)
    return fmap(x -> synchronize!(x; root_rank), ps)
end

function synchronize!(x::AbstractArray{T}; root_rank::Integer=0) where {T<:Number}
    bcast!(x, root_rank, MPI.COMM_WORLD)
    return x
end

function synchronize!(x::ComponentArray; root_rank::Integer=0)
    d = bcast!(getdata(x), root_rank, MPI.COMM_WORLD)
    return ComponentArray(d, getaxes(x))
end

# Ideally these things should be Tuples and not arrays
function synchronize!(x::AbstractArray; root_rank::Integer=0)
    synchronize!.(x; root_rank)
end

function synchronize!(l::Leaf; root_rank::Integer=0)
    @set! l.state = synchronize!(l.state; root_rank)
end

function synchronize!(x::Number; root_rank::Integer=0)
    return bcast!([x], root_rank, MPI.COMM_WORLD)[1]
end

# If we don't know what to synchronize, we don't do it
# -- Symbols, Nothing, Missing, Val, etc.
synchronize!(x; kwargs...) = x


text:
"""
    impRelax_f(f::Function,h::Function,hj::Function,X::Vector{Interval{T}},
               P::Vector{Interval{T}},p::Vector{T},pmid::Vector{T},
               mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}})

Relaxes the function `f(x,p)` by relaxation the state variable `x` using the implicit
function determined by `h(x,p)` with `x` in `X` and `p` in `P`. The reference
point for the affine relaxations is `pmid`. The parameters generated from the
relaxation at `pmid` are `param` and the basic parameters of the fixed point
method are `mc_opt`.
"""
function impRelax_f(f::Function,h::Function,hj::Function,X::Vector{Interval{T}},
                    P::Vector{Interval{T}},p::Vector{T},pmid::Vector{T},
                    mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}}) where {N,T<:AbstractFloat}
  np::Int64 = length(P)
  sone::SVector{np,T} = @SVector ones(np)
  p_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  xpMC::Vector{SMCg{np,T}} = MC_impRelax(h,hj,p_mc,pmid,X,P,mc_opt,param)
  return f(xpMC,p_mc)
end

"""
    impRelax_fg(f::Function,g::Function,h::Function,hj::Function,
                X::Vector{Interval{T}},P::Vector{Interval{T}},
                p::Vector{T},pmid::Vector{T},
                mc_opt::mc_opts,param::Vector{Vector{SMCg{N,T}}})

Relaxes the functions `f(x,p)` and `g(x,p)` by relaxation the state variable `x`
using the implicit function determined by `h(x,p)` with `x` in `X` and `p` in `P`.
The reference point for the affine relaxations is `pmid`. The parameters generated
from the relaxation at `pmid` are `param` and the basic parameters of the fixed point
method are `mc_opt`.
"""
function impRelax_fg(f::Function,g::Function,h::Function,hj::Function,
                     X::Vector{Interval{T}},P::Vector{Interval{T}},
                     p::Vector{T},pmid::Vector{T},
                     mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}}) where {N,T<:AbstractFloat}
  np::Int64 = length(P)
  sone::SVector{np,T} = @SVector ones(np)
  p_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  xpMC::Vector{SMCg{np,T}} = MC_impRelax(h,hj,p_mc,pmid,X,P,mc_opt,param)
  return f(xpMC,p_mc),g(xpMC,p_mc)
end

"""
    NimpRelax_f(f::Function,h::Function,hj::Function,X::Vector{Interval{T}},
               P::Vector{Interval{T}},p::Vector{T},pmid::Vector{T},
               mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}})

Relaxes the function `f(x,p)` by relaxation the state variable `x` using the implicit
function determined by `h(x,p)` with `x` in `X` and `p` in `P`. The reference
point for the affine relaxations is `pmid`. The parameters generated from the
relaxation at `pmid` are `param` and the basic parameters of the fixed point
method are `mc_opt`. Preconditioning is using a sparse LDU with full pivoting approach,
the jacobian with respect to X and the preconditioner are stored in a sparse format,
sparse contractors are applied, and h & hj are evaluated in place.
"""
function NimpRelax_f(f::Function,h!::Function,hj!::Function,X::Vector{Interval{T}},
                    P::Vector{Interval{T}},p::Vector{T},pmid::Vector{T},
                    mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}}) where {N,T<:AbstractFloat}
  np::Int64 = length(P)
  sone::SVector{np,T} = @SVector ones(np)
  p_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  xpMC::Vector{SMCg{np,T}} = MC_NimpRelax(h!,hj!,p_mc,pmid,X,P,mc_opt,param)
  return f(xpMC,p_mc)
end

"""
    NimpRelax_fg(f::Function,g::Function,h::Function,hj::Function,
                X::Vector{Interval{T}},P::Vector{Interval{T}},
                p::Vector{T},pmid::Vector{T},
                mc_opt::mc_opts,param::Vector{Vector{SMCg{N,T}}})

Relaxes the functions `f(x,p)` and `g(x,p)` by relaxation the state variable `x`
using the implicit function determined by `h(x,p)` with `x` in `X` and `p` in `P`.
The reference point for the affine relaxations is `pmid`. The parameters generated
from the relaxation at `pmid` are `param` and the basic parameters of the fixed point
method are `mc_opt`. Preconditioning is using a sparse LDU with full pivoting approach,
the jacobian with respect to X and the preconditioner are stored in a sparse format,
sparse contractors are applied, and h & hj are evaluated in place.
"""
function NimpRelax_fg(f::Function,g::Function,h!::Function,hj!::Function,
                     X::Vector{Interval{T}},P::Vector{Interval{T}},
                     p::Vector{T},pmid::Vector{T},
                     mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}}) where {N,T<:AbstractFloat}
  np::Int64 = length(P)
  sone::SVector{np,T} = @SVector ones(np)
  p_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  xpMC::Vector{SMCg{np,T}} = MC_NimpRelax(h!,hj!,p_mc,pmid,X,P,mc_opt,param)
  return f(xpMC,p_mc),g(xpMC,p_mc)
end

"""
    NdimpRelax_f(f::Function,h::Function,hj::Function,X::Vector{Interval{T}},
               P::Vector{Interval{T}},p::Vector{T},pmid::Vector{T},
               mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}})

Relaxes the function `f(x,p)` by relaxation the state variable `x` using the implicit
function determined by `h(x,p)` with `x` in `X` and `p` in `P`. The reference
point for the affine relaxations is `pmid`. The parameters generated from the
relaxation at `pmid` are `param` and the basic parameters of the fixed point
method are `mc_opt`. Preconditioning is using a LU with full pivoting approach,
the jacobian with respect to X and the preconditioner are stored, and h & hj are evaluated in place.
"""
function NdimpRelax_f(f::Function,h!::Function,hj!::Function,X::Vector{Interval{T}},
                    P::Vector{Interval{T}},p::Vector{T},pmid::Vector{T},
                    mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}}) where {N,T<:AbstractFloat}
  np::Int64 = length(P)
  sone::SVector{np,T} = @SVector ones(np)
  p_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  xpMC::Vector{SMCg{np,T}} = MC_NdimpRelax(h!,hj!,p_mc,pmid,X,P,mc_opt,param)
  return f(xpMC,p_mc)
end

"""
    NdimpRelax_fg(f::Function,g::Function,h::Function,hj::Function,
                X::Vector{Interval{T}},P::Vector{Interval{T}},
                p::Vector{T},pmid::Vector{T},
                mc_opt::mc_opts,param::Vector{Vector{SMCg{N,T}}})

Relaxes the functions `f(x,p)` and `g(x,p)` by relaxation the state variable `x`
using the implicit function determined by `h(x,p)` with `x` in `X` and `p` in `P`.
The reference point for the affine relaxations is `pmid`. The parameters generated from the
relaxation at `pmid` are `param` and the basic parameters of the fixed point
method are `mc_opt`. Preconditioning is using a LU with full pivoting approach,
the jacobian with respect to X and the preconditioner are stored, and h & hj are evaluated in place.
"""
function NdimpRelax_fg(f::Function,g::Function,h!::Function,hj!::Function,
                     X::Vector{Interval{T}},P::Vector{Interval{T}},
                     p::Vector{T},pmid::Vector{T},
                     mc_opt::mc_opts{T},param::Vector{Vector{SMCg{N,T}}}) where {N,T<:AbstractFloat}
  np::Int64 = length(P)
  sone::SVector{np,T} = @SVector ones(np)
  p_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  xpMC::Vector{SMCg{np,T}} = MC_NdimpRelax(h!,hj!,p_mc,pmid,X,P,mc_opt,param)
  return f(xpMC,p_mc),g(xpMC,p_mc)
end


text:
module Diffusion

export
  Problem,
  updatevars!,
  set_c!

using Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!

"""
    Problem(; parameters...)

Construct a constant diffusivity problem.
"""
function Problem(;
            nx = 128,
            Lx = 2π,
             κ = 0,
            dt = 0.01,
       stepper = "RK4",
             T = Float64,
           dev = CPU()
  )

      grid = OneDGrid(dev, nx, Lx; T=T)
    params = Params(dev, κ)
      vars = Vars(dev, grid)
  equation = Equation(dev, κ, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

struct Params{T} <: AbstractParams
  κ :: T
end

Params(dev, κ::Number) = Params(κ)
Params(dev, κ::AbstractArray) = Params(ArrayType(dev)(κ))

"""
    Equation(dev, κ, grid)

Returns the equation for a constant diffusivity problem on `grid` with diffusivity κ.
"""
function Equation(dev::Device, κ::T, grid) where T<:Number
  L = zeros(dev, T, grid.nkr)
  @. L = - κ * grid.kr^2

  return FourierFlows.Equation(L, calcN!, grid)
end

Equation(dev::Device, κ::T, grid::AbstractGrid{Tg}) where {T<:AbstractArray, Tg} =
  FourierFlows.Equation(0, calcN!, grid; dims=(grid.nkr,), T=cxtype(Tg))

struct Vars{Aphys, Atrans} <: AbstractVars
    c :: Aphys
   cx :: Aphys
   ch :: Atrans
  cxh :: Atrans
end

"""
    Vars(dev, grid)

Returns the variables `vars` for a constant diffusivity problem on `grid`.
"""
function Vars(::Dev, grid::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T grid.nx c cx
  @devzeros Dev Complex{T} grid.nkr ch cxh
  Vars(c, cx, ch, cxh)
end

"""
    calcN!(N, sol, t, clock, vars, params, grid)

Calculate the nonlinear term for the 1D heat equation.
"""
function calcN!(N, sol, t, clock, vars, params::Params{T}, grid) where T<:Number
  @. N = 0

  return nothing
end

function calcN!(N, sol, t, clock, vars, params::Params{T}, grid) where T<:AbstractArray
  @. vars.cxh = im * grid.kr * sol
  ldiv!(vars.cx, grid.rfftplan, vars.cxh)
  @. vars.cx *= params.κ
  mul!(vars.cxh, grid.rfftplan, vars.cx)
  @. N = im * grid.kr * vars.cxh

  return nothing
end

"""
    updatevars!(vars, grid, sol)

Update the variables in `vars` on the `grid` with the solution in `sol`.
"""
function updatevars!(vars, grid, sol)
  ldiv!(vars.c, grid.rfftplan, deepcopy(sol))
  @. vars.ch = sol

  return nothing
end

updatevars!(prob) = updatevars!(prob.vars, prob.grid, prob.sol)

"""
    set_c!(prob, c)

Set the solution as the transform of `c` and update `vars`.
"""
function set_c!(prob, c)
  T = typeof(prob.vars.c)
  prob.vars.c .= T(c) # this makes sure that c is converted to the ArrayType used in prob.vars.c (e.g., convert to CuArray if user gives c as Arrray)
  mul!(prob.sol, prob.grid.rfftplan, prob.vars.c)
  updatevars!(prob)

  return nothing
end

end # module


text:
module UpdateWindow

using Gtk, Gtk.ShortNames, Pkg
using ..Ahorn

function do_ask_dialog(text::String; no::String="Cancel", yes::String="Update")
    if ask_dialog(text, no, yes, Ahorn.window)
        info_dialog("Ahorn will now be updated.\nThis will freeze the window during the process.\nDo not close Ahorn, you will get notified when the update completes.", Ahorn.window)
        Pkg.update(["Maple", "Ahorn"])
        if ask_dialog("Ahorn updated!\nBut the old version of Ahorn still lingers around.\nWe can run gc for you to clean up installed package versions,\nbut that is going to affect all your julia environments.\nIf you do not have any other Julia packages installed,\nthis should be safe to do and save you some disk space.", Ahorn.window)
            Pkg.gc()
        end
        exit()
    end
end

# Modified from stream.jl
function take_stdout(f::Function)
    local oldout = Base.stdout
    local (rd, wr) = redirect_stdout()
    f()
    redirect_stdout(oldout)
    close(wr)
    local res = readavailable(rd)
    close(rd)
    return String(res)
end

pkgpath(pkg) = pathof(pkg) |> dirname |> dirname

function pkghash()
    try
        local ctx = Pkg.Types.Context()
        local pkg = PackageSpec(name="Ahorn", uuid=ctx.env.project.deps["Ahorn"])
        Pkg.Operations.resolve_versions!(ctx, [pkg])
        return string(pkg.repo.tree_sha)[1:7]
    catch e
        return nothing
    end
end

function updateAhorn(widget::Union{Ahorn.MenuItemsTypes, Nothing}=nothing)
    ask_dialog("""Would you like to check for updates?
            This will download files required for the update if there is one.
            The window might also freeze for a bit.
            Make sure not to close Ahorn during this!""",
            "Cancel", "Check for Updates", Ahorn.window) || return
    try
        sim = take_stdout() do
            # Simulate an update. This is still going to download the repo.
            Pkg.update(["Maple", "Ahorn"], preview=true)
        end

        if occursin(r"(Ahorn|Maple)[^\n]+⇒", sim) || occursin(r"~ (Ahorn|Maple)", sim)
            do_ask_dialog("A new version is available!\nDo you wish to update Ahorn?\nThis will close the program afterwards and you will have to rerun it.")
        else
            h = pkghash()
            do_ask_dialog("Ahorn seems to be up-to-date$(h !== nothing ? " at hash $h" : "").\nDo you wish to try updating anyway?\nThis will close the program afterwards and you will have to rerun it.", yes="Try Updating Anyway")
        end
    catch e
        do_ask_dialog("The update check failed for some reason.\nDo you wish to update Ahorn anyway?\nThis will close the program afterwards and you will have to rerun it.")
    end
end

end


text:
function PearsonsCorrelation()
    return PearsonsCorrelation(())
end

function PearsonsCorrelation(arg0::Covariance)
    return PearsonsCorrelation((Covariance,), arg0)
end

function PearsonsCorrelation(arg0::RealMatrix)
    return PearsonsCorrelation((RealMatrix,), arg0)
end

function PearsonsCorrelation(arg0::RealMatrix, arg1::jint)
    return PearsonsCorrelation((RealMatrix, jint), arg0, arg1)
end

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

function compute_correlation_matrix(obj::PearsonsCorrelation, arg0::RealMatrix)
    return jcall(obj, "computeCorrelationMatrix", RealMatrix, (RealMatrix,), arg0)
end

function compute_correlation_matrix(obj::PearsonsCorrelation, arg0::Vector{Vector{jdouble}})
    return jcall(obj, "computeCorrelationMatrix", RealMatrix, (Vector{Vector{jdouble}},), arg0)
end

function correlation(obj::PearsonsCorrelation, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "correlation", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function covariance_to_correlation(obj::PearsonsCorrelation, arg0::RealMatrix)
    return jcall(obj, "covarianceToCorrelation", RealMatrix, (RealMatrix,), arg0)
end

function get_correlation_matrix(obj::PearsonsCorrelation)
    return jcall(obj, "getCorrelationMatrix", RealMatrix, ())
end

function get_correlation_p_values(obj::PearsonsCorrelation)
    return jcall(obj, "getCorrelationPValues", RealMatrix, ())
end

function get_correlation_standard_errors(obj::PearsonsCorrelation)
    return jcall(obj, "getCorrelationStandardErrors", RealMatrix, ())
end


text:
################################
## Non-scalar linear indexing ##
################################

@inline function Base.getindex(a::TupleVector{N}, ::Colon) where N
    _getindex(a::TupleVector, Val(N), :)
end

@generated function _getindex(a::TupleVector, s::Val{N}, ::Colon) where N
    exprs = [:(a[$i]) for i = 1:N]
    return quote
        @_inline_meta
        @inbounds return similar_type(a,s)(tuple($(exprs...)))
    end
end

@propagate_inbounds function Base.getindex(a::TupleVector, inds::TupleVector{N,Int}) where N
    _getindex(a, Val(N), inds)
end

@generated function _getindex(a::TupleVector, s::Val{N}, inds::TupleVector{N, Int}) where N
    exprs = [:(a[inds[$i]]) for i = 1:N]
    return quote
        Base.@_propagate_inbounds_meta
        similar_type(a, s)(tuple($(exprs...)))
    end
end

@inline function Base.setindex!(a::TupleVector{N}, v, ::Colon) where N
    _setindex!(a::TupleVector, v, Val(N), :)
    return v
end

@generated function _setindex!(a::TupleVector, v, ::Val{L}, ::Colon) where {L}
    exprs = [:(a[$i] = v) for i = 1:L]
    return quote
        @_inline_meta
        @inbounds $(Expr(:block, exprs...))
    end
end

@generated function _setindex!(a::TupleVector, v::AbstractVector, ::Val{L}, ::Colon) where {L}
    exprs = [:(a[$i] = v[$i]) for i = 1:L]
    return quote
        Base.@_propagate_inbounds_meta
        @boundscheck if length(v) != L
            throw(DimensionMismatch("tried to assign $(length(v))-element array to length-$L destination"))
        end
        @inbounds $(Expr(:block, exprs...))
    end
end

@generated function _setindex!(a::TupleVector, v::TupleVector{M}, ::Val{L}, ::Colon) where {M,L}
    exprs = [:(a[$i] = v[$i]) for i = 1:L]
    return quote
        Base.@_propagate_inbounds_meta
        @boundscheck if M != L
            throw(DimensionMismatch("tried to assign $M-element array to length-$L destination"))
        end
        $(Expr(:block, exprs...))
    end
end

@propagate_inbounds function Base.setindex!(a::TupleVector, v, inds::TupleVector{N,Int}) where N
    _setindex!(a, v, Val(N), inds)
    return v
end

@generated function _setindex!(a::TupleVector, v, ::Val{N}, inds::TupleVector{N,Int}) where N
    exprs = [:(a[inds[$i]] = v) for i = 1:N]
    return quote
        Base.@_propagate_inbounds_meta
        similar_type(a, s)(tuple($(exprs...)))
    end
end

@generated function _setindex!(a::TupleVector, v::AbstractVector, ::Val{N}, inds::TupleVector{N,Int}) where N
    exprs = [:(a[inds[$i]] = v[$i]) for i = 1:N]
    return quote
        Base.@_propagate_inbounds_meta
        @boundscheck if length(v) != $N
            throw(DimensionMismatch("tried to assign $(length(v))-element array to length-$N destination"))
        end
        $(Expr(:block, exprs...))
    end
end

@generated function _setindex!(a::TupleVector, v::TupleVector{M}, ::Val{N}, inds::TupleVector{N,Int}) where {N,M}
    exprs = [:(a[inds[$i]] = v[$i]) for i = 1:N]
    return quote
        Base.@_propagate_inbounds_meta
        @boundscheck if M != N
            throw(DimensionMismatch("tried to assign $M-element array to length-$N destination"))
        end
        $(Expr(:block, exprs...))
    end
end


text:
module TestNormiid

using Distributions
using TimeData
using Base.Test

import Econometrics

## reload(joinpath(Pkg.dir(), "Econometrics", "src", "Econometrics.jl"))

#############
## NormIID ##
#############

## constructors / display
##-----------------------

issubtype(Econometrics.NormIID, Econometrics.AbstrUnivarModel)

@test_throws Exception Econometrics.NormIID(0.3, -0.2)

mod = Econometrics.NormIID(0.3, 1.2)
mod = Econometrics.NormIID([0.3, 1.2])

println("\nTesting display output:")
display(mod)
Econometrics.description(mod)

## getParams
##----------

mod = Econometrics.NormIID(0.3, 1.2)
@test Econometrics.getParams(mod) == (0.3, 1.2)

## simulate
##---------

nSim = 5000
simVals = Econometrics.simulate(mod, nSim)

## estimate
##---------

normIID_fit = Econometrics.estimate(Econometrics.NormIID, simVals)

## compare parameters
[[Econometrics.getParams(mod)...]'
[Econometrics.getParams(normIID_fit)...]']

end


text:
include(joinpath("manifolds", "SN.jl"))
include(joinpath("manifolds", "SO3.jl"))
include(joinpath("manifolds", "RN.jl"))
include(joinpath("manifolds", "product_manifolds.jl"))

embdim(::Type{M}) where M <: Manifold = dim(M)
basedim(::Type{M}) where M <: Manifold = dim(M)

dim(::M) where M <: Manifold = dim(M)
embdim(::M) where M <: Manifold = embdim(M)
basedim(::M) where M <: Manifold = dim(M)

isglobal(::Chart) = true    # Charts are considered global by default
default_chart(::Type{M}) where M <: Manifold = Chart{1,M}()
default_chart(::M) where M <: Manifold = default_chart(M)

function chart_transition(p1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J}
    p_emb = chart_to_emb(p1, C1)
    p2 = emb_to_chart(p_emb, C2)
end

choose_chart_emb(pe, C::Chart) =
    isglobal(C) ? C : choose_chart_emb(chart_choice_coord_rep(C), pe, C)
choose_chart_chart(p, C::Chart) =
    isglobal(C) ? C : choose_chart_chart(chart_choice_coord_rep(C), p, C)

choose_chart_emb(::ChartRep, pe, C) = choose_chart_chart(ChartRep(), emb_to_chart(pe, C), C)
choose_chart_chart(::EmbRep, p, C) = choose_chart_emb(EmbRep(), chart_to_emb(p, C), C)

choose_chart(::EmbRep, pe, C) = choose_chart_emb(p, C)
choose_chart(::ChartRep, p, C) = choose_chart_chart(p, C)

### Tangent bundles
T(::Chart{I,M}) where {M<:Manifold,I} = Chart{I,T{M}}()
base(::Chart{I,T{M}}) where {M<:Manifold,I} = Chart{I,M}()

dim(::Type{TangentBundle{M}}) where M <: Manifold = 2*dim(M)
embdim(::Type{TangentBundle{M}}) where M <: Manifold = 2*embdim(M)
basedim(::Type{TangentBundle{M}}) where M <: Manifold = dim(M)
baseembdim(::Type{TangentBundle{M}}) where M <: Manifold = embdim(M)

function tangent_vector_emb_splitview(pe, ::Type{TM}) where TM <: TangentBundle
    m, m_tb = baseembdim(TM), embdim(TM)
    xe, ve = pe[SVector{m}(1:m)], pe[SVector{m}(m+1:m_tb)]
end

function tangent_vector_chart_splitview(p, ::Type{TM}) where TM <: TangentBundle
    m, m_tb = basedim(TM), dim(TM)
    x, v = p[SVector{m}(1:m)], p[SVector{m}(m+1:m_tb)]
end

function chart_to_emb_differential(x, v, C::Chart{I,M}) where {M,I}
    xe = chart_to_emb(x, C)
    ∂xe_∂x = chart_to_emb_jacobian(x, C)
    ve = ∂xe_∂x*v
    xe, ve
end

function chart_to_emb_differential(x, v, a, C::Chart{I,M}) where {M,I}
    xe = chart_to_emb(x, C)
    p = [x; v]
    w = [v; a]
    ∂pe_∂p = chart_to_emb_jacobian(p, T(C))
    we = ∂pe_∂p*w
    ve, ae = tangent_vector_emb_splitview(we, T{M})
    xe, ve, ae
end

function emb_to_chart_differential(xe, ve, C::Chart{I,M}) where {M,I}
    x = emb_to_chart(xe, C)
    ∂x_∂xe = emb_to_chart_jacobian(xe, C)
    v = ∂x_∂xe*ve
    x, v
end

function emb_to_chart_differential(xe, ve, ae, C::Chart{I,M}) where {M,I}
    x = emb_to_chart(xe, C)
    pe = [xe; ve]
    we = [ve; ae]
    ∂p_∂pe = emb_to_chart_jacobian(pe, T(C))
    w = ∂p_∂pe*we
    v, a = tangent_vector_chart_splitview(w, T{M})
    x, v, a
end

function emb_to_chart(pe, C::Chart{I,T{M}}) where {M,I}
    xe, ve = tangent_vector_emb_splitview(pe, T{M})
    x, v = emb_to_chart_differential(xe, ve, base(C))
    p = [x; v]
end

function chart_to_emb(p, C::Chart{I,T{M}}) where {M,I}
    x, v = tangent_vector_chart_splitview(p, T{M})
    xe, ve = chart_to_emb_differential(x, v, Chart{I,M}())
    [xe; ve]
end

emb_to_chart_jacobian(pe, C::Chart{I,M}) where {M,I} =
    ForwardDiff.jacobian(pe -> emb_to_chart(pe, C), pe)
chart_to_emb_jacobian(pe, C::Chart{I,M}) where {M,I} =
    ForwardDiff.jacobian(pe -> chart_to_emb(pe, C), pe)
chart_transition_jacobian(p1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J} =
    ForwardDiff.jacobian(p1 -> chart_transition(p1, C1, C2), p1)

function chart_transition_jacobian(p1, C1::Chart{I,T{M}}, C2::Chart{J,T{M}}) where {M,I,J}
    m = dim(M)
    x1, v1 = tangent_vector_chart_splitview(p1, T{M})
    ∂x2_∂x1 = chart_transition_jacobian(x1, base(C1), base(C2))
    Hx2_x1 =  chart_transition_hessian(x1, base(C1), base(C2))
    ∂v2_∂v1 = SMatrix{dim(T{M}),dim(T{M})}([∂x2_∂x1 @SMatrix zeros(m,m)
               sum(Hx2_x1[k,:,:]*v1[k] for k in 1:m) ∂x2_∂x1])
end

function chart_transition_differential(x1, v1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J}
    x2 = chart_transition(x1, C1, C2)
    ∂x2_∂x1 = chart_transition_jacobian(x1, C1, C2)
    v2 = ∂x2_∂x1*v1
    x2, v2
end

function chart_transition_differential(x1, v1, a1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J}
    x2 = chart_transition(x1, C1, C2)
    p1 = [x1; v1]
    w1 = [v1; a1]
    ∂p2_∂p1 = chart_transition_jacobian(p1, T(C1), T(C2))
    w2 = ∂p2_∂p1*w1
    v2, a2 = tangent_vector_chart_splitview(w2, T{M})
    x2, v2, a2
end

function chart_transition(p1, C1::Chart{I,T{M}}, C2::Chart{J,T{M}}) where {M,I,J}
    x1, v1 = tangent_vector_chart_splitview(p1, T{M})
    x2, v2 = chart_transition_differential(x1, v1, base(C1), base(C2))
    p2 = [x2; v2]
end

# Print H[i,:,:] to see symmetric matrices
# index convention: ∂p2k_∂p1i∂p1j -> H[k,i,j]
function chart_transition_hessian(p1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J}
    m, S = dim(M), eltype(p1)
    inds = static(1):static(m)
    SArray{Tuple{m,m,m},S}([ForwardDiff.hessian(p1 -> chart_transition(p1, C1, C2)[k], p1)[i,j]
        for k in inds, i in inds, j in inds])
end

# Frames and coframes
function global_frame_coeffs_chart_transition(p1, B1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J}
    ∂p2_∂p1 = chart_transition_jacobian(p1, C1, C2)
    B2 = B1*∂p2_∂p1'
end

function global_coframe_coeffs_chart_transition(p1, β1, C1::Chart{I,M},
        C2::Chart{J,M}) where {M,I,J}
    ∂p2_∂p1 = chart_transition_jacobian(p1, C1, C2)
    β2 = inv(∂p2_∂p1)'*β1
end


text:
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Architectures: AbstractArchitecture

using KernelAbstractions: NoneEvent

import Oceananigans.BoundaryConditions:
    apply_x_bcs!,
    apply_y_bcs!,
    apply_z_bcs!,
    apply_x_east_bc!,
    apply_x_west_bc!,
    apply_y_south_bc!,
    apply_y_north_bc!,
    apply_z_top_bc!,
    apply_z_bottom_bc!

# Bunch o' shortcuts for halo communication bcs
apply_x_bcs!(Gc, ::AbstractGrid, c, ::HaloCommunicationBC, ::HaloCommunicationBC, ::AbstractArchitecture, args...) = NoneEvent()
apply_y_bcs!(Gc, ::AbstractGrid, c, ::HaloCommunicationBC, ::HaloCommunicationBC, ::AbstractArchitecture, args...) = NoneEvent()
apply_z_bcs!(Gc, ::AbstractGrid, c, ::HaloCommunicationBC, ::HaloCommunicationBC, ::AbstractArchitecture, args...) = NoneEvent()

@inline apply_x_east_bc!(  Gc, loc, ::HaloCommunicationBC, args...) = nothing
@inline apply_x_west_bc!(  Gc, loc, ::HaloCommunicationBC, args...) = nothing
@inline apply_y_north_bc!( Gc, loc, ::HaloCommunicationBC, args...) = nothing
@inline apply_y_south_bc!( Gc, loc, ::HaloCommunicationBC, args...) = nothing
@inline apply_z_top_bc!(   Gc, loc, ::HaloCommunicationBC, args...) = nothing
@inline apply_z_bottom_bc!(Gc, loc, ::HaloCommunicationBC, args...) = nothing


text:
using stoosc
rc("font", family="arial")

## Gemeral settings
params = Dict(
    "I_max" => 0.3,
    "forcing_period_arr" => [0.85, 1.0, 1.15],
    "forcing_amplitude" => 0.1,
    "output_figure_folder" => "arnold_tongues_traces"
)

## Arnold tongues
function params_plot_arnold_tongue(params, problem_type; n=nothing, σ=nothing)

    I_max = params["I_max"]
    forcing_period_arr = params["forcing_period_arr"]
    forcing_amplitude = params["forcing_amplitude"]
    output_figure_folder = params["output_figure_folder"]

    if problem_type == "ode"
        filename = "ode.csv"
    elseif problem_type == "sde"
        filename = "sde_n=$(n)_σ=$(σ).csv"
    end

    arnold_tongue = load_data("arnold_tongues", "model01", filename)
    entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
    println("entrainment area = $(entrainment_area_ode)")

    fig, ax = subplots(figsize=(2.5, 2.25))
    plot_arnold_tongue(arnold_tongue, I_max, ax=ax)

    ax.plot(
        forcing_period_arr,
        fill(forcing_amplitude, length(forcing_period_arr)),
        "o",
        color="red"
    )

    ax.set_title(replace(filename, "_" => ", ", ".csv" => ""))
    fig.tight_layout()

    if !isnothing(output_figure_folder)
        figure_filename = "arnold_" * replace(filename, ".csv" => ".svg")
        save_figure(fig, output_figure_folder, figure_filename)
    end

end

params_plot_arnold_tongue(params, "ode")
params_plot_arnold_tongue(params, "sde"; n=10, σ=0.005)
params_plot_arnold_tongue(params, "sde"; n=1000, σ=0.005)

## Time traces
function params_plot_time_traces(params, problem_type; n=nothing, σ=nothing)

    forcing_amplitude = params["forcing_amplitude"]
    forcing_period_arr = params["forcing_period_arr"]
    output_figure_folder = params["output_figure_folder"]

    model_name = "model01"
    i_variable = 3

    if problem_type == "ode"
        parameter_names = ["A", "τ", "I"]
        parameter_values = [0.1, 3.66, forcing_amplitude]
        trajectories = 1
        figure_filename = "time_ode.svg"
    elseif problem_type == "sde"
        parameter_names = ["A", "τ", "I", "σ"]
        parameter_values = [0.1, 3.66, forcing_amplitude, σ]
        trajectories = n
        figure_filename = "time_sde_n=$(n)_σ=$(σ).svg"
    end

    model = load_model(model_name, problem_type)
    model = set_solver(model, saveat=0.01)
    model = set_timespan(model, (0.0, 200.0))
    outfun = sol -> hcat(sol[i_variable, :])
    model = set_output(model, outfun)
    model = set_parameter_value(model, parameter_names, parameter_values)

    fig, ax_arr = subplots(length(forcing_period_arr), figsize=(7, 5))
    for (i, forcing_period) = enumerate(forcing_period_arr)
        end_time = model.problem.tspan[end]
        nevents = end_time / forcing_period + 1
        events = create_events([(:DD, forcing_period*0.5),
            (:LD, nevents, forcing_period*0.5, forcing_period*0.5)])
        model = set_input(model, events, "I")

        solution = simulate_population(model, trajectories; save_trajectories=false,
            show_progress=true)

        solution = select_time(solution, mintime=end_time-10)
        t = solution.time
        m = solution.mean[:, 1]

        ax_arr[i].plot(t, m, color="black")
        plot_events(solution.events, ax=ax_arr[i])
        ax_arr[i].set_xlabel("Time")
        ax_arr[i].set_ylabel("Z")
        if problem_type == "ode"
            ax_arr[i].set_title("ode, I = $(forcing_amplitude), \
            T = $forcing_period", loc="left")
        elseif problem_type == "sde"
            ax_arr[i].set_title("sde, n = $(n), I = $(forcing_amplitude), \
            σ = $(σ), T = $forcing_period", loc="left")
        end
        ax_arr[i].set_xlim(0, 10)
        ax_arr[i].set_ylim(minimum(m), maximum(m))
    end

    fig.tight_layout()

    if !isnothing(output_figure_folder)
        save_figure(fig, output_figure_folder, figure_filename)
    end

end

params_plot_time_traces(params, "ode")
params_plot_time_traces(params, "sde", n=10, σ=0.005)
params_plot_time_traces(params, "sde", n=1000, σ=0.005)

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
@testset "core" begin
  mp = 100
  a = fv(mp, .1, years=10, init=123)
  @test monthly_payment(a, .1, years=10, init=123) ≈ mp atol=.001
end


text:
module RocketFirstOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: first()" begin

    println("Testing: operator first()")

    run_testset([

        (
            source = from(1:42) |> first(),
            values = @ts([ 1, c ])
        ),
        (
            source = timer(50) |> first(),
            values = @ts(50 ~ [ 0, c ])
        ),
        (
            source = completed() |> first(),
            values = @ts(e(FirstNotFoundException()))
        ),
        (
            source      = completed(Int) |> first(default = "String"),
            values      = @ts([ "String", c ]),
            source_type = Union{Int, String}
        ),
        (
            source      = faulted(Int, "e") |> first(),
            values      = @ts(e("e")),
            source_type = Union{Int}
        ),
        (
            source      = faulted(Int, "e") |> first(default = "String"),
            values      = @ts(e("e")),
            source_type = Union{Int, String}
        ),
        (
            source = never() |> first(),
            values = @ts()
        )
    ])

end

end


text:
using COSMO, Test, SparseArrays, Random, LinearAlgebra

rng = Random.MersenneTwister(1872381)

P_int = 4
q_int = 2

P_Float = 4.
q_Float = 2.

P_UInt = UInt64(4)
q_UInt = UInt64(2)

P_mat = rand(rng, -10.:10., 10, 10) #float64 types
q_mat = rand(rng, -10.:10., 10, 1)

P_sparse = sprand(rng, 10, 10, 0.4) #float64 types
q_sparse = rand(rng, 1.:10.10, 10, 1)

A_vec = [1 2 3 4]
b_vec = 1

scalar_c = COSMO.Constraint(1., 2., COSMO.ZeroSet)
constr = COSMO.Constraint(Matrix(1.0I, 10, 10), rand(rng, 10), COSMO.ZeroSet)

@testset "Model" begin
     model = COSMO.Model()
     @test assemble!(model,P_int,q_int,[scalar_c]) == nothing
    #
     model = COSMO.Model{typeof(P_Float)}()
    @test assemble!(model,P_Float,q_Float,[scalar_c]) == nothing

    model = COSMO.Model()       #default type
    @test assemble!(model, P_Float, q_Float, [scalar_c]) == nothing

    # model = COSMO.Model{typeof(P_UInt)}()
    # @test assemble!(model,P_UInt,q_UInt,[scalar_c]) == nothing

    model = COSMO.Model()
    @test assemble!(model, P_mat, q_mat, [constr]) == nothing

    model = COSMO.Model()
    @test assemble!(model, P_sparse, q_sparse, [constr]) == nothing

end
nothing


text:
#=
  Badges problem

  See my analysis of the Badges problem written:
  http://www.hakank.org/data_mining/badge_problem.html
  """
  One of the best collections of machine learning data sets includes this "recreational" data set called "Badge Problem". I haven't seen any solutions to the problem anywhere,
  so here is my take.

  [This info is from the datafile:]

  1. Title: ML94/COLT94 Badge Problem

  2. Source Information
   -- Creator: Haym Hirsh, after an idea by Rob Schapire
   -- Donor: Haym Hirsh (hirsh@cs.rutgers.edu)
   -- Date: September, 1994

  3. Past Usage:
    Every pre-registered attendee at the 1994 Machine Learning
    Conference and 1994 Computational Learning Theory Conference received
    a badge labeled with a "+" or "-".  The labeling was due to some
    function known only to the badge generator (Haym Hirsh), and it
    depended only on the attendee's name.  The goal for conference
    attendees was to identify the unknown function used to generate the
    +/- labeling.

  4. Relevant Information:
    Part of the problem in using an automated program to discover the
    unknown target function is to decide how to encode names such that
    the program can be used.  The data below are presented in the form
    of a +/- label followed by the person's name.  It is up to the
    learning-system user to decide how to convert this data into something
    usable by the system (e.g., what attributes to use if your favorite
    learner requires feature-vector data).

  """

  Note: The ARFF file http://www.hakank.org/data_mining/badges_plain.arff
  includes all the preprocessed dataset, but it's more fun to do everything
  in Julia.
  I have saved the converted CSV file to badges.csv for easier handling.

  Here are the (derived) attributes:
    even_odd         {0,1} 	length of name even?
    first_char_vowel {0,1} 	is first character a vowel?
    second_char_vowel {0,1} 	is second character a vowel?
    vowels numeric 	number of vowels in the name
    consonants numeric 	number of consonants
    vowel_consonant_ratio numeric 	the ratio of vowels / consonant (this is not farfetched???)
    spaces numeric 	number of spaces
    dots: number of "." in the name, i.e. name initials
    words: number of words, i.e number of names, including name initials..

    class {1,0} 	the badge labels

  The best attribute is the params[3] (second_char_vowel) which happens to be correct.

   Here's the output of the model:

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           a    0.0452    0.0407     0.0013    0.0015   1052.1756    0.9990      227.9410
   params[1]    0.0580    0.0463     0.0015    0.0011    779.7220    0.9995      168.9172
   params[2]    0.4932    0.2944     0.0093    0.0082   1112.2659    0.9995      240.9588
   params[3]    0.7670    0.1014     0.0032    0.0046    602.2420    1.0003      130.4684
   params[4]    0.0088    0.0082     0.0003    0.0004    591.0964    1.0011      128.0538
   params[5]    0.0045    0.0040     0.0001    0.0002    863.8360    1.0053      187.1395
   params[6]    0.0808    0.0623     0.0020    0.0019    814.4369    1.0009      176.4378
   params[7]    0.4963    0.2794     0.0088    0.0049   1090.7765    0.9991      236.3034
   params[8]    0.4994    0.3009     0.0095    0.0134    891.6018    0.9999      193.1546
   params[9]    0.0207    0.0186     0.0006    0.0004    930.5087    0.9999      201.5833
       sigma    0.4964    0.2985     0.0094    0.0088    796.7233    0.9991      172.6004

...

Predictions:

Confusion matrix:
2×2 Matrix{Float64}:
 84.0    0.0
  0.0  209.0

Bad classfied:
Any[]
num bad: 0

Testing: "+ Hakan Kjellerstrand" give this (correct) prediction, i.e. that it's a "+":

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64

        y[1]    1.0067    0.9834     0.0311    0.0363   793.7538    1.0008

As we can see, the model does a pretty nice work identifying the proper class.

=#
using Turing, StatsPlots, Distributions
using LinearAlgebra
using DataFrames
using HTTP, CSV
# pyplot()
include("jl_utils.jl")

#=
# Read the dataset (from my site)
# Saved as a CSV file for easier handling in this program.
badges_url = "http://www.hakank.org/data_mining/badges.data"
dataset = CSV.read(Download.download(badges_url))
CSV.write("badges.csv",dataset)
=#

function convert_name(str)
  xcat,name = String.(split(str," ",limit=2))
  len = length(name)
  even_len = len % 2 == 0 ? 1 : 0 # Is length even?
  first_char_vowel = name[1] in vowels ? 1 : 0 # is first char in name a vowel?
  second_char_vowel = name[2] in vowels ? 1 : 0 # is second char in name a vowel?
  num_vowels = filter(c -> c in vowels, name) |> length # number of vowels
  num_consonants = filter(c -> !(c in vowels), name) |> length # number of consonants
  vowel_consonant_ratio = num_vowels / num_consonants > 0.5 ? 1 : 0  # vowels / consonants
  num_spaces = filter(c -> c == " ", name) |> length # number of spaces
  num_dots = filter(c -> c == ".", name) |> length # number of dots
  num_words = split(name, " ") |> length
          #    1        2                3                 4          5              6                     7          8        9
  return xcat, [even_len first_char_vowel second_char_vowel num_vowels num_consonants vowel_consonant_ratio num_spaces num_dots num_words]

end

df = CSV.read("badges.csv",DataFrame)
vowels = "aeiou" # note: y is not a vowel
# println(data)
x = []
y = []
for d in eachrow(df)
    cat, v = convert_name(d[1])
    push!(x,v)
    if cat == "+"
        push!(y,1)
    else
        push!(y,0)
    end
end
y = 1.0*y
x = 1.0 .* x
x = vcat(x...)
println("size(x): $(size(x))")

println("3 random instances of x:")
display(x[rand(1:size(x)[1],3),:])

@model function badges_model(x,y=missing, ::Type{T} = Float64) where {T}
    n,p = size(x)
    if y === missing
      y = Vector{T}(undef, n)
    end

    a ~ Uniform(0,1)
    params ~ filldist(Uniform(0,1),p)
    sigma ~ Uniform(0,1)

    ## Very slow, and give bad rhats
    # for i in 1:n
    #     mu = a + dot(x[i,:], params)
    #     # mu = a .+ x[i,:]' * params
    #     # y[i] ~ Normal(mu,sigma) # very slowe, bad rhat
    #     y[i] ~ Logistic(mu,sigma) # very slow, nice rhat
    # end

    # Better, much faster and neater
    mu = a .+ x*params
    y ~ MvNormal(mu, 1)

    # Return the classes for prediction.
    return y

end

model = badges_model(x,y)

# println("Prior:")
# chns = sample(model, Prior(), 10_000)
# display(chns)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 1_000)

# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)
display(chns)
# display(plot(chns))

println("\nPredictions:")

mpred = badges_model(x,missing)
pred = predict(mpred,chns)
display(pred)

cats = Dict(0=>1,1=>2)
confusion, bads = confusion_matrix(y,pred,cats,true)
println("Confusion matrix:")
display(confusion)
println("\nBad classfied:")
display(bads)
println("num bad: $(length(bads))")
if length(bads) > 0
  println("\nbads: credible_intervals")
  for b in bads
    println("b:",b)
    ix = b[1]
    println("x[$ix]: $(x[ix,:])  true cat (y[$ix]): $(y[ix])  (pred mean: ", mean(pred["y[$ix]"]))
    credible_interval(pred,"y[$(b[1])]",0.3)
  end
end

println("Testing an unknown name (I was not at the ML94/COLT94 conference)")
xcat, x = convert_name("+ Hakan Kjellerstrand")
println("xcat:$xcat x:$x")
mpred = badges_model(x,missing)
pred = predict(mpred,chns)
display(pred)
println("prediction: ", round(Int64,mean(pred["y[1]"])))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# BelgianDataset.analyze_contact_data(
#     "resources/2010_Willem_BELGIUM_participant_common.csv",
#     "resources/2010_Willem_BELGIUM_participant_extra.csv",
#     "resources/2010_Willem_BELGIUM_contact_common.csv",
#     "resources/2010_Willem_BELGIUM_contact_extra.csv",
# )

# # df_partecipant = BelgianDataset.join_partecipant_common_and_extra_datasets(
# #     "resources/2010_Willem_BELGIUM_participant_common.csv",
# #     "resources/2010_Willem_BELGIUM_participant_extra.csv",
# # )

# # occupation_distribution = BelgianDataset.get_occupation_distribution(df_partecipant, x -> x != 6 ? true : false)

# # 	filter = (age,occupation,distribution) -> occupation == 6 ?  BelgianDataset.get_occupation_single_person(age,distribution) : occupation

# # 	transform!(df_partecipant,[:part_age,:part_occupation] => ((x,y) -> filter.(x,y,occupation_distribution)) => :part_occupation)

# # df_extra = BelgianDataset.join_contact_common_and_extra_datasets("resources/2010_Willem_BELGIUM_contact_common.csv",
# # 																 "resources/2010_Willem_BELGIUM_contact_extra.csv")

# # df = CSV.File("resources/processed_contact_with_occupation.csv") |> DataFrame

# # gdf = groupby(df,:cnt_occupation)
# # length(occupation_distribution[1])

# # df_extra = BelgianDataset.join_contact_common_and_extra_datasets(
# #     "resources/2010_Willem_BELGIUM_contact_common.csv",
# #     "resources/2010_Willem_BELGIUM_contact_extra.csv",
# # )

# # df = CSV.File("resources/processed_contact_with_occupation.csv") |> DataFrame

# # gdf = groupby(df, :cnt_occupation)

# # combine(gdf, nrow)

# # ######################
# # Check-in definition

# using BelgianDataset, Dates, SimpleHypergraphs, CSV, DataFrames

# df = CSV.File("resources/processed_partecipant_and_contact_normalized.csv") |> DataFrame
# gdf = DataFrames.groupby(df, :hh_id)

# f = (x, y) -> if x == 6
#     return x + 1
# else
#     return 0
# end

# df2 = transform!(df, [:part_age, :part_occupation] => ByRow(f) => :part_occupation)

# #TODO: Check number of school nelle fiandre
# parameter = Contact_simulation_options(700, 200, 8000 , 900)

# households = BelgianDataset.start_simulation(
#     "resources/processed_partecipant_and_contact_normalized.csv",
#     parameter
# )

# # ########### Validation ################
# using JuliaDB, BelgianDataset, Dates, SimpleHypergraphs, CSV, DataFrames, Query, Statistics, DataFramesMeta

# df_generated_model, intervals, user2vertex, loc2he = BelgianDataset.generate_model_data(
#                                         "resources/generated_contact.csv",
#                                         [:Time, :Id, :position], #Column list
#                                         :Id, #userid column
#                                         :position, #venue column
#                                         :Time, #Check in column
#                                         "yyyy-mm-ddTHH:MM:SS", #Date format
#                                         Δ = Dates.Millisecond(86400000), # 24 hours
#                                         δ = Dates.Millisecond(600000), # minutes
#                                         # limit = 10000 # limit
# )
# intervals = sort(intervals)
# df_partecipant = CSV.File("resources/processed_partecipant_and_contact_normalized.csv") |> DataFrame

# df_contact = CSV.File("resources/generated_contact.csv") |> DataFrame

# #### Param config ####
# δ = Dates.Millisecond(600000)
# users = unique(df_generated_model.userid)

# lk = ReentrantLock()
# contact_per_interval = []
# for interval in intervals
#     contact_num_dict = Dict{Int,Int}()
#     current_user_count = Threads.Atomic{Int}(0)
#     t = table(filter(r -> r.timestamp > interval[2][1] && r.timestamp < interval[2][2],df_generated_model),pkey = :userid)

# Threads.@threads for user in users # For all users
#         seen = Set{Int}()
#         Threads.atomic_add!(current_user_count,1)
#         println("Current user: $current_user_count")

#         cont = 0
#         uid = user
#         ucontacts = filter(r -> r.userid == uid, t)

#         for uc in ucontacts #For each user check in
#             for r in t #For each check in of this day
#                 if uc.userid != r.userid # If I'm not checking with myself
#                     if uc.venueid == r.venueid # If we have been in the same place
#                         if abs(uc.timestamp - r.timestamp) <= δ  && r.userid ∉ seen # If the interval is small and if I haven't seen him yet
#                             # Direct contact
#                             if rand() > 0.6 #Only in 40% of the cases, they actually meet
#                                 cont = cont + 1
#                                 push!(seen,r.userid)
#                             end
#                         end
#                     end
#                 end
#             end
#         end

#         lock(lk) do
#             push!(contact_num_dict, user => cont) # Update contact num of this user
#         end
#     end

#     push!(contact_per_interval,contact_num_dict)
# end

# function normalize_age(age::Int)
#     if 18 <= age <= 29
#         return 1
#     elseif 30 <= age <= 39
#         return 2
#     elseif 40 <= age <= 49
#         return 3
#     elseif 50 <= age <= 59
#         return 4
#     elseif 60 <= age <= 69
#         return 5
#     else
#         return 6
#     end
# end

# function get_age_interval_from_normalized(age::Int)
#     if age == 1
#         return "18-29"
#     elseif age == 2
#         return "30-39"
#     elseif age == 3
#         return "40-49"
#     elseif age == 4
#         return "50-59"
#     elseif age == 5
#         return "60-69"
#     else
#         return "70+"
#     end
# end

# for i in 1:length(contact_per_interval)
#     current_dict = contact_per_interval[i]
#     sorted_contact_num_dict = sort(current_dict)

#     filtered_df_partecipant = filter(r ->  in(r.part_id,collect(keys(sorted_contact_num_dict))),df_partecipant)
#     sorted_df_partecipant = sort(filtered_df_partecipant)

#     contact_num = DataFrame(id = collect(keys(sorted_contact_num_dict)), num_contact = collect(values(sorted_contact_num_dict)), age = sorted_df_partecipant.part_age, occupation = sorted_df_partecipant.part_occupation)
#     CSV.write("contact_num_$(i)_$(parameter.work_place_num)_$(parameter.school_place_num)_$(parameter.transport_place_num)_$(parameter.leisure_place_num)_.csv",contact_num)

#     contact_num = @linq contact_num |> where(:age .> 17)
# transform!(contact_num,:age => (age -> normalize_age.(age)) => :age_normalized)
# contact_mean_per_age = DataFrame(age = [], contact_mean = [])
# for g in DataFrames.groupby(contact_num,:age_normalized)
#     push!(contact_mean_per_age,[get_age_interval_from_normalized(g.age_normalized[1]), mean(g.num_contact)])
# end

# CSV.write("contact_mean_per_age_$(i)_$(parameter.work_place_num)_$(parameter.school_place_num)_$(parameter.transport_place_num)_$(parameter.leisure_place_num).csv",sort(contact_mean_per_age))

# end

# contact_num = CSV.File("contact_num_1_700_200_8000_900_.csv") |> DataFrame

# @time distr = BelgianDataset.evaluate_direct_contacts_distribution(intervals,df,Dates.Millisecond(600000))

# df = CSV.File("resources/processed_partecipant_and_contact_normalized.csv") |> DataFrame

# subset_df = get_dataframe_subset(df,5000,[70,10,20],1)
# parameter = Contact_simulation_options(20, 3, 8000 , 100,31)

# BelgianDataset.start_contact_generation(
#     subset_df,
#     parameter
# )

# h = BelgianDataset.generatehg!(
#         nothing, # Generate a new Hypergraph
#         df, # Check in dataset
#         intervals[1][1], # mindate
#         intervals[1][2], # maxdate
#         user2vertex,
#         loc2he,
#         zeros(Int, length(unique(df.userid))), # 1 for each user: indicate if a new position for the given user must be evaluated
#         0
# )

# h

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
## Data handling

"""
    get_data(g, label)

Retrieve vertex-level metadata for the vertex with label `label`.
"""
function get_data(g::AbstractDataGraph{T,VL}, label::VL) where {T,VL}
    v = get_vertex(g, label)
    return get_data(g, v)
end

"""
    get_data(g, label_s, label_d)

Retrieve edge-level metadata for the edge between the vertices with labels `label_s` and `label_d`.
"""
function get_data(g::AbstractDataGraph{T,VL}, label_s::VL, label_d::VL) where {T,VL}
    s, d = get_vertex(g, label_s), get_vertex(g, label_d)
    return get_data(g, s, d)
end

"""
    set_data!(g, label, data)

Set vertex-level metadata to value `data` for the vertex with label `label`.
"""
function set_data!(g::AbstractDataGraph{T,VL,VD}, label::VL, data::VD) where {T,VL,VD}
    v = get_vertex(g, label)
    return set_data!(g, v, data)
end

"""
    sett_data!(g, label_s, label_d, data)

Set edge-level metadata to value `data` for the edge between the vertices with labels `label_s` and `label_d`.
"""
function set_data!(
    g::AbstractDataGraph{T,VL,VD,ED}, label_s::VL, label_d::VD, data::ED
) where {T,VL,VD,ED}
    s, d = get_vertex(g, label_s), get_vertex(g, label_d)
    return set_data!(g, s, d, data)
end


text:
using MPI, Test

include("../testhelpers.jl")

@testset "Ocean" begin
    runmpi(joinpath(@__DIR__, "HydrostaticBoussinesq/test_ocean_gyre_short.jl"))
    runmpi(joinpath(@__DIR__, "SplitExplicit/test_spindown_short.jl"))
    include(joinpath("OceanProblems", "test_initial_value_problem.jl"))
    include(joinpath("Domains", "test_rectangular_domain.jl"))
    include(joinpath("Fields", "test_rectangular_element.jl"))
    include(joinpath(
        "HydrostaticBoussinesqModel",
        "test_hydrostatic_boussinesq_model.jl",
    ))
end

max_forks_repo_name: jm-c/CLIMA


text:
# this is the benchmark "winner".
# see ../benchmarks/run-benchmarks.jl
export move_up!
function move_up!(x)
    pts_scored = 0
    updated = false

    if all(==(0), x)
        return x, pts_scored, updated
    end

    # combine phase
    @inbounds for i in 1:3
        if x[i] > 0
            for j in i+1:4
                if x[j] > 0
                    if x[i] == x[j]
                        updated = true
                        x[i] += 1
                        x[j] = 0
                        pts_scored += 1 << x[i]
                    end
                    break # break out of j loop; only combine once
                end
            end
        end
    end

    # compact phase i.e. remove 0 in between
    @inbounds for i in 1:3
        while (x[i] == 0) && any(!=(0), @view x[i+1:4])
            updated = true
            x[i:3] .= @view x[i+1:4]
            x[4] = 0
        end
    end
    x, pts_scored, updated
end