text:
import YAML
using CSV
using DataFrames
using Images: load
using StatsBase: sample
using Knet: Data

STOPWORDS = ["this", "is", "a", "an", "the",
             '.', '?', '!', ';']

function slicematrix(A::AbstractMatrix{T}) where T
    m, n = size(A)
    B = Vector{T}[Vector{T}(undef, n) for _ in 1:m]
    for i in 1:m
        B[i] .= A[i, :]
    end

    if n == 1
        return permutedims(A[:, 1])
    end
    return B
end

function readimage(imagepath; dirpath="")
    image = load(joinpath(DATADIR, dirpath, imagepath))
    image3d(image)
end

function image3d(image)
    imagergb = [float([pixel.r, pixel.g, pixel.b]) for pixel in image]
    image3d = reshape(hcat(imagergb...)', (size(image, 1), size(image, 2),
                                           length(image[1])))
end

function readtext(filename, dlm=",")
    data = split.(readlines(filename), dlm)
    features = data[1]
    rows = data[2:end]
    numfeatures = length(features)

    fdict = Dict{String, Array{T,1} where T}(
        String(features[i]) => [String.(row[i]) for row in rows]
        for i in 1:numfeatures)
end

function readyaml(filename)
    yamldict = YAML.load(open(filename))
    inputfeatures = [(dict["name"], dict["type"])
                     for dict in yamldict["input_features"]]
    outputfeatures = [(dict["name"], dict["type"])
                      for dict in yamldict["output_features"]]
    return inputfeatures, outputfeatures
end

function tokenize(text)
    punctuations = ['.', '?', '!', ';', '(', ')']
    text = replace(text, "<br />" => "")

    for p in punctuations
        text = replace(text, "$p" => " $p")
    end

    tokenized = split(text)
end

function skipwords(x; stopwords=STOPWORDS)
    removed = [[word for word in doc if !in(word, stopwords)]
               for doc in x]
end

function doc2ids(data::Array{String,1}; dict=nothing)
    dict = dict == nothing ? Dict{String, Int}() : dict
    pid = length(dict) + 1

    for x in data
        if !haskey(dict, x)
            dict[x] = pid
            pid += 1
        end
    end
    ids = [dict[x] for x in data]
    return ids, dict
end

function frequencies(x)
    freqs = Dict{String, Integer}()

    for doc in x
        for word in doc
            freqs[word] = 1 + get(freqs, word, 0)
        end
    end
    freqs
end

function vocabulary(freqs; voc=nothing, threshold=1)
    vocabulary = voc == nothing ? Dict{String, Integer}() : voc
    vocabulary["<pad>"] = 1
    vocabulary["<unk>"] = 2
    pid = length(vocabulary) + 1

    frequency_order = sort([(freqs[word], word)
                            for word in keys(freqs)], rev=true)

    for (f, w) in frequency_order
        if f >= threshold && !in(w, keys(vocabulary))
            vocabulary[w] = pid
            pid += 1
        end
    end

    vocabulary
end

maxsentence(x) = maximum(length(doc) for doc in x)

function addpadding!(x; pos="post", lenpadding=nothing)
    max_length = lenpadding != nothing ? lenpadding : maxsentence(x)

    for doc in x
        while length(doc) < max_length
            pos == "post" ? push!(doc, 1) : pushfirst!(doc, 1)
        end
    end
    x
end

function truncate!(x, th=300; pos="post")
    for doc in x
        while length(doc) > th
            pos == "post" ? pop!(doc) : popfirst!(doc)
        end
    end
    x
end

function doc2ids(x, voc)
    x = [map(lowercase, doc) for doc in x]
    ids = [map(x -> haskey(voc, x) ? voc[x] : voc["<unk>"], doc)
           for doc in x]
end

function preprocesstext(x; padding="pre", freqthreshold=10, lensentence=50,
                        voc=nothing, changevoc=false)
    x = map(doc->map(lowercase, doc), tokenize.(x))

    if voc == nothing || changevoc
        freqs = frequencies(skipwords(x))
        voc = vocabulary(freqs; threshold=freqthreshold, voc=voc)
    end

    ids = addpadding!(doc2ids(x, voc), pos=padding;
                      lenpadding=lensentence)
    ids = truncate!(ids, lensentence, pos=padding)
    ids, voc
end

function csv2data(csvpath::String)
    df = CSV.read(csvpath)
    df = cleandata(df)
    fnames = names(df)
    data = Dict{String, Array{T,1} where T}(String(fname) => Array(replace(
        df[!, fname], missing => "?")) for fname in fnames)
end

function csv2data(df::DataFrames.DataFrame)
    fnames = names(df)
    data = Dict{String, Array{T,1} where T}(String(fname) => Array(replace(
        df[!, fname], missing => "?")) for fname in fnames)
end

function cleandata(df::DataFrames.DataFrame)
    df = dropmissing(df)
end

function splitdata(df::DataFrames.DataFrame, outputcolumn=nothing; trainprop=0.8)
    categorical = outputcolumn != nothing

    if !categorical
        examplesize = size(df, 1)
        trainsize = Int(round(examplesize * trainprop))

        trnindices = sample(1:examplesize, trainsize, replace=false)
        tstindices = [i for i=1:examplesize if !in(i, trnindices)]

        trn, tst = df[trnindices, :], df[tstindices, :]
    else
        classes = unique(outputcolumn)
        divided = [df[outputcolumn .== c, :] for c in classes]
        trainframes, testframes = [], []
        for rows in divided
            examplesize = size(rows, 1)
            trainsize = Int(round(examplesize * trainprop))

            trnindices = sample(1:examplesize, trainsize, replace=false)
            tstindices = [i for i=1:examplesize if !in(i, trnindices)]

            push!(trainframes, rows[trnindices, :])
            push!(testframes, rows[tstindices, :])
        end
        trn = vcat(trainframes...)
        tst = vcat(testframes...)
    end
    trn, tst
end

function splitdata(data::Dict{String, Array{T, 1} where T},
                   outputcolumn=nothing; trainprop=0.8)
    categorical = outputcolumn != nothing

    if !categorical
        examplesize = length(iterate(values(data))[1])
        trainsize = Int(round(examplesize * trainprop))

        trnindices = sample(1:examplesize, trainsize, replace=false)
        tstindices = [i for i=1:examplesize if !in(i, trnindices)]

        trn = Dict(fname => value[trnindices] for (fname, value) in data)
        tst = Dict(fname => value[tstindices] for (fname, value) in data)
    else
        classes = unique(outputcolumn)
        examplesize = length(iterate(values(data))[1])
        classindices = [[i for i=1:examplesize if outputcolumn[i] == c]
                        for c in classes]
        train, test = [], []

        for indices in classindices
            examplesize = length(indices)
            trainsize = Int(round(examplesize * trainprop))

            trnindices = sample(indices, trainsize, replace=false)
            tstindices = [x for x in indices if !in(x, trnindices)]

            newtrn = Dict(fname => value[trnindices] for (fname, value) in data)
            newtst = Dict(fname => value[tstindices] for (fname, value) in data)

            push!(train, newtrn)
            push!(test, newtst)
        end
        trn = Dict{String, Array{T, 1} where T}(fname => vcat(
            [dict[fname] for dict in train]...) for fname in keys(data))
        tst = Dict{String, Array{T, 1} where T}(fname => vcat(
            [dict[fname] for dict in test]...) for fname in keys(data))
    end
    trn, tst
end

function splitdata(d::Data; trainprop=0.8)
    categorical = eltype(d.y) == Int

    if !categorical
        examplesize = size(d.x, 2)
        trainsize = Int(round(examplesize * trainprop))

        trnindices = sample(1:examplesize, trainsize, replace=false)
        tstindices = [i for i in 1:examplesize if !in(i, trnindices)]

        xtrn = d.x[:, trnindices]
        ytrn = d.y[:, trnindices]
        dtrn = minibatch(xtrn, ytrn, d.batchsize; shuffle=true)

        xtst = d.x[:, tstindices]
        ytst = d.y[:, tstindices]
        dtst = minibatch(xtst, ytst, d.batchsize; shuffle=true)

        return dtrn, dtst
    else
        examplesize = size(d.x, 2)
        classes = unique(d.y)
        classindices = [[i for i=1:examplesize if d.y[i] == c]
                        for c in classes]
        train, test = [], []

        for indices in classindices
            examplesize = length(indices)
            trainsize = Int(round(examplesize * trainprop))

            trnindices = sample(indices, trainsize, replace=false)
            tstindices = [i for i in indices if !in(i, trnindices)]

            xtrn = d.x[:, trnindices]
            ytrn = d.y[:, trnindices]

            xtst = d.x[:, tstindices]
            ytst = d.y[:, tstindices]

            push!(train, (xtrn, ytrn))
            push!(test, (xtst, ytst))
        end

        xtrn = hcat([x[1] for x in train]...)
        ytrn = hcat([x[2] for x in train]...)
        dtrn = minibatch(xtrn, ytrn, d.batchsize; shuffle=true)

        xtst = hcat([x[1] for x in test]...)
        ytst = hcat([x[2] for x in test]...)
        dtst = minibatch(xtst, ytst, d.batchsize; shuffle=true)

        return dtrn, dtst
    end
end

function kfolds(x, k::Int)
    n = size(x, 2)
    s = n / k
    folds = [x[:, round(Int64, (i-1)*s)+1:min(n,round(Int64, i*s))] for i=1:k]
end

function onehot(i, dims)
    v = zeros(dims)
    v[i] = 1

    return v
end

function onehotencode(data)
    dims = length(unique(data))
    ids = doc2ids(data)
    vectors = [onehot(i,  dims) for i in ids]
end

function frequentvalue(array)
    counts = [(count(x->x==i, array), i) for i in unique(array)]
    mostfrequent = maximum(counts)[2]
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
module Flux3D

using Flux, NearestNeighbors, LinearAlgebra, Makie
using Flux: @functor
using Zygote: @nograd

# utilities
include("utils.jl")

# representation
include("rep/utils.jl")
include("rep/pcloud.jl")

# visualization
include("visualize.jl")

# Dataset module
include("datasets/Dataset.jl")
using .Dataset
export ModelNet10, ModelNet40, AbstractCustomDataset, Dataset

# models
include("models/utils.jl")
include("models/dgcnn.jl")
include("models/pointnet.jl")

end # module

ext: jl


text:
assignObj(scene, superObjType, obj) = nothing

function assignObj(scene::Scene, superObjType::SuperObjCollision, obj::Object3D)
    if canCollide(obj)
      push!(superObjType.superObj, obj)
    end
end

function assignObj(scene::Scene, superObjType::SuperObjMass, obj::Object3D)
    if dataHasMass(obj)
      push!(superObjType.superObj, obj)
    end
end

function assignObj(scene::Scene, superObjType::SuperObjForce, obj::Object3D)
    if canCollide(obj) || hasJoint(obj) || hasForceElement(obj)
      push!(superObjType.superObj, obj)
    end
end

function assignObj(scene::Scene, superObjType::SuperObjVisu, obj::Object3D)
    renderer = Modia3D.renderer[1]
    if ( isVisible(obj, renderer) || (typeof(obj.visualizationFrame) != NOTHING) ) && !hasJoint(obj) && !canCollide(obj) && !hasForceElement(obj) && !hasCutJoint(obj) && !hasChildJoint(obj) #&& !dataHasMass(obj)
      push!(scene.treeVisu, obj)
    end
end

function assignAll(scene::Scene, superObj::SuperObjsRow, obj::Object3D, world::Object3D, actPos::Int64)
    names = fieldnames(typeof(superObj))
    for val in names
        tmp = getfield(superObj,val)
        assignObj(scene,tmp,obj)
    end

    if hasCutJoint(obj)
      for cutJoint in obj.twoObject3Dobject
        if typeof(cutJoint) <: Modia3D.AbstractJoint
          if cutJoint.visited
            push!(scene.noCPairsHelp[cutJoint], actPos)
    end; end; end; end

#=
    if hasForceElement(obj)
        if !(obj.forceElement in scene.forceElements)
            push!(scene.forceElements,obj.forceElement)
        end
    end
=#
  return nothing
end

function assignDynamics!(obj::Object3D)
  obj.dynamics = Object3Ddynamics()
end

function createCutJoints!(scene::Scene, obj::Object3D)
  for cutJoint in obj.twoObject3Dobject
    if typeof(cutJoint) <: Modia3D.AbstractJoint
      if !cutJoint.visited
        println("\n... Cut-joint ", ModiaMath.instanceName(cutJoint), " pushed on scene.cutJoints vector")
        push!(scene.cutJoints, cutJoint)
        push!(scene.noCPairsHelp, cutJoint => [])
        cutJoint.visited = true
  end; end; end
  return nothing
end

function fillVisuElements!(scene::Scene, obj::Object3D, world::Object3D)
  renderer            = Modia3D.renderer[1]
  visualizeFrames     = scene.options.visualizeFrames
  enableVisualization = scene.options.enableVisualization
  if enableVisualization
    if isNotCoordinateSystem(obj) && obj.visualizeFrame != Modia3D.False &&
       ((visualizeFrames && obj.visualizeFrame == Modia3D.Inherited) || obj.visualizeFrame == Modia3D.True)
      if obj != world
        obj.visualizationFrame = copyObject3D(obj, scene.autoCoordsys)
      else
        obj.visualizationFrame = copyObject3D(obj, Graphics.CoordinateSystem(2*scene.options.defaultFrameLength))
      end
      push!(scene.allVisuElements, obj.visualizationFrame)
    end
    if isVisible(obj, renderer) # isVisible(obj.data, renderer)
       push!(scene.allVisuElements, obj)
    end
  end

  scene.visualize = length(scene.allVisuElements) > 0
  return nothing
end

function assign_Visu_CutJoint_Dynamics!(scene::Scene, obj::Object3D, world::Object3D)
  assignDynamics!(obj)
  createCutJoints!(scene, obj)
  fillVisuElements!(scene, obj, world)
  return nothing
end

function assignAccVelo(tree::Vector{Object3D}, obj::Object3D)
  # compute velocity of this object
  if canCollide(obj) && !hasChildJoint(obj)
    push!(tree, obj)
  end
  # compute acceleration of this object
  if hasChildJoint(obj) || hasCutJoint(obj)
    obj.computeAcceleration = true
    push!(tree, obj)
  end
end
assignAccVelo(tree, obj) = nothing


text:
using DelimitedFiles
depths = readdlm("depths_test.txt",Int32)

@show length(depths)

prev_depth = 1000000
depth_increases = 0

for di in eachindex(depths)
    if (di+2) > length(depths)
        break
    end

    d3 = depths[di] + depths[di+1] + depths[di+2]
    if d3 > prev_depth
        global depth_increases += 1
    end
    global prev_depth = d3
end

@show depth_increases


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

using Base: @propagate_inbounds, @_inline_meta
import Base: length, size, axes, IndexStyle, getindex, setindex!, parent, vec, convert, similar

### basic definitions (types, aliases, constructors, abstractarray interface, sundry similar)

# note that Adjoint and Transpose must be able to wrap not only vectors and matrices
# but also factorizations, rotations, and other linear algebra objects, including
# user-defined such objects. so do not restrict the wrapped type.
"""
    Adjoint

Lazy wrapper type for an adjoint view of the underlying linear algebra object,
usually an `AbstractVector`/`AbstractMatrix`, but also some `Factorization`, for instance.
Usually, the `Adjoint` constructor should not be called directly, use [`adjoint`](@ref)
instead. To materialize the view use [`copy`](@ref).

This type is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims).

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> adjoint(A)
2×2 Adjoint{Complex{Int64},Array{Complex{Int64},2}}:
 3-2im  8-7im
 9-2im  4-6im
```
"""
struct Adjoint{T,S} <: AbstractMatrix{T}
    parent::S
    function Adjoint{T,S}(A::S) where {T,S}
        checkeltype_adjoint(T, eltype(A))
        new(A)
    end
end
"""
    Transpose

Lazy wrapper type for a transpose view of the underlying linear algebra object,
usually an `AbstractVector`/`AbstractMatrix`, but also some `Factorization`, for instance.
Usually, the `Transpose` constructor should not be called directly, use [`transpose`](@ref)
instead. To materialize the view use [`copy`](@ref).

This type is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims).

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> transpose(A)
2×2 Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 3+2im  8+7im
 9+2im  4+6im
```
"""
struct Transpose{T,S} <: AbstractMatrix{T}
    parent::S
    function Transpose{T,S}(A::S) where {T,S}
        checkeltype_transpose(T, eltype(A))
        new(A)
    end
end

function checkeltype_adjoint(::Type{ResultEltype}, ::Type{ParentEltype}) where {ResultEltype,ParentEltype}
    Expected = Base.promote_op(adjoint, ParentEltype)
    ResultEltype === Expected || error(string(
        "Element type mismatch. Tried to create an `Adjoint{", ResultEltype, "}` ",
        "from an object with eltype `", ParentEltype, "`, but the element type of ",
        "the adjoint of an object with eltype `", ParentEltype, "` must be ",
        "`", Expected, "`."))
    return nothing
end

function checkeltype_transpose(::Type{ResultEltype}, ::Type{ParentEltype}) where {ResultEltype, ParentEltype}
    Expected = Base.promote_op(transpose, ParentEltype)
    ResultEltype === Expected || error(string(
        "Element type mismatch. Tried to create a `Transpose{", ResultEltype, "}` ",
        "from an object with eltype `", ParentEltype, "`, but the element type of ",
        "the transpose of an object with eltype `", ParentEltype, "` must be ",
        "`", Expected, "`."))
    return nothing
end

# basic outer constructors
Adjoint(A) = Adjoint{Base.promote_op(adjoint,eltype(A)),typeof(A)}(A)
Transpose(A) = Transpose{Base.promote_op(transpose,eltype(A)),typeof(A)}(A)

Base.dataids(A::Union{Adjoint, Transpose}) = Base.dataids(A.parent)
Base.unaliascopy(A::Union{Adjoint,Transpose}) = typeof(A)(Base.unaliascopy(A.parent))

# wrapping lowercase quasi-constructors
"""
    A'
    adjoint(A)

Lazy adjoint (conjugate transposition). Note that `adjoint` is applied recursively to
elements.

For number types, `adjoint` returns the complex conjugate, and therefore it is equivalent to
the identity function for real numbers.

This operation is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims).

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> adjoint(A)
2×2 Adjoint{Complex{Int64},Array{Complex{Int64},2}}:
 3-2im  8-7im
 9-2im  4-6im

julia> x = [3, 4im]
2-element Array{Complex{Int64},1}:
 3 + 0im
 0 + 4im

julia> x'x
25 + 0im
```
"""
adjoint(A::AbstractVecOrMat) = Adjoint(A)

"""
    transpose(A)

Lazy transpose. Mutating the returned object should appropriately mutate `A`. Often,
but not always, yields `Transpose(A)`, where `Transpose` is a lazy transpose wrapper. Note
that this operation is recursive.

This operation is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims), which is non-recursive.

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> transpose(A)
2×2 Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 3+2im  8+7im
 9+2im  4+6im
```
"""
transpose(A::AbstractVecOrMat) = Transpose(A)

# unwrapping lowercase quasi-constructors
adjoint(A::Adjoint) = A.parent
transpose(A::Transpose) = A.parent
adjoint(A::Transpose{<:Real}) = A.parent
transpose(A::Adjoint{<:Real}) = A.parent

# some aliases for internal convenience use
const AdjOrTrans{T,S} = Union{Adjoint{T,S},Transpose{T,S}} where {T,S}
const AdjointAbsVec{T} = Adjoint{T,<:AbstractVector}
const TransposeAbsVec{T} = Transpose{T,<:AbstractVector}
const AdjOrTransAbsVec{T} = AdjOrTrans{T,<:AbstractVector}
const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix}

# for internal use below
wrapperop(A::Adjoint) = adjoint
wrapperop(A::Transpose) = transpose

# AbstractArray interface, basic definitions
length(A::AdjOrTrans) = length(A.parent)
size(v::AdjOrTransAbsVec) = (1, length(v.parent))
size(A::AdjOrTransAbsMat) = reverse(size(A.parent))
axes(v::AdjOrTransAbsVec) = (Base.OneTo(1), axes(v.parent)...)
axes(A::AdjOrTransAbsMat) = reverse(axes(A.parent))
IndexStyle(::Type{<:AdjOrTransAbsVec}) = IndexLinear()
IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian()
@propagate_inbounds getindex(v::AdjOrTransAbsVec, i::Int) = wrapperop(v)(v.parent[i-1+first(axes(v.parent)[1])])
@propagate_inbounds getindex(A::AdjOrTransAbsMat, i::Int, j::Int) = wrapperop(A)(A.parent[j, i])
@propagate_inbounds setindex!(v::AdjOrTransAbsVec, x, i::Int) = (setindex!(v.parent, wrapperop(v)(x), i-1+first(axes(v.parent)[1])); v)
@propagate_inbounds setindex!(A::AdjOrTransAbsMat, x, i::Int, j::Int) = (setindex!(A.parent, wrapperop(A)(x), j, i); A)
# AbstractArray interface, additional definitions to retain wrapper over vectors where appropriate
@propagate_inbounds getindex(v::AdjOrTransAbsVec, ::Colon, is::AbstractArray{Int}) = wrapperop(v)(v.parent[is])
@propagate_inbounds getindex(v::AdjOrTransAbsVec, ::Colon, ::Colon) = wrapperop(v)(v.parent[:])

# conversion of underlying storage
convert(::Type{Adjoint{T,S}}, A::Adjoint) where {T,S} = Adjoint{T,S}(convert(S, A.parent))
convert(::Type{Transpose{T,S}}, A::Transpose) where {T,S} = Transpose{T,S}(convert(S, A.parent))

# Strides and pointer for transposed strided arrays — but only if the elements are actually stored in memory
Base.strides(A::Adjoint{<:Real, <:AbstractVector}) = (stride(A.parent, 2), stride(A.parent, 1))
Base.strides(A::Transpose{<:Any, <:AbstractVector}) = (stride(A.parent, 2), stride(A.parent, 1))
# For matrices it's slightly faster to use reverse and avoid calling stride twice
Base.strides(A::Adjoint{<:Real, <:AbstractMatrix}) = reverse(strides(A.parent))
Base.strides(A::Transpose{<:Any, <:AbstractMatrix}) = reverse(strides(A.parent))

Base.unsafe_convert(::Type{Ptr{T}}, A::Adjoint{<:Real, <:AbstractVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)
Base.unsafe_convert(::Type{Ptr{T}}, A::Transpose{<:Any, <:AbstractVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)

Base.elsize(::Type{<:Adjoint{<:Real, P}}) where {P<:AbstractVecOrMat} = Base.elsize(P)
Base.elsize(::Type{<:Transpose{<:Any, P}}) where {P<:AbstractVecOrMat} = Base.elsize(P)

# for vectors, the semantics of the wrapped and unwrapped types differ
# so attempt to maintain both the parent and wrapper type insofar as possible
similar(A::AdjOrTransAbsVec) = wrapperop(A)(similar(A.parent))
similar(A::AdjOrTransAbsVec, ::Type{T}) where {T} = wrapperop(A)(similar(A.parent, Base.promote_op(wrapperop(A), T)))
# for matrices, the semantics of the wrapped and unwrapped types are generally the same
# and as you are allocating with similar anyway, you might as well get something unwrapped
similar(A::AdjOrTrans) = similar(A.parent, eltype(A), axes(A))
similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, axes(A))
similar(A::AdjOrTrans, ::Type{T}, dims::Dims{N}) where {T,N} = similar(A.parent, T, dims)

# sundry basic definitions
parent(A::AdjOrTrans) = A.parent
vec(v::TransposeAbsVec) = parent(v)
vec(v::AdjointAbsVec{<:Real}) = parent(v)

### concatenation
# preserve Adjoint/Transpose wrapper around vectors
# to retain the associated semantics post-concatenation
hcat(avs::Union{Number,AdjointAbsVec}...) = _adjoint_hcat(avs...)
hcat(tvs::Union{Number,TransposeAbsVec}...) = _transpose_hcat(tvs...)
_adjoint_hcat(avs::Union{Number,AdjointAbsVec}...) = adjoint(vcat(map(adjoint, avs)...))
_transpose_hcat(tvs::Union{Number,TransposeAbsVec}...) = transpose(vcat(map(transpose, tvs)...))
typed_hcat(::Type{T}, avs::Union{Number,AdjointAbsVec}...) where {T} = adjoint(typed_vcat(T, map(adjoint, avs)...))
typed_hcat(::Type{T}, tvs::Union{Number,TransposeAbsVec}...) where {T} = transpose(typed_vcat(T, map(transpose, tvs)...))
# otherwise-redundant definitions necessary to prevent hitting the concat methods in sparse/sparsevector.jl
hcat(avs::Adjoint{<:Any,<:Vector}...) = _adjoint_hcat(avs...)
hcat(tvs::Transpose{<:Any,<:Vector}...) = _transpose_hcat(tvs...)
hcat(avs::Adjoint{T,Vector{T}}...) where {T} = _adjoint_hcat(avs...)
hcat(tvs::Transpose{T,Vector{T}}...) where {T} = _transpose_hcat(tvs...)
# TODO unify and allow mixed combinations

### higher order functions
# preserve Adjoint/Transpose wrapper around vectors
# to retain the associated semantics post-map/broadcast
#
# note that the caller's operation f operates in the domain of the wrapped vectors' entries.
# hence the adjoint->f->adjoint shenanigans applied to the parent vectors' entries.
map(f, avs::AdjointAbsVec...) = adjoint(map((xs...) -> adjoint(f(adjoint.(xs)...)), parent.(avs)...))
map(f, tvs::TransposeAbsVec...) = transpose(map((xs...) -> transpose(f(transpose.(xs)...)), parent.(tvs)...))
quasiparentt(x) = parent(x); quasiparentt(x::Number) = x # to handle numbers in the defs below
quasiparenta(x) = parent(x); quasiparenta(x::Number) = conj(x) # to handle numbers in the defs below
broadcast(f, avs::Union{Number,AdjointAbsVec}...) = adjoint(broadcast((xs...) -> adjoint(f(adjoint.(xs)...)), quasiparenta.(avs)...))
broadcast(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs...) -> transpose(f(transpose.(xs)...)), quasiparentt.(tvs)...))
# Hack to preserve behavior after #32122; this needs to be done with a broadcast style instead to support dotted fusion
Broadcast.broadcast_preserving_zero_d(f, avs::Union{Number,AdjointAbsVec}...) = adjoint(broadcast((xs...) -> adjoint(f(adjoint.(xs)...)), quasiparenta.(avs)...))
Broadcast.broadcast_preserving_zero_d(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs...) -> transpose(f(transpose.(xs)...)), quasiparentt.(tvs)...))
# TODO unify and allow mixed combinations with a broadcast style

### linear algebra

(-)(A::Adjoint)   = Adjoint(  -A.parent)
(-)(A::Transpose) = Transpose(-A.parent)

## multiplication *

# Adjoint/Transpose-vector * vector
*(u::AdjointAbsVec{T}, v::AbstractVector{T}) where {T<:Number} = dot(u.parent, v)
*(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = dot(u.parent, v)
*(u::AdjOrTransAbsVec, v::AbstractVector) = sum(uu*vv for (uu, vv) in zip(u, v))

# vector * Adjoint/Transpose-vector
*(u::AbstractVector, v::AdjOrTransAbsVec) = broadcast(*, u, v)
# Adjoint/Transpose-vector * Adjoint/Transpose-vector
# (necessary for disambiguation with fallback methods in linalg/matmul)
*(u::AdjointAbsVec, v::AdjointAbsVec) = throw(MethodError(*, (u, v)))
*(u::TransposeAbsVec, v::TransposeAbsVec) = throw(MethodError(*, (u, v)))

# AdjOrTransAbsVec{<:Any,<:AdjOrTransAbsVec} is a lazy conj vectors
# We need to expand the combinations to avoid ambiguities
(*)(u::TransposeAbsVec, v::AdjointAbsVec{<:Any,<:TransposeAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))
(*)(u::AdjointAbsVec,   v::AdjointAbsVec{<:Any,<:TransposeAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))
(*)(u::TransposeAbsVec, v::TransposeAbsVec{<:Any,<:AdjointAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))
(*)(u::AdjointAbsVec,   v::TransposeAbsVec{<:Any,<:AdjointAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))

## pseudoinversion
pinv(v::AdjointAbsVec, tol::Real = 0) = pinv(v.parent, tol).parent
pinv(v::TransposeAbsVec, tol::Real = 0) = pinv(conj(v.parent)).parent

## left-division \
\(u::AdjOrTransAbsVec, v::AdjOrTransAbsVec) = pinv(u) * v

## right-division /
/(u::AdjointAbsVec, A::AbstractMatrix) = adjoint(adjoint(A) \ u.parent)
/(u::TransposeAbsVec, A::AbstractMatrix) = transpose(transpose(A) \ u.parent)
/(u::AdjointAbsVec, A::Transpose{<:Any,<:AbstractMatrix}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent)
/(u::TransposeAbsVec, A::Adjoint{<:Any,<:AbstractMatrix}) = transpose(conj(A.parent) \ u.parent) # technically should be transpose(copy(transpose(copy(A))) \ u.parent)

## complex conjugate
conj(A::Transpose) = adjoint(A.parent)
conj(A::Adjoint) = transpose(A.parent)

max_stars_repo_name: martincornejo/julia

max_issues_repo_name: martincornejo/julia

max_forks_repo_name: martincornejo/julia


text:
"""
    idx_for(i::Int)

Simple helper function for converting 0-based indices to Julia's 1-based indices.
"""
function idx_for(i::Int)
    return i + 1
end

"""
    heap_sort!(arr::Vector{T}, gt = >, N::Int = length(arr)) where {T}

Sort the given vector (in-place) using the Heapsort algorithm.

Heapsort consists of two stages:
  1. Building a (max) heap of the array
  2. Repeatedly extracting the largest element and inserting it at the front of the sorted part of the array

After the largest element has been extracted, the tree is updated to maintain the heap property via a "sifting" operation.

  Storing a heap in an array is pretty straightforward - for every node with index n, its children are stored at indices
  2n + 1 and 2n + 2 (for 0-based indices). Index 0 contains the root node.
  Since Julia's indices are 1-based, we need to change this a little bit. We're using a trivial helper function
  idx_for to convert from 0-based to 1-based.

  See https://en.wikipedia.org/wiki/Heapsort for a complete explanation of Heapsort.

  Contributed By:- [Frank Schmitt](https://github.com/frankschmitt)
"""
function heap_sort!(arr::Vector{T}, gt = >, N::Int = length(arr)) where {T}
    n = N
    i = div(n, 2)
    t = -1

    while true
        if i > 0    # First stage - Sorting the heap
            i -= 1    # Save its index to i
            t = arr[idx_for(i)]    # Save parent value to t
        else          # Second stage - Extracting elements in-place
            n -= 1    # Make the new heap smaller
            if (n == 0)
                return # When the heap is empty, we are done
            end
            t = arr[idx_for(n)]               # Save last value (it will be overwritten)
            arr[idx_for(n)] = arr[idx_for(0)] # Save largest value at the end of arr
        end

        parent = i # We will start pushing down t from parent
        child = i * 2 + 1 # parent's left child

        # Sift operation - pushing the value of t down the heap
        while child < n
            if (
                child + 1 < n &&
                gt(arr[idx_for(child + 1)], arr[idx_for(child)])
            )
                child += 1 # Choose the largest child
            end
            if (gt(arr[idx_for(child)], t))                   # If any child is bigger than the parent
                arr[idx_for(parent)] = arr[idx_for(child)]    # Move the largest child up
                parent = child                                # Move parent pointer to this child
                child = parent * 2 + 1                            #  find the next child
            else
                break # t's place is found
            end
        end
        arr[idx_for(parent)] = t
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
    AbsoluteValue()

Ensures that the first variable is the absolute value of the second one.
"""
struct AbsoluteValue <: MOI.AbstractVectorSet
end

MOI.dimension(set::AbsoluteValue) = 2
copy(set::AbsoluteValue) = set
Base.:(==)(::AbsoluteValue, ::AbsoluteValue) = true

"""
    Modulo()

Ensures that the three variables are related as:

``x \\mod y \\equiv z``
"""
struct Modulo <: MOI.AbstractVectorSet
end

MOI.dimension(set::Modulo) = 3
copy(set::Modulo) = set
Base.:(==)(::Modulo, ::Modulo) = true


text:
## Helper functions
"""
    randID(k::Int)

return random number between 1 and k
"""
function randID(k::Int)
    return Int(ceil(k * rand()))
end

"""
    likelihoodratio(m::Array{Float32}, mtest::Array{Float32}, d::Array{Float32})

Generic Likelihood calculation using Poisson noise statistics
"""
function likelihoodratio(m::Array{Float32}, mtest::Array{Float32}, d::Array{Float32})
    nelem=prod(size(m))
    minmodel=1f-6 #to avoid division by zero and log of zero
    LLR = 0;
    for ii = 1:nelem
        LLR += m[ii] - mtest[ii] + d[ii] * log(max(mtest[ii],minmodel) / max(m[ii],minmodel));
    end
    L = exp(LLR)
    if L < 0
        println(L, LLR)
    end
    return exp(LLR)
end

""""
    poissrnd!(d::Array{Float32})
In-place Poisson noise corruptor
"""
function poissrnd!(d::Array{Float32})
    nelem=prod(size(d))
    for nn = 1:nelem
        d[nn] = Float32(rand(Poisson(Float64(d[nn]))))
    end
end

# "An inexpensive Normal random number generator for CUDA"
function curandn()
    r=0;
    for ii=1:12
        r+=rand()
    end
    return r-6f0
end

"""
    calcintialstate(rjs::RJStruct, seed::Int32)

calculates an intial state by finding the maximum of the data and places an emitter
at that location with an intensity selected from the prior. Calculates the model for the
initial state and stores in rjs.modeldata

allows the input of an optional non-negative, integer random seed that generates from
the Mersenne Twister library to produce deterministic initial states. If no integer is input,
the function uses the OS provided entropy.

"""
function calcintialstate(rjs::RJStruct, seed::Int32=Int32(-1)) # find initial state for direct detection data
    if seed == -1
    else
        Random.seed!(seed);
    end
    d = rjs.data.data
    state1 = StateFlatBg()
    state1.n = 1
    state1.bg = minimum(d)
    roimax = maximum(d)
    coords = findall(x -> x == roimax, d)
    state1.y[1] = coords[1].I[1]
    state1.x[1] = coords[1].I[2]
    state1.photons[1] = rand(rjs.prior_photons)

    genmodel!(state1,rjs.psf,rjs.modeldata)
    return state1
end


text:
# To be used instead of _PMACDC.variable_dc_converter_ne() - supports deduplication of variables
function variable_dc_converter_ne(pm::_PM.AbstractDCPModel; kwargs...)
    variable_ne_converter_indicator(pm; kwargs..., relax=true) # FlexPlan version: replaces _PMACDC.variable_converter_ne().
    variable_ne_converter_investment(pm; kwargs...)
    _PMACDC.variable_converter_active_power_ne(pm; kwargs...)
    _PMACDC.variable_dcside_power_ne(pm; kwargs...)
    _PMACDC.variable_converter_filter_voltage_ne(pm; kwargs...)
    _PMACDC.variable_converter_internal_voltage_ne(pm; kwargs...)
    _PMACDC.variable_converter_to_grid_active_power_ne(pm; kwargs...)

    _PMACDC.variable_conv_transformer_active_power_to_ne(pm; kwargs...)
    _PMACDC.variable_conv_reactor_active_power_from_ne(pm; kwargs...)
end


text:
# Tudor Berariu <tudor.berariu@gmail.com>, 2015

module Bullseye

export bullseye

function bullseye{T<:FloatingPoint}(N::Integer, D::Integer,
                                    radii::Array{T,1}, widths::Array{T,1})
    if length(radii) != length(widths)
        error("Parameters radii and widths must have same length");
    end
    const M = length(radii);
    data = rand(T, D, N);
    labels = rand(1:M, N);
    for i = 1:N
        len = radii[labels[i]] + randn() * widths[labels[i]] * 0.5;
        angle = rand(D-1)*2*pi;
        sins = one(T);
        for j = 1:D-1
            data[j,i] = len * sins * cos(angle[j]);
            sins *= sin(angle[j]);
        end
        data[D,i] = len * sins;
    end
    return data, labels;
end

function bullseye2{T<:FloatingPoint}(N::Integer, D::Integer,
                                    radii::Array{T,1}, widths::Array{T,1})
    if length(radii) != length(widths)
        error("Parameters radii and widths must have same length");
    end
    const M = length(radii);
    data = rand(T, D, N);
    labels = rand(1:M, N);
    len = radii[labels] + randn(N) .* widths[labels] * 0.5;
    angle = rand(N, D-1) * 2 * pi;
    for i = 1:N
        len = radii[labels[i]] + randn() * widths[labels[i]] * 0.5;
        angle = rand(D-1)*2*pi;
        sins = one(T);
        for j = 1:D-1
            data[j,i] = len * sins * cos(angle[j]);
            sins *= sin(angle[j]);
        end
        data[D,i] = len * sins;
    end
    return data, labels;
end

end

max_stars_repo_name: tudor-berariu/ml-hacks

max_issues_repo_name: tudor-berariu/ml-hacks

max_forks_repo_name: tudor-berariu/ml-hacks


text:
struct StaticTraitLayer end

struct Trait{V} <: Measurement
    val::V
end

Base.rand(::Type{Trait}, x::Int, y::Int) = Base.rand(Trait, 0.5, x, y)
Base.rand(::Type{Trait}, dist::Distribution, x::Int, y::Int) = begin
    z = zeros(Trait, x, y)
    for i in eachindex(z)
        z[i] = rand(dist)
    end
    return z
end

Base.zero(::Trait{T}) where {T} = Trait{T}(zero(T))
Base.zero(::Type{Trait}) = Trait{Float64}(0)
Base.zero(::Type{<:Trait{T}}) where {T} = Trait(zero(T))

Base.fill(Trait, v, dims) = Base.fill(Trait(v), dims...)


text:
# Good old CICR with

using Parameters
include("common.jl")

@with_kw struct LCCParams
    A = 2.0
    B = 2.0
    ω = 0.01
    f = 0.3
    g = 2.0
    P_CA = 1.24E-3
    P_K = 1.11E-11
    I_CA_HALF = -0.4583
end

@with_kw struct RYRParams
    R_RYR = 3.6
    N = 4
    M = 3
    KA_P = 1.125E10
    KA_M = 0.576
    KB_P = 4.05E6
    KB_M = 1.93
    KC_P = 0.1
    KC_M = 8E-4
end

# LCC system
function lcc_system(c0, c1, c2, c3, c4, o, cca0, cca1, cca2, cca3, ca_ss, vm, pLCC::LCCParams)
    @unpack A, B, ω, f, g = pLCC
    cca4 = one(c0) - c0 - c1 - c2 - c3 - c4 - o - cca0 - cca1 - cca2 - cca3
    α = 0.4 * exp((vm + 2) / 10)
    β = 0.05 * exp(-(vm+2) / 13)
    α′ = A * α
    β′ = β / B
    γ = 0.1875 * ca_ss
    v01 = 4 * α * c0 - β * c1
    v12 = 3 * α * c1 - 2 * β * c2
    v23 = 2 * α * c2 - 3 * β * c3
    v34 = α * c3 - 4 * β * c4
    v45 = f * c4 - g * o
    v67 = 4 * α′ * cca0 - β′ * cca1
    v78 = 3 * α′ * cca1 - 2 * β′ * cca2
    v89 = 2 * α′ * cca2 - 3 * β′ * cca3
    v910 = α′ * cca3 - 4 * β′ * cca4
    v06 = γ * c0 - ω * cca0
    v17 = A * γ * c1 - ω / B * cca1
    v28 = (A^2) * γ * c2 - ω / (B^2)  * cca2
    v39 = (A^3) * γ * c3 - ω  / (B^3) * cca3
    v410 = (A^4) * γ * c4 - ω / (B^4) * cca4
    d_c0 = -v01 - v06
    d_c1 = v01 - v12 - v17
    d_c2 = v12 - v23 - v28
    d_c3 = v23 - v34 - v39
    d_c4 = v34 - v45 - v410
    d_o = v45
    d_cca0 = v06 - v67
    d_cca1 = v17 + v67 - v78
    d_cca2 = v28 + v78 - v89
    d_cca3 = v39 + v89 - v910
    # d_cca4 = v410 + v910
    return (d_c0, d_c1, d_c2, d_c3, d_c4, d_o, d_cca0, d_cca1, d_cca2, d_cca3)
end

# Voltage inactivated factor of LCC
function dyca(y_ca, vm)
    y_inf = inv(1 + exp((vm+55) / 7.5)) + 0.5 / (1 + exp((-vm + 21) / 6))
    τ = 20.0 + 600.0 / (1 + exp((vm + 30) / 9.5))
    return (y_inf - y_ca) / τ
end

function ical(y_ca, o, vfrt, evfrtm1, k_i, k_o, ca_o, pLCC::LCCParams)
    @unpack P_CA, P_K, I_CA_HALF = pLCC
    ev2frtm1 = evfrtm1 * (evfrtm1 + 2)
    iCaMax = _ghk(2 * vfrt, ev2frtm1, P_CA, 0.001, 0.341 * ca_o, 2)
    iCaL = 6 * y_ca * o * iCaMax
    iCaK = _mm(I_CA_HALF, iCaMax) * y_ca * o * _ghk(vfrt, evfrtm1, P_K, k_i, k_o)
    return (iCaL, iCaK)
end

function ryr_system(po1, po2, pc2, ca_ss, pRYR::RYRParams, CA_SS_THR=0.095)
    @unpack N, M, KA_P, KA_M, KB_P, KB_M, KC_P, KC_M = pRYR
    pc1 = one(po1) - (po1 + po2 + pc2)
    # Switch formulation (rapid equlibrium assumption) by Plank et al. (2008)
    KA_P_CA = KA_P * ca_ss^N
    if ca_ss >= CA_SS_THR
        po1 = (po1 + pc1) * _mm(KA_P_CA, KA_M)
        d_po1 = zero(po1)
    else
        d_po1 = -KA_M * po1 + KA_P_CA * pc1
    end
    vo1o2 = KB_P * ca_ss^M * po1 - KB_M * po2
    vo1c2 = KC_P * po1 - KC_M * pc2
    d_po1 += -vo1o2 - vo1c2
    d_po2 = vo1o2
    d_pc2 = vo1c2
    return (d_po1, d_po2, d_pc2)
end

jrel(po1, po2, ca_jsr, ca_ss, R_RYR) = R_RYR * (po1 + po2) * (ca_jsr - ca_ss)


text:
#####
##### Naive Implementation
#####

function reuse_naive(v::AbstractVector{T}) where {T}
    lastuse = Dict{T,Int}()
    counts = Set{T}()
    histogram = Dict{Int,Int}()

    ProgressMeter.@showprogress 1 for (t, i) in enumerate(v)
        lasttime = get(lastuse, i, t)
        d = -1
        if lasttime != t
            empty!(counts)
            for u in view(v, (lasttime + 1):(t - 1))
                push!(counts, u)
            end
            d = length(counts)
        end
        histogram[d] = 1 + get(histogram, d, 0)
        lastuse[i] = t
    end
    return histogram
end

function gramsequal(a, b)
    passed = true
    for (k, v) in a
        if haskey(b, k)
            passed &= (b[k] == v)
        else
            passed = false
        end
    end

    for (k, v) in b
        if haskey(a, k)
            passed &= (a[k] == v)
        else
            passed = false
        end
    end
    return passed
end


text:
@testset "Saddle-Point Systems" begin

  @testset "Matrix tests" begin

    A1 = Float64[1 2; 2 1]
    B2 = Float64[2 3;-1 -1]
    B1 = B2'
    C = Matrix{Float64}(undef,2,2)
    C .= [5 -2; 3 -4];

    A = SaddleSystem(A1,B2,B1,C)

    # test inputs and outputs as tuples
    rhs = ([1.0,2.0],[3.0,4.0])
    sol = (zeros(2),zeros(2))

    ldiv!(sol,A,rhs)

    @test norm((A*sol)[1]-rhs[1]) < 1e-14
    @test norm((A*sol)[2]-rhs[2]) < 1e-14

    sol2 = A\rhs

    @test norm((A*sol2)[1]-rhs[1]) < 1e-14
    @test norm((A*sol2)[2]-rhs[2]) < 1e-14

    # test inputs and outputs as vectors
    rhs1, rhs2 = rhs
    rhsvec = [rhs1;rhs2]
    solvec = similar(rhsvec)

    solvec = A\rhsvec

    @test norm(A*solvec-rhsvec) < 1e-14

    # test with the other constructor
    Aother = SaddleSystem(A1,B2,B1,C,zeros(Float64,2),zeros(Float64,2))

    sol_other = Aother\rhs

    @test norm((Aother*sol_other)[1]-rhs[1]) < 1e-14
    @test norm((Aother*sol_other)[2]-rhs[2]) < 1e-14

    # Test with SaddleVector
    rhs = SaddleVector(rhs1,rhs2)
    As = SaddleSystem(A1,B2,B1,C,rhs)

    sol = As\rhs

    @test norm(state(As*sol)-state(rhs)) < 1e-14
    @test norm(constraint(As*sol)-constraint(rhs)) < 1e-14

  end

  nx = 130; ny = 130
  Lx = 2.0
  dx = Lx/(nx-2)
  w = Nodes(Dual,(nx,ny))
  q = Edges(Primal,w)

  L = plan_laplacian(size(w),with_inverse=true)

  n = 128
  θ = range(0,stop=2π,length=n+1)
  R = 0.5
  xb = 1.0 .+ R*cos.(θ[1:n])
  yb = 1.0 .+ R*sin.(θ[1:n])
  ds = (2π/n)*R
  X = VectorData(xb,yb)
  f = ScalarData(X)

  E = Regularize(X,dx;issymmetric=true)
  Hmat,Emat = RegularizationMatrix(E,f,w)

  @testset "Construction of linear maps" begin

    u = similar(w)
    wvec = vec(w)

    w .= rand(size(w)...)

    uvec = zeros(length(w))

    Lop = ConstrainedSystems.linear_map(L,w)

    u = L*w

    uvec = Lop*wvec

    @test ConstrainedSystems._wrap_vec(uvec,u) == u

    Linv = ConstrainedSystems.linear_inverse_map(L,w)

    yvec = Linv*wvec

    y = L\w

    @test ConstrainedSystems._wrap_vec(yvec,y) == y

    # point-wise operators
    fvec = vec(f)

    f[10] = 1.0
    Hop = ConstrainedSystems.linear_map(Hmat,f,w);

    y = Hmat*f

    yvec = Hop*fvec

    @test ConstrainedSystems._wrap_vec(yvec,y) == y

    Eop = ConstrainedSystems.linear_map(Emat,w,f)

    g = Emat*w

    gvec = Eop*vec(w)

    @test ConstrainedSystems._wrap_vec(gvec,g) == g

    # Test in-place function wrapping
    my_ip_fcn!(outp::Vector,inp::Vector) = outp .= 2.0*inp
    inp = rand(1000)
    outp = zero(inp)
    my_ip_fcn!(outp,inp)

    my_ip_lm! = ConstrainedSystems._create_vec_function!(my_ip_fcn!,outp,inp)
    outp2 = zero(inp)
    my_ip_lm!(outp2,inp)
    @test outp2 == outp

    curl_lm! = ConstrainedSystems._create_vec_function!(CartesianGrids.curl!,q,w)
    CartesianGrids.curl!(q,w)

    uvec .= wvec
    qvec = zeros(length(q))
    curl_lm!(qvec,uvec)

    @test ConstrainedSystems._wrap_vec(qvec,q) == q

    curl_lm! = ConstrainedSystems.linear_map(CartesianGrids.curl!,w,q)
    mul!(qvec,curl_lm!,uvec)

    @test ConstrainedSystems._wrap_vec(qvec,q) == q

  end

  ψb = ScalarData(X)
  w = Nodes(Dual,(nx,ny))
  ψb .= -(xb .- 1)
  f .= ones(Float64,n)*ds
  ψ = Nodes(Dual,w)

  nada = empty(f)

  @testset "Field operators" begin

    rhs = SaddleVector(w,ψb)

    A = SaddleSystem(L,Emat,Hmat,rhs)

    A.S*vec(f)

    A.S⁻¹*vec(f)

    sol = deepcopy(rhs)

    ldiv!(sol,A,rhs)

    sol2 = A\rhs

    @test state(sol2) == state(sol)
    @test constraint(sol2) == constraint(sol)

    rhs2 = A*sol2

    @test norm(constraint(rhs2)-ψb) < 1e-14

    fex = -2*cos.(θ[1:n])
    @test norm(constraint(sol)-fex*ds) < 0.02
    @test ψ[nx,65] ≈ -ψ[1,65]
    @test ψ[65,ny] ≈ ψ[65,1]

  end

  fv = VectorData(X)
  q = Edges(Primal,w)
  Hvmat,Evmat = RegularizationMatrix(E,fv,q)

  @testset "Vector force data" begin

    B₁ᵀ(f) = curl(Hvmat*f)
    B₂(w) = -(Evmat*(curl(L\w)))

    rhsv = SaddleVector(w,fv)
    A = SaddleSystem(I,B₂,B₁ᵀ,rhsv)

    sol = A\rhsv
    @test state(sol) == zero(w) && constraint(sol) == zero(fv)

  end

  Ẽ = Regularize(X,dx;weights=ds,filter=true)
  H̃mat = RegularizationMatrix(Ẽ,ψb,w)
  Ẽmat = InterpolationMatrix(Ẽ,w,ψb)
  Pmat = Ẽmat*H̃mat

  @testset "Filtering" begin

    rhs = SaddleVector(w,ψb)
    Afilt = SaddleSystem(L,Emat,Hmat,rhs,filter=Pmat)

    sol = Afilt\rhs

    fex = -2*cos.(θ[1:n])
    @test norm(constraint(sol)-fex*ds) < 0.01

  end

  @testset "Reduction to unconstrained system" begin

    op = linear_map(nothing,nada)

    @test size(op) == (0,0)
    @test op*nada == ()

    op = linear_map(nothing,w,nada)

    @test size(op) == (0,length(w))
    @test op*vec(w) == ()

    op = linear_map(nothing,nada,w)

    @test size(op) == (length(w),0)
    @test op*nada == vec(zero(w))

    rhsnc = SaddleVector(w,nada)

    Anc = SaddleSystem(L,rhsnc)

    sol = Anc\rhsnc
    ψ = state(sol)
    fnull = constraint(sol)

    @test ψ == L\w
    @test fnull == nada

    q.u .= 3;
    rhsqnc = SaddleVector(q,nada)
    Aqnc = SaddleSystem(I,rhsqnc)

    sol = Aqnc\rhsqnc

    q2 = state(sol)
    fnull = constraint(sol)

    @test all(q2.u .== 3.0) && all(q2.v .== 0.0)

  end

  @testset "Tuple of saddle point systems" begin

    rhs = SaddleVector(w,ψb)
    A = SaddleSystem(L,Emat,Hmat,rhs)
    sol = A\rhs
    ψ = state(sol)
    f = constraint(sol)

    rhsnc = SaddleVector(w,nada)
    Anc = SaddleSystem(L,rhsnc)

    sys = (A,Anc)
    sol1, sol2 = sys\(rhs,rhsnc)

    @test state(sol1) == ψ && constraint(sol1) == f && state(sol2) == zero(w) && constraint(sol2) == nada

    newrhs1,newrhs2 = sys*(sol1,sol2)

    @test newrhs1 == A*sol1 && newrhs2 == Anc*sol2

  end

  @testset "Recursive saddle point with vectors" begin
    A1 = Float64[1 2; 2 1]
    B21 = Float64[2 3]
    B11 = B21'
    C1 = Matrix{Float64}(undef,1,1)
    C1.= 5

    B22 = Float64[-1 -1 3]
    B12 = Float64[-1 -1 -2]'
    C2 = Matrix{Float64}(undef,1,1)
    C2.= -4

    rhs11 = [1.0,2.0]
    rhs12 = Vector{Float64}(undef,1)
    rhs12 .= 3.0
    #rhs1 = (rhs11,rhs12)
    rhs1 = [rhs11;rhs12]

    rhs2 = Vector{Float64}(undef,1)
    rhs2 .= 4.0

    rhs = (rhs1,rhs2)

    A = SaddleSystem(A1,B21,B11,C1)

    Abig = SaddleSystem(A,B22,B12,C2,rhs1,rhs2)

    sol = Abig\rhs

    out = Abig*sol

    @test norm(out[1]-rhs1) < 1e-14

  end

end


text:
### Fetch required packages and reaction networks ###
using DiffEqBase, Catalyst, Random, Test

using StableRNGs
rng = StableRNG(12345)

### Tests various cutom made functions ###
new_hill(x, v, k, n) = v*x^n/(k^n+x^n)
new_poly(x,p1,p2) = .5*p1*x^2
new_exp(x,p) = exp(-p*x)

custom_function_network_1 = @reaction_network begin
    hill(X,v1,K1,2), X + Y --> Z1
    mm(X,v2,K2), X + Y --> Z2
    p1*X^2+p2, X + Y --> Z3
    exp(-p3*Y), X + Y --> Z4
    hillR(X,v3,K3,2), X + Y --> Z5
    mmR(X,v4,K4), X + Y --> Z6
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4

custom_function_network_2 = @reaction_network begin
    new_hill(X,v1,K1,2), X + Y --> Z1
    v2*X/(X+K2), X + Y --> Z2
    2*new_poly(X,p1,p2)+p2, X + Y --> Z3
    new_exp(Y,p3), X + Y --> Z4
    v3*(K3^2)/(K3^2+X^2), X + Y --> Z5
    v4*K4/(X+K4), X + Y --> Z6
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4

f1 = ODEFunction(convert(ODESystem,custom_function_network_1),jac=true)
f2 = ODEFunction(convert(ODESystem,custom_function_network_2),jac=true)
g1 = SDEFunction(convert(SDESystem,custom_function_network_1))
g2 = SDEFunction(convert(SDESystem,custom_function_network_2))
for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
    u0 = factor*rand(rng,length(custom_function_network_1.states))
    p = factor*rand(rng,length(custom_function_network_2.ps))
    t = rand(rng)
    @test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 100*eps())
    @test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p,t)) .< 100*eps())
    @test all(abs.(g1(u0,p,t) .- g2(u0,p,t)) .< 100*eps())
end


text:
using HTTP2
import HTTP2: bytearr
using HTTP2.Frame
using Test
using Dates
using Sockets
using MbedTLS

function sslaccept(server, certfile, keyfile)
    println("Expecting a SSL connection ...")
    sslconfig = MbedTLS.SSLConfig(certfile, keyfile)
    buffer = accept(server)
    sslbuffer = MbedTLS.SSLContext()
    MbedTLS.setup!(sslbuffer, sslconfig)
    MbedTLS.associate!(sslbuffer, buffer)
    MbedTLS.handshake!(sslbuffer)
    return sslbuffer
end

# test serve method
function test_serve(port, body, certfile=nothing, keyfile=nothing)
    server = listen(port)

    println("Server started.")
    connid = 1
    while connid < 4
        println("Waiting for a connection ...")
        buffer = ((certfile === nothing) || (keyfile === nothing)) ? accept(server) : sslaccept(server, certfile, keyfile)
        println("Processing a connection ...")

        connection = HTTP2.Session.new_connection(buffer; isclient=false)
        @info("Connected", connection)
        ## Recv the client preface, and send an empty SETTING frame.

        while true
            headers_evt = HTTP2.Session.take_evt!(connection)
            @info("Received headers", headers_evt)
            if isa(headers_evt, HTTP2.Session.EvtGoaway)
                close(buffer)
                break
            end

            stream_identifier = headers_evt.stream_identifier
            @info("Stream ", stream_identifier)

            sending_headers = [(":status", "200"),
                               ("server", "HTTP2.jl"),
                               ("date", Dates.format(now(Dates.UTC), Dates.RFC1123Format)),
                               ("content-type", "text/html; charset=UTF-8")]
            sending_body = isa(body, String) ? convert(Vector{UInt8}, codeunits(body)) : body
            @info("Resopnding", sending_headers, sending_body)

            HTTP2.Session.put_act!(connection, HTTP2.Session.ActSendHeaders(stream_identifier, sending_headers, false))
            @info("sent headers")
            HTTP2.Session.put_act!(connection, HTTP2.Session.ActSendData(stream_identifier, sending_body, true))
            @info("sent body")
        end
        ## We are done!
        connid += 1
    end
end

# A server example
if length(ARGS) == 2
    certfile = ARGS[1]
    keyfile = ARGS[2]
    test_serve(8000, bytearr("
Hello, world!
"), certfile, keyfile)
else
    test_serve(8000, bytearr("
Hello, world!
"))
end


text:
# CircleMaps

@compat abstract type AbstractCircleMap{D<:PeriodicDomain,R<:PeriodicDomain} <: AbstractMarkovMap{D,R} end

@compat struct FwdCircleMap{D<:PeriodicDomain,R<:PeriodicDomain,ff,gg,T} <: AbstractCircleMap{D,R}
  f::ff
  dfdx::gg
  domain::D
  rangedomain::R
  cover::Int
  fa::T
  fb::T
end
function FwdCircleMap(f::ff,domd::D,randm::R,dfdx::gg=autodiff(f,domd)) where {D<:PeriodicDomain,R<:PeriodicDomain,ff,gg}
  # @assert isempty(∂(randm)) isempty(∂(domd))
  fa = f(leftendpoint(domd)); fb = f(rightendpoint(domd))
  cover_est = (fb-fa)/arclength(randm)
  cover_integer = round(Int,cover_est)
  cover_est ≈ cover_integer || error("Circle map lift does not have integer covering number.")
  FwdCircleMap{D,R,ff,gg,typeof(fa)}(f,dfdx,domd,randm,abs(cover_integer),fa,fb)
end
function FwdCircleMap(f::ff,dom,ran,dfdx::gg=autodiff(f,convert(PeriodicDomain,dom))) where {ff,gg}
  domd = convert(PeriodicDomain,dom); randm = convert(PeriodicDomain,ran)
  FwdCircleMap{typeof(domd),typeof(randm),ff,gg}(f,domd,randm,dfdx)
end

(m::FwdCircleMap)(x) = domain_mod(m.f(x),m.rangedomain)
mapD(m::FwdCircleMap,x) = m.dfdx(x)
mapP(m::FwdCircleMap,x) = (m(x),m.dfdx(x))

mapinv(m::FwdCircleMap,i::Integer,y) = domain_mod(domain_newton(m.f,m.dfdx,interval_mod(y+i*arclength(m.rangedomain),m.fa,m.fb),
    m.domain),m.domain)
mapinvD(m::FwdCircleMap,i::Integer,y) = m.dfdx(mapinv(m,i,y))
function mapinvP(m::FwdCircleMap,i::Integer,y)
  x = mapinv(m,i,y)
  (x,inv(m.dfdx(x)))
end

@compat struct RevCircleMap{D<:Domain,R<:Domain,ff,gg,T} <: AbstractCircleMap{D,R}
  v::ff
  dvdx::gg
  domain::D
  rangedomain::R
  cover::Int
  va::T
  vb::T
end

function RevCircleMap(v::ff,domd::D,randm::R,dvdx::gg=autodiff(v,randm), maxcover=10000) where {D<:Domain,R<:Domain,ff,gg}
  # @assert isempty(∂(ran)) isempty(∂(dom))
  ra = leftendpoint(randm); va = v(ra)
  dr = arclength(randm); dd = arclength(domd)
  cover = 1; vr = v(ra+dr)-va
  while (abs(vr) <= dd && cover < maxcover)
    vr = v(ra+cover*dr)-va
    abs(vr) ≈ dd && break
    abs(vr) > dd && error("Inverse lift doesn't appear to have an inverse")
    cover += 1
  end
  cover == maxcover && error("Can't get to the end of the inverse lift after $maxcover steps")

  RevCircleMap{D,R,ff,gg,typeof(va)}(v,dvdx,domd,randm,cover,va,v(rightendpoint(randm)))
end
function RevCircleMap(v::ff,dom,ran=dom,dvdx::gg=autodiff(v,ran)) where {ff,gg}
  domd = convert(PeriodicDomain,dom); randm = convert(PeriodicDomain,ran)
  RevCircleMap{typeof(domd),typeof(randm),ff,gg}(v,domd,randm,dvdx)
end

mapL(m::RevCircleMap,x) = domain_newton(m.v,m.dvdx,interval_mod(x,m.va,m.vb),m.rangedomain)
(m::RevCircleMap)(x) = domain_mod(mapL(m,x),m.domain)
mapD(m::RevCircleMap,x) = inv(m.dvdx(mapL(m,x)))
function mapP(m::RevCircleMap,x)
  y = mapL(m,x)
  (interval_mod(y,m.domain),inv(m.dvdx(y)))
end

mapinv(m::RevCircleMap,i::Integer,y) = m.v(y+i*arclength(m.rangedomain))
mapinvD(m::RevCircleMap,i::Integer,y) = m.dvdx(y+i*arclength(m.rangedomain))
function mapinvP(m::RevCircleMap,i::Integer,y)
  ys = y+i*arclength(m.rangedomain)
  (m.v(ys),m.dvdx(ys))
end

# autodiff
# FwdCircleMap(f,d,r=d;diff=autodiff(f,d)) = FwdCircleMap(f,diff,d,r);
# RevCircleMap(f,d,r=d;diff=autodiff(f,r)) = RevCircleMap(f,diff,d,r);
CircleMap(f,d,r=d;dir=Forward,diff=autodiff(f,dir==Forward ? d : r)) = dir == Forward ?
  FwdCircleMap(f,d,r,diff) : RevCircleMap(f,d,r,diff)

ncover(m::AbstractCircleMap) = m.cover
eachbranchindex(m::AbstractCircleMap) = 1:m.cover

for TYP in (:FwdCircleMap,:RevCircleMap)
  @eval ApproxFun.domain(m::$TYP) = m.domain
  @eval rangedomain(m::$TYP) = m.rangedomain
end

# Transfer function

function transferfunction(x,m::AbstractCircleMap,f)
  y = zero(eltype(x));
  for b = 1:ncover(m)
    (v,dvdx) = mapinvP(m,b,x)
    y += abs(det(dvdx))*f(v)
  end;
  y
end

function transferfunction_int(x,y,m::AbstractCircleMap,f)
  q = zero(eltype(x));
  csf = cumsum(f)
  for b = 1:ncover(m)
    vy = mapinv(m,b,y); vx = mapinv(m,b,x)
    sgn = sign((vy-vx)/(y-x))
    q += sgn*(csf(vy)-csf(vx))
  end;
  q
end

function transferfunction(x,m::AbstractCircleMap,sk::BasisFun)
  y = zero(eltype(x));
  for b = 1:ncover(m)
    (v,dvdx) = mapinvP(m,b,x)
    y += abs(det(dvdx)).*getbasisfun(v,sk)
  end
  y
end
function transferfunction_int(x,y,m::AbstractCircleMap,sk::BasisFun)
  q = zero(eltype(x));
  for b = 1:ncover(m)
    vy = mapinv(m,b,y); vx = mapinv(m,b,x)
    sgn = sign((vy-vx)/(y-x))
    q += sgn*(getbasisfun_int(vy,sk)-getbasisfun_int(vx,sk))
  end
  q
end


text:
#=
  Copyright 2013 Artur Mkrtchyan

  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License. =#

import Random.randstring
import Base.length

using Memoize

#=
#############################################################

 EntrySpec spells out the specification for an entry .
 For example EntrySpec(:N, 5) is numeric of lenth 5
 and EntrySpec(:C, 16) is alphanumeric of lenth 16

############################################################# =#

const alphabets = Dict(
    :A => ['A':'Z';],
    :N => ['0':'9';],
    :C => vcat(['A':'Z';], ['0':'9';], ['a':'z';])
)

struct EntrySpec
    type::Symbol
    alphabet::Vector{Char}
    length::Int
    EntrySpec(type::Symbol, length::Int) = new(type, alphabets[type], length)
end

_sample(spec::EntrySpec)::String = randstring(spec.alphabet, spec.length)

#=
#############################################################

 EntryType is a classification of IBAN entries as EntrySpecs
 according to domain types, e.g.: BankCode, BranchCode, etc...

############################################################# =#
abstract type EntryType end

struct BankCode <: EntryType
    spec::EntrySpec
end
struct BranchCode <: EntryType
    spec::EntrySpec
end
struct AccountNumber <: EntryType
    spec::EntrySpec
end
struct NationalCheckDigit <: EntryType
    spec::EntrySpec
end
struct AccountType <: EntryType
    spec::EntrySpec
end
struct OwnerAccountType <: EntryType
    spec::EntrySpec
end
struct IdentificationNumber <: EntryType
    spec::EntrySpec
end

stringify(iban_entry) = split(string(typeof(iban_entry.entry_type)), ".")[2]

@memoize _symname(e::EntryType) = Symbol(typeof(e))

#=
#############################################################

BbanStructure is an orderd collection of EntryTypes. Different
countries would have diffrenet structures

############################################################# =#
struct BbanStructure
    entries::Vector{EntryType}
end

##
## Return the length of the BbanStructure in characters
##
@memoize function _length(structure::BbanStructure)::Int
    mapreduce(entry -> entry.spec.length, +, structure.entries)
end

##
## [4 4 12] => [1:4, 5:8, 9:20]
##
@memoize function _slicing(structure::BbanStructure)
    # println("XXX memoize_cache(_slicing): $(length(memoize_cache(_slicing)))")
    result = Tuple{Symbol,UnitRange}[]
    offset = 1

    map(structure.entries) do ent
        symbol = _symname(ent)
        push!(result, (symbol, UnitRange(offset, ent.spec.length + offset - 1)))
        offset += ent.spec.length
    end
    result
end

function parse_bban!(countrycode, partup, allow_random_values=false)

    theiban = TheIban(countrycode)

    bban_structure = for_country(theiban.country_code)
    for entry in bban_structure.entries
        parkey = Symbol(typeof(entry))

        if allow_random_values && ( partup[parkey] === nothing )
            parval = _sample(entry.spec)
        else
            parval = partup[parkey]
        end

        push!(theiban.bban, IbanEntry(parval, entry))
    end

    theiban.bban_str = mapreduce(entry -> entry.value, *, values(theiban.bban))
    theiban
end

for_country(country_code)::Maybe{BbanStructure} =
   get(bban_structures_by_country, country_code, nothing)

"""
    is_supported_country(country_code)::Bool

Return a boolean indicating if the country identified by `country_code` is supported.

# Examples
```julia
julia> is_supported_country("DE")
true

julia> is_supported_country("ZZ")
false
```
"""
is_supported_country(country_code::AbstractString) =
    haskey(bban_structures_by_country, country_code)

"""
    supported_countries()::Array{String,1}

Return an array of supported country codes.
"""
@memoize supported_countries() = collect(keys(bban_structures_by_country))

# *******************
# https://www.swift.com/resource/iban-registry-pdf
# ::Dict{String,BbanStructure}

const bban_structures_by_country = Dict(
    "AD" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:C, 12)),
    ]),
    "AE" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "AL" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 4)),
        NationalCheckDigit(EntrySpec(:N, 1)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "AT" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:N, 11)),
    ]),
    "AZ" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 20)),
    ]),
    "BA" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 8)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "BE" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 7)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "BG" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountType(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:C, 8)),
    ]),
    "BH" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 14)),
    ]),
    "BR" => BbanStructure([
        BankCode(EntrySpec(:N, 8)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:N, 10)),
        AccountType(EntrySpec(:A, 1)),
        OwnerAccountType(EntrySpec(:C, 1)),
    ]),
    "BY" => BbanStructure([
        BankCode(EntrySpec(:C, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "CH" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 12)),
    ]),
    "CR" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 14)),
    ]),
    "CY" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "CZ" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 16)),
    ]),
    "DE" => BbanStructure([
        BankCode(EntrySpec(:N, 8)),
        AccountNumber(EntrySpec(:N, 10)),
    ]),
    "DK" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 10)),
    ]),
    "DO" => BbanStructure([
        BankCode(EntrySpec(:C, 4)),
        AccountNumber(EntrySpec(:N, 20)),
    ]),
    "EE" => BbanStructure([
        BankCode(EntrySpec(:N, 2)),
        BranchCode(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:N, 11)),
        NationalCheckDigit(EntrySpec(:N, 1)),
    ]),
    "ES" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        BranchCode(EntrySpec(:N, 4)),
        NationalCheckDigit(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:N, 10)),
    ]),
    "FI" => BbanStructure([
        BankCode(EntrySpec(:N, 6)),
        AccountNumber(EntrySpec(:N, 7)),
        NationalCheckDigit(EntrySpec(:N, 1)),
    ]),
    "FO" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 9)),
        NationalCheckDigit(EntrySpec(:N, 1)),
    ]),
    "FR" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 11)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "GB" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        BranchCode(EntrySpec(:N, 6)),
        AccountNumber(EntrySpec(:N, 8)),
    ]),
    "GE" => BbanStructure([
        BankCode(EntrySpec(:A, 2)),
        AccountNumber(EntrySpec(:N, 16)),
    ]),
    "GI" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 15)),
    ]),
    "GL" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 10)),
    ]),
    "GR" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "GT" => BbanStructure([
        BankCode(EntrySpec(:C, 4)),
        AccountNumber(EntrySpec(:C, 20)),
    ]),
    "HR" => BbanStructure([
        BankCode(EntrySpec(:N, 7)),
        AccountNumber(EntrySpec(:N, 10)),
    ]),
    "HU" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 16)),
        NationalCheckDigit(EntrySpec(:N, 1)),
    ]),
    "IE" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        BranchCode(EntrySpec(:N, 6)),
        AccountNumber(EntrySpec(:N, 8)),
    ]),
    "IL" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 13)),
    ]),
    "IR" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 19)),
    ]),
    "IS" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        BranchCode(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:N, 6)),
        IdentificationNumber(EntrySpec(:N, 10)),
    ]),
    "IT" => BbanStructure([
        NationalCheckDigit(EntrySpec(:A, 1)),
        BankCode(EntrySpec(:N, 5)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 12)),
    ]),
    "JO" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:C, 18)),
    ]),
    "KW" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 22)),
    ]),
    "KZ" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:C, 13)),
    ]),
    "LB" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:C, 20)),
    ]),
    "LC" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 24)),
    ]),
    "LI" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 12)),
    ]),
    "LT" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:N, 11)),
    ]),
    "LU" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:C, 13)),
    ]),
    "LV" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 13)),
    ]),
    "MC" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 11)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "MD" => BbanStructure([
        BankCode(EntrySpec(:C, 2)),
        AccountNumber(EntrySpec(:C, 18)),
    ]),
    "ME" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 13)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "MK" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:C, 10)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "MR" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:N, 11)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "MT" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 18)),
    ]),
    "MU" => BbanStructure([
        BankCode(EntrySpec(:C, 6)),
        BranchCode(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:C, 18)),
    ]),
    "NL" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:N, 10)),
    ]),
    "NO" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 6)),
        NationalCheckDigit(EntrySpec(:N, 1)),
    ]),
    "PK" => BbanStructure([
        BankCode(EntrySpec(:C, 4)),
        AccountNumber(EntrySpec(:N, 16)),
    ]),
    "PL" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        BranchCode(EntrySpec(:N, 4)),
        NationalCheckDigit(EntrySpec(:N, 1)),
        AccountNumber(EntrySpec(:N, 16)),
    ]),
    "PS" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 21)),
    ]),
    "PT" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 11)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "QA" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 21)),
    ]),
    "RO" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "RS" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 13)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "SA" => BbanStructure([
        BankCode(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:C, 18)),
    ]),
    "SC" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 16)),
        AccountType(EntrySpec(:A, 3)),
    ]),
    "SE" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 17)),
    ]),
    "SI" => BbanStructure([
        BankCode(EntrySpec(:N, 2)),
        BranchCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 8)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "SK" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 16)),
    ]),
    "SM" => BbanStructure([
        NationalCheckDigit(EntrySpec(:A, 1)),
        BankCode(EntrySpec(:N, 5)),
        BranchCode(EntrySpec(:N, 5)),
        AccountNumber(EntrySpec(:C, 12)),
    ]),
    "ST" => BbanStructure([
        BankCode(EntrySpec(:N, 4)),
        BranchCode(EntrySpec(:N, 4)),
        AccountNumber(EntrySpec(:N, 13)),
    ]),
    "SV" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:N, 20)),
    ]),
    "TL" => BbanStructure([
        BankCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:N, 14)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
    "TN" => BbanStructure([
        BankCode(EntrySpec(:N, 2)),
        BranchCode(EntrySpec(:N, 3)),
        AccountNumber(EntrySpec(:C, 15)),
    ]),
    "TR" => BbanStructure([
        BankCode(EntrySpec(:N, 5)),
        NationalCheckDigit(EntrySpec(:C, 1)),
        AccountNumber(EntrySpec(:C, 16)),
    ]),
    "UA" => BbanStructure([
        BankCode(EntrySpec(:N, 6)),
        AccountNumber(EntrySpec(:N, 19)),
    ]),
    "VG" => BbanStructure([
        BankCode(EntrySpec(:A, 4)),
        AccountNumber(EntrySpec(:N, 16)),
    ]),
    "XK" => BbanStructure([
        BankCode(EntrySpec(:N, 2)),
        BranchCode(EntrySpec(:N, 2)),
        AccountNumber(EntrySpec(:N, 10)),
        NationalCheckDigit(EntrySpec(:N, 2)),
    ]),
)


text:
include("input.jl")

o2input = split(input, "\n")
c2input = split(input, "\n")

for bit in 1:length(o2input[1])
    num0 = 0
    num1 = 0
    for byte in o2input
        if byte[bit] == '1'
            num1+=1
        else
            num0+=1
        end
    end
    thing = string(Int(num1>=num0))
    if length(o2input) == 1
        break
    end
    for byte in length(o2input):-1:1
        if string(o2input[byte][bit]) != thing
            deleteat!(o2input, byte)
        end
    end
end

for bit in 1:length(c2input[1])
    num0 = 0
    num1 = 0
    for byte in c2input
        if byte[bit] == '1'
            num1+=1
        else
            num0+=1
        end
    end
    thing = string(Int(num0>num1))
    if length(c2input) == 1
        break
    end
    for byte in length(c2input):-1:1
        if string(c2input[byte][bit]) != thing
            deleteat!(c2input, byte)
        end
    end
end

println(parse(Int,o2input[1],base=2)*parse(Int,c2input[1],base=2))


text:
import FLOWNoise
using Test
include("testingfunctions.jl")

@testset "Post-processing: Time to frequency domain" begin
    @testset "Case 1" begin # currently only converts time-domain output to frequency-domain output and compares to WOPWOP frequency domain

        filepath = "./data/case1/"
        windowtype = "None"

        thickness_dB_calc, loading_dB_calc, total_dB_calc,
            thickness_dBA_calc, loading_dBA_calc, total_dBA_calc,
            thickness_dB_actual, loading_dB_actual, total_dB_actual,
            thickness_dBA_actual, loading_dBA_actual, total_dBA_actual,
            thickness_OASPL_calc, loading_OASPL_calc, total_OASPL_calc,
            thickness_OASPLA_calc, loading_OASPLA_calc, total_OASPLA_calc,
            thickness_OASPL_actual, loading_OASPL_actual, total_OASPL_actual,
            thickness_OASPLA_actual, loading_OASPLA_actual, total_OASPLA_actual = prepare_test(filepath, windowtype)

        tol_case1 = 1e-2 # 1% error

        @testset "OASPL" begin
            @test isapprox(percent_error(thickness_OASPL_actual[1], thickness_OASPL_calc), 0.0, atol=tol_case1)
            @test isapprox(percent_error(loading_OASPL_actual[1], loading_OASPL_calc), 0.0, atol=tol_case1)
            @test isapprox(percent_error(total_OASPL_actual[1], total_OASPL_calc), 0.0, atol=tol_case1)
        end
        @testset "OASPLA" begin
            @test isapprox(percent_error(thickness_OASPLA_actual[1], thickness_OASPLA_calc), 0.0, atol=tol_case1)
            @test isapprox(percent_error(loading_OASPLA_actual[1], loading_OASPLA_calc), 0.0, atol=tol_case1)
            @test isapprox(percent_error(total_OASPLA_actual[1], total_OASPLA_calc), 0.0, atol=tol_case1)
        end
        @testset "Thickness SPL" begin #for loops are split up to make the tests more distinct, easier to see what's wrong instead of lumping together
            for i in 1:length(thickness_dB_actual)
                @test isapprox(percent_error(thickness_dB_actual[i], thickness_dB_calc[i]), 0.0, atol=tol_case1)
            end
        end
        @testset "Loading SPL" begin
            for i in 1:length(thickness_dB_actual)
                @test isapprox(percent_error(loading_dB_actual[i], loading_dB_calc[i]), 0.0, atol=tol_case1)
            end
        end
        @testset "Total SPL" begin
            for i in 1:length(thickness_dB_actual)
                @test isapprox(percent_error(total_dB_actual[i], total_dB_calc[i]), 0.0, atol=tol_case1)
            end
        end
        @testset "Thickness SPLA" begin
            for i in 1:length(thickness_dBA_actual)
                @test isapprox(percent_error(thickness_dBA_actual[i], thickness_dBA_calc[i]), 0.0, atol=tol_case1)
            end
        end
        @testset "Loading SPLA" begin
            for i in 1:length(thickness_dBA_actual)
                @test isapprox(percent_error(loading_dBA_actual[i], loading_dBA_calc[i]), 0.0, atol=tol_case1)
            end
        end
        @testset "Total SPLA" begin
            for i in 1:length(thickness_dBA_actual)
                @test isapprox(percent_error(total_dBA_actual[i], total_dBA_calc[i]), 0.0, atol=tol_case1)
            end
        end
    end

    # @testset "Case 2" begin #- Look at low/high pass filtering

    #     filepath = "./data/case2/"
    #     windowtype = "Hann"

    #     thickness_OASPL_calc, loading_OASPL_calc, total_OASPL_calc,
    #         thickness_OASPLA_calc, loading_OASPLA_calc, total_OASPLA_calc,
    #         thickness_OASPL_actual, loading_OASPL_actual, total_OASPL_actual,
    #         thickness_OASPLA_actual, loading_OASPLA_actual, total_OASPLA_actual = prepare_test(filepath, windowtype)

    #     tol_case1 = 1e-2

    #     @test isapprox(OASPL_thickness_dB, thickness_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_loading_dB, loading_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_total_dB, total_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_thickness_dBA, thickness_OASPLA[1], atol=tol_case1)
    #     @test isapprox(OASPL_loading_dBA, loading_OASPLA[1], atol=tol_case1)
    #     @test isapprox(OASPL_total_dBA, total_OASPLA[1], atol=tol_case1)

    # end

    # @testset "Case 4" begin - need to look at filtering

    # filepath = "./data/case5/"
    #     windowtype = "None"

    #     thickness_OASPL_calc, loading_OASPL_calc, total_OASPL_calc,
    #         thickness_OASPLA_calc, loading_OASPLA_calc, total_OASPLA_calc,
    #         thickness_OASPL_actual, loading_OASPL_actual, total_OASPL_actual,
    #         thickness_OASPLA_actual, loading_OASPLA_actual, total_OASPLA_actual = prepare_test(filepath, windowtype)

    #     tol_case1 = 1e-2

    #     @test isapprox(OASPL_thickness_dB, thickness_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_loading_dB, loading_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_total_dB, total_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_thickness_dBA, thickness_OASPLA[1], atol=tol_case1)
    #     @test isapprox(OASPL_loading_dBA, loading_OASPLA[1], atol=tol_case1)
    #     @test isapprox(OASPL_total_dBA, total_OASPLA[1], atol=tol_case1)

    # end

    # @testset "Case 5" begin - SPL found only for two specific ranges, not each frequency in the spectrum

    # filepath = "./data/case5/"
    #     windowtype = "None"

    #     thickness_OASPL_calc, loading_OASPL_calc, total_OASPL_calc,
    #         thickness_OASPLA_calc, loading_OASPLA_calc, total_OASPLA_calc,
    #         thickness_OASPL_actual, loading_OASPL_actual, total_OASPL_actual,
    #         thickness_OASPLA_actual, loading_OASPLA_actual, total_OASPLA_actual = prepare_test(filepath, windowtype)

    #     tol_case1 = 1e-2

    #     @test isapprox(OASPL_thickness_dB, thickness_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_loading_dB, loading_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_total_dB, total_OASPL[1], atol=tol_case1)
    #     @test isapprox(OASPL_thickness_dBA, thickness_OASPLA[1], atol=tol_case1)
    #     @test isapprox(OASPL_loading_dBA, loading_OASPLA[1], atol=tol_case1)
    #     @test isapprox(OASPL_total_dBA, total_OASPLA[1], atol=tol_case1)

    # end

end

max_stars_repo_name: byuflowlab/FLOWNoise

max_issues_repo_name: byuflowlab/FLOWNoise

max_forks_repo_name: byuflowlab/FLOWNoise


text:
"""
	doc(fname, title="", width=800, height=600)

Display the contents of file name `fname` in a browser (webview) window.
Alternatively, `fname` can be a URL address.

Example, display the contents of the README.md file located in current dir.

   doc("README.md")

to display the contents of the `examples.md` file located in the docs/src dir of the module GMT, do

   doc(GMT, "examples.md")

To display the doc strings of any documented function, just do (i.e to show the *sin* docs)

   doc(sin)
"""
function doc(fname::AbstractString; title::AbstractString="", width::Integer=0, height::Integer=0)

	if (width == 0)  width  = 800  end
	if (height == 0) height = 600  end
	if (title == "") title, = splitext(basename(fname))  end
	if (startswith(fname, "http") || startswith(fname, "ftp") || startswith(fname, "file"))
		webview(title, fname, width, height)
	else
		MD = Markdown.parse_file(fname);
		webview(title, "data:text/html," * html(MD) * "", width, height)
	end
end

function doc(mod::Module, fname::AbstractString; title::AbstractString="", width::Integer=0, height::Integer=0)
	if (fname == "README" || fname == "README.md")
		fn = joinpath(dirname(pathof(mod)), "..", "README.md")
	else
		fn = joinpath(dirname(pathof(mod)), "..", joinpath("docs/src", fname))
	end
	doc(fn, title=title, width=width, height=height)
end

function doc(fun; width::Integer=0, height::Integer=0)
	# Method to show the @doc string of a function
	if (width == 0)  width  = 800  end
	if (height == 0) height = 600  end
	webview(string(fun), "data:text/html," * html(Base.Docs.doc(fun)) * "", width, height)
end


text:
function ImageKeypointRegression(
        blocks::Tuple{<:Image{N},<:Keypoints{N}},
        data=nothing;
        size=ntuple(i -> 128, N),
		aug_projections=DataAugmentation.Identity(),
		aug_image=DataAugmentation.Identity(),
        C=RGB{N0f8},
        computestats=false,
    ) where N
    return SupervisedTask(
        blocks,
        (
			ProjectiveTransforms(size; augmentations=aug_projections),
			getimagepreprocessing(data, computestats; C=C, augmentations=aug_image),
            KeypointPreprocessing(size),
        )
    )
end

"""
    ImageKeypointRegression(size, nkeypoints; kwargs...)

Learning task for regressing a set of `nkeypoints` keypoints from
images. Images are resized to `size` and a class is predicted for every pixel.
"""
function ImageKeypointRegression(size::NTuple{N,Int}, nkeypoints::Int; kwargs...) where N
    blocks = (Image{N}(), Keypoints{N}((nkeypoints,)))
    return ImageKeypointRegression(blocks; size = size, kwargs...)
end

registerlearningtask!(FASTAI_METHOD_REGISTRY, ImageKeypointRegression, (Image, Keypoints))

# ## Tests

@testset "ImageKeypointRegression [task]" begin
    task = ImageKeypointRegression((16, 16), 10)
    FastAI.checktask_core(task)
    @testset "Show backends" begin
        @testset "ShowText" begin
            FastAI.test_task_show(task, ShowText(Base.DevNull()))
        end
    end
end


text:
#=

  Thermostat problem.

  OMEGA: Fast, casual Inference from Simple Parts
  From  https://www.youtube.com/watch?v=oCvbqKE2tWA
  @ ~ 21min

  The Omega (Julia) model is here:
  https://github.com/zenna/Omega.jl/blob/master/OmegaModels/src/causal/thermostat.jl

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

            time    1.9912    0.8130     0.0081    0.0108   6852.7308    0.9999      574.6525
  is_window_open    0.4929    0.5000     0.0050    0.0065   6729.5331    1.0001      564.3214
        is_ac_on    0.2550    0.4359     0.0044    0.0056   6931.3351    1.0000      581.2440
    outside_temp   20.7561    9.2221     0.0922    0.1171   5860.1765    1.0002      491.4194
       room_temp   23.6955    2.9221     0.0292    0.0377   7120.8508    0.9999      597.1363
      thermostat   22.6500    4.0236     0.0402    0.0527   6486.0439    1.0000      543.9031

   Distributions of variable time
   afternoon  =>    3390  (0.339000)
   morning    =>    3349  (0.334900)
   evening    =>    3261  (0.326100)

   Distributions of variable is_window_open (num:0)
   0.00000 =>    5071  (0.507100)
   1.00000 =>    4929  (0.492900)

   Distributions of variable is_ac_on (num:0)
   0.00000 =>    7450  (0.745000)
   1.00000 =>    2550  (0.255000)

  Cf ~/cplint/thermostat.pl
     ~/blog/thermostat.blog
     ~/webppl/thermostat.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function thermostat()
    morning = 1
    afternoon = 2
    evening = 3
    times = [morning, afternoon,evening];

    time ~ DiscreteUniform(morning,evening)

    is_window_open ~ flip(0.5)

    # a.c. is off when window is closed
    is_ac_on ~ is_window_open ? Dirac(false) : flip(0.5)

    # hottest at noon, cool at night
    outside_temp ~
        (time == morning)   ? Normal(20,2) :
        (time == afternoon) ? Normal(32,2) :
        (time == evening)   ? Normal(10,2) : Normal(25,2)

    # a.c. chills the room
    room_temp ~ is_ac_on ? Normal(20,2) : Normal(25,2)

    # great insulation
    thermostat ~ Dirac(is_window_open ? (outside_temp + room_temp) / 2.0 : room_temp)

    # true ~ Dirac(is_ac_on==true)
    # true ~ Dirac(is_window_open == false)
    # true ~ Dirac(room_temp < 20.0)
    # true ~ Dirac(room_temp > 20.0)
    # true ~ Dirac(time == evening)
    # true ~ Dirac(outside_temp > 10.0)

end

model = thermostat()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:time,["morning","afternoon","evening"])
show_var_dist_pct(chns,:is_window_open)
show_var_dist_pct(chns,:is_ac_on)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
push!(LOAD_PATH, dirname(@__DIR__))

using Documenter
using StructArrays

makedocs(
    sitename = "StructArrays",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [StructArrays]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/JuliaArrays/StructArrays.jl.git",
    push_preview = true,
)


text:
# This file was generated, do not modify it. # hide
for _ in 1:10^4
    push!(G, rand(Float32, dim))  # push! inserts one item at a time
end

append!(G, MatrixDatabase(rand(Float32, dim, 10^4))) # append! inserts many items at once

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
# Discrete Uniform Stein Distribution
#
# Represents a uniform distribution with p independent components taking on
# integer values between 0 and ulim

type SteinDiscreteUniform <: SteinDistribution
    # Marginals are DUnif(0,ulim)
    ulim::Int64
    # Dimensionality
    p::Int64
    # Stein factors
    c1::Float64
    c2::Float64
    c3::Float64
end

# Constructor with (c1,c2,c3) = (1,1,1)
SteinDiscreteUniform(ulim::Int64, p::Int64) = SteinDiscreteUniform(
    ulim,
    p,
    1.0,
    1.0,
    1.0
);

# Draw n independent samples from distribution
function rand(d::SteinDiscreteUniform, n::Int64)
    p = d.p;
    # Sample DUnif(0,ulim) variables
    rand(0:d.ulim, n, p)
end

# Joint distribution evaluated at a point y (unnormalized)
function jointdistribution(d::SteinDiscreteUniform, y::Array{Float64})
    1.0
end

# Conditional distribution of jth coordinate of y, given all other coordinates,
# represented as a probability vector
function conddistribution(d::SteinDiscreteUniform, y::Array{Float64,1}, j::Int64)
    ones(d.ulim+1) / (d.ulim+1)
end

function supportlowerbound(d::SteinDiscreteUniform, j::Int64)
    0
end

function supportupperbound(d::SteinDiscreteUniform, j::Int64)
    d.ulim
end

# Cumulative distribution function (only valid when X is univariate)
function cdf(d::SteinDiscreteUniform, t)
    if t >= d.ulim
        return 1.0
    elseif t < 0
        return 0.0
    else
        return (floor(t)+1.0)/(d.ulim+1.0)
    end
end

function numdimensions(d::SteinDiscreteUniform)
    d.p
end


text:
"""
    namefilter(ids::Vector{T}) where {T <: Integer}

Returns a subset of the names table where only the given taxids are present.
"""
function namefilter(ids::Vector{T}) where {T <: Integer}
    return leftjoin(DataFrame(tax_id=ids), NCBITaxonomy.names_table; on=:tax_id)
end

"""
    namefilter(taxa::Vector{NCBITaxon})

Returns a subset of the names table dataset, where the taxids of the taxa are
present. This includes all names, not only the scientific names.
"""
function namefilter(taxa::Vector{NCBITaxon})
    return namefilter([t.id for t in taxa])
end

"""
    namefilter(division::Symbol)

Returns a subset of the names table for all names under a given NCBI division.
"""
function namefilter(division::Symbol)
    ids = findall(isequal(division), NCBITaxonomy.nodes_table.division_code)
    return namefilter(NCBITaxonomy.nodes_table.tax_id[ids])
end

"""
    namefilter(division::Vector{Symbol})

Returns a subset of the names table for all names under a number of multiple
NCBI divisions.
"""
function namefilter(division::Vector{Symbol})
    ids = findall(x -> x in division, NCBITaxonomy.nodes_table.division_code)
    return namefilter(NCBITaxonomy.nodes_table.tax_id[ids])
end


text:
using MicroLogging

@info "# Logging macros"
@debug "A message for debugging (filtered out by default)"
@info "Information about normal program operation"
@warn "A potential problem was detected"
@error "Something definitely went wrong"
x = [1 2;3 4]
@info "Support for key value pairs" x a=1 b="asdf"

#-------------------------------------------------------------------------------
@info "# Progress logging"
for i=1:100
    sleep(0.01)
    @info "algorithm1" progress=i/100
end

#-------------------------------------------------------------------------------
@info "# Log record filtering"
@debug begin
    error("Should not be executed unless logging at debug level")
    "A message"
end
configure_logging(min_level=:debug)
@debug "Logging enabled at debug level and above"
for i=1:10
    @warn "Log suppression iteration $i (max_log=2)" max_log=2
end
module LogTest
    using MicroLogging
    function f()
        @debug "Message from LogTest"
        @info  "Message from LogTest"
        @warn  "Message from LogTest"
        @error "Message from LogTest"
    end
end
LogTest.f()
configure_logging(LogTest, min_level=:error)
@info "Set log filtering to error level for LogTest module"
LogTest.f()

#-------------------------------------------------------------------------------
@info "# Task-based log dispatch using dynamic scoping"
function some_operation()
    @info "Dispatches to the current task logger, or the global logger"
end
logstream = IOBuffer()
with_logger(SimpleLogger(logstream)) do
    @info "Logging redirected"
    some_operation()
end
@info "Logs, captured separately in the with_logger() block" logstring=strip(String(take!(logstream)))

#-------------------------------------------------------------------------------
@info "# Formatting logs can't crash the application"
@info "Blah $(error("An intentional error"))"

configure_logging(min_level=:info)
nothing


text:
function Generate(M::LGSSM, n::Int, N::Int)
    ydim, qdim = size(M.B)
    xdim, pdim = size(M.R)
    L = cholfact(M.Σᵥ)[:L]
    # RR = M.R*M.R'
    # SS = M.S*M.S'
    X = zeros(xdim, n, N)
    Y = zeros(ydim, n, N)
    for t in 1:n
        # Generate States
        if t == 1
            X[:, 1, :] = L*randn(xdim, N)
        else
            X[:,t,:] = M.A*X[:,t-1,:] + M.R*randn(pdim ,N)
        end

        # Generate Measurements
        Y[:,t,:] = M.B*X[:,t,:] + M.S*randn(ydim, N)
    end
    return Dict(:X => X, :Y => Y)
end


text:
"""
    A = JopLaplacian(sp)

Build a 2D (d^2/dx^2+d^2/dy^2) or 3D (d^2/dx^2+d^2/dy^2+d^2/dz^2) Laplacian
operator with three point centered finite difference stencil, with
`sp::JetSpace` and `A::JopLaplacian`.

For example,
```julia
A = JopLaplacian(JetSpace(Float64,128,256))
m = rand(domain(A))
d = A*m
```
"""
JopLaplacian(spc::JetSpace) = JopLn(dom = spc, rng = spc, df! = JopLaplacian_df!, df′! = JopLaplacian_df′!)

export JopLaplacian

function JopLaplacian_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; kwargs...) where {T}
    nz,nx=size(m)
    d .= 0
    for ix=2:nx-1, iz=2:nz-1
        d[iz,ix] -= 4 * m[iz+0,ix+0]
        d[iz,ix] += m[iz-1,ix+0]
        d[iz,ix] += m[iz+1,ix+0]
        d[iz,ix] += m[iz+0,ix-1]
        d[iz,ix] += m[iz+0,ix+1]
    end
    d
end

function JopLaplacian_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; kwargs...) where {T}
    nz,nx=size(m)
    m .= 0
    for ix=2:nx-1, iz=2:nz-1
        m[iz+0,ix+0] -= 4 * d[iz,ix]
        m[iz-1,ix+0] += d[iz,ix]
        m[iz+1,ix+0] += d[iz,ix]
        m[iz+0,ix-1] += d[iz,ix]
        m[iz+0,ix+1] += d[iz,ix]
    end
    m
end

function JopLaplacian_df!(d::AbstractArray{T,3}, m::AbstractArray{T,3}; kwargs...) where {T}
    nz,ny,nx=size(m)
    d .= 0
    for ix=2:nx-1,iy=2:ny-1,iz=2:nz-1
        d[iz,iy,ix] -= 6 * m[iz,iy,ix]
        d[iz,iy,ix] += m[iz-1,iy+0,ix+0]
        d[iz,iy,ix] += m[iz+1,iy+0,ix+0]
        d[iz,iy,ix] += m[iz+0,iy-1,ix+0]
        d[iz,iy,ix] += m[iz+0,iy+1,ix+0]
        d[iz,iy,ix] += m[iz+0,iy+0,ix-1]
        d[iz,iy,ix] += m[iz+0,iy+0,ix+1]
    end
    d
end

function JopLaplacian_df′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; kwargs...) where {T}
    nz,ny,nx=size(m)
    m .= 0
    for ix=2:nx-1,iy=2:ny-1,iz=2:nz-1
        m[iz+0,iy+0,ix+0] -= 6 * d[iz,iy,ix]
        m[iz-1,iy+0,ix+0] += d[iz,iy,ix]
        m[iz+1,iy+0,ix+0] += d[iz,iy,ix]
        m[iz+0,iy-1,ix+0] += d[iz,iy,ix]
        m[iz+0,iy+1,ix+0] += d[iz,iy,ix]
        m[iz+0,iy+0,ix-1] += d[iz,iy,ix]
        m[iz+0,iy+0,ix+1] += d[iz,iy,ix]
    end
    m
end