text:
struct Point{Dim}
    coords::Vec{Dim,Float}
    Point(vec::Vec{Dim,T})  where {Dim,T} = new{Dim}(Vec{Dim,Float}(vec))
end

Point(vec::AbstractVector) = Point(Vec{length(vec)}(vec))
Point(coords...) = Point(Vec(coords...))
Point(p::Point) = p
coordinates(p::Point) = p.coords

# point  + vector ---> point
Base.:(+)(p::Point{Dim},v::Vec{Dim}) where Dim = Point(coordinates(p) + v)
# vector + point ---> point
Base.:(+)(v::Vec{Dim},p::Point{Dim}) where Dim =   Point(p + v)
# point  - point ---> vector
Base.:(-)(a::Point{Dim},b::Point{Dim}) where Dim = Vec(coordinates(a) - coordinates(b))
# point  - vector ---> point
Base.:(-)(a::Point{Dim},v::Vec{Dim}) where Dim = Point(coordinates(a) - v)
# vector - point ---> point
Base.:(-)(v::Vec{Dim},p::Point{Dim}) where Dim = Point(p - v)


text:
include("NKOM.jl")
using .NKOM

zs = - collect(Float64, range(0.0, stop=1000.0, length=101))

mid_zs = (zs[1:end-1] + zs[2:end]) / 2.0

Δhs = zs[1:end-1] - zs[2:end]
Δzs = (Δhs[1:end-1] + Δhs[2:end]) / 2.0

qs = zeros(Float64, length(zs)-1)

for k = 1:length(qs)
    qs[k] = exp(- ((mid_zs[k] + 500.0) / 100.0)^2.0 )
end

for k = 1:length(qs)

    if mid_zs[k] > -50.0
        qs[k] = 1.0
    elseif qs[k] > -500.0
        qs[k] = 1.0 - (-50 - mid_zs[k]) * 0.0001
    else
        qs[k] = qs[k-1]
    end

end

#qs .= 1

ws = copy(zs) * 0.0 .+ 1e-4
ws[1] = 0.0

Δt = 1e4 # 3600.0 * 3

Nz = length(qs)

qstmp  = copy(qs)
flxtmp = copy(qs)

steps = 5

qs_rec = zeros(Float64, length(qs), steps + 1)
qs_rec[:, 1] = qs

for t = 1:steps

    NKOM.doZAdvection_MacCormack!(
        Nz = Nz,
        qs = qs,
        ws = ws,
        Δzs = Δzs,
        Δt = Δt,
        qstmp = qstmp,
        flxtmp=flxtmp,
    )

    qs_rec[:, t+1] = qs
    println(t, ": ", sum(qs .* Δhs))
end

using PyCall

plt = pyimport("matplotlib.pyplot")

fig, ax = plt.subplots(1, 1)

for t = 1:1:steps+1

    c = ["r", "g", "b", "k", "gray", "purple"][mod(t-1, 5)+1]
    for k = 1:Nz

        ax.plot([ qs_rec[k, t], qs_rec[k,   t]], [zs[k], zs[k+1]] , color=c)

        if k < Nz
            ax.plot([ qs_rec[k, t], qs_rec[k+1, t]], [zs[k+1], zs[k+1]] , color=c)
        end

    end

end

plt.show()

ext: jl

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function getmixturemodel(Q, mu, C)

  MixtureModel([MvNormal(mu[kk,:],C[kk]) for kk=1:length(Q)], Q)

end


text:
function _handle_base_mva(pwf_data::Dict)
    baseMVA = 100.0 # Default value for this field in .pwf
    if haskey(pwf_data, "DCTE")
        if haskey(pwf_data["DCTE"], "BASE")
            baseMVA = pwf_data["DCTE"]["BASE"]
        end
    end
    return baseMVA
end

function _create_dict_dctr(data::Dict)
    dctr_dict = Dict{Tuple, Any}()

    for (k,v) in data
        sub_data = Dict{String, Any}()

        voltage_control = haskey(v, "VOLTAGE CONTROL") ? v["VOLTAGE CONTROL"] : true
        phase_control = haskey(v, "PHASE CONTROL") ? v["PHASE CONTROL"] : true
        if voltage_control && !phase_control
            sub_data["TYPE OF CONTROL"] = "VOLTAGE CONTROL"
            sub_data["MINIMUM VOLTAGE"] = v["MINIMUM VOLTAGE"]
            sub_data["MAXIMUM VOLTAGE"] = v["MAXIMUM VOLTAGE"]
            sub_data["SPECIFIED VALUE"] = v["SPECIFIED VALUE"]
            sub_data["MEASUREMENT EXTREMITY"] = v["MEASUREMENT EXTREMITY"]
        elseif phase_control && !voltage_control
            sub_data["TYPE OF CONTROL"] = "PHASE CONTROL"
            sub_data["MINIMUM PHASE"] = v["MINIMUM PHASE"]
            sub_data["MAXIMUM PHASE"] = v["MAXIMUM PHASE"]
            sub_data["SPECIFIED VALUE"] = v["SPECIFIED VALUE"]
            sub_data["MEASUREMENT EXTREMITY"] = v["MEASUREMENT EXTREMITY"]
            sub_data["CONTROL TYPE"] = v["CONTROL TYPE"]
        else
            @error("DCTR data should refer either to voltage or phase control")
        end
        current_value = get(dctr_dict, (v["FROM BUS"], v["TO BUS"]), Dict{Int, Any}())
        current_value[v["CIRCUIT"]] = sub_data
        dctr_dict[(v["FROM BUS"], v["TO BUS"])] = current_value
    end
    return dctr_dict
end

function _pwf2pm_transformer!(pm_data::Dict, pwf_data::Dict, branch::Dict; add_control_data::Bool=false) # Two-winding transformer
    sub_data = Dict{String,Any}()

    sub_data["f_bus"] = pop!(branch, "FROM BUS")
    sub_data["t_bus"] = pop!(branch, "TO BUS")
    sub_data["br_r"] = pop!(branch, "RESISTANCE") / 100
    sub_data["br_x"] = pop!(branch, "REACTANCE") / 100
    sub_data["g_fr"] = 0.0
    sub_data["g_to"] = 0.0
    sub_data["tap"] = pop!(branch, "TAP")
    sub_data["shift"] = -pop!(branch, "LAG")
    sub_data["angmin"] = -60.0 # PowerModels.jl standard
    sub_data["angmax"] = 60.0 # PowerModels.jl standard
    sub_data["transformer"] = true

    if branch["STATUS"] == branch["OPENING FROM BUS"] == branch["OPENING TO BUS"] == 'L'
        sub_data["br_status"] = 1
    else
        sub_data["br_status"] = 0
    end

    n = 0 # count(x -> x["f_bus"] == sub_data["f_bus"] && x["t_bus"] == sub_data["t_bus"], values(pm_data["branch"]))
    sub_data["source_id"] = ["transformer", sub_data["f_bus"], sub_data["t_bus"], 0, "0$(n + 1)", 0]
    sub_data["index"] = length(pm_data["branch"]) + 1

    dict_dshl = haskey(pwf_data, "DSHL") ? _create_dict_dshl(pwf_data["DSHL"]) : nothing
    b = _handle_b(pwf_data, sub_data["f_bus"], sub_data["t_bus"], branch["SHUNT SUSCEPTANCE"], branch["CIRCUIT"], dict_dshl)
    sub_data["b_fr"] = b[1]
    sub_data["b_to"] = b[2]

    sub_data["rate_a"] = pop!(branch, "NORMAL CAPACITY")
    sub_data["rate_b"] = pop!(branch, "EMERGENCY CAPACITY")
    sub_data["rate_c"] = pop!(branch, "EQUIPAMENT CAPACITY")

    if sub_data["rate_a"] >= 9999
        delete!(sub_data, "rate_a")
    end
    if sub_data["rate_b"] >= 9999
        delete!(sub_data, "rate_b")
    end
    if sub_data["rate_c"] >= 9999
        delete!(sub_data, "rate_c")
    end

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()
        sub_data["control_data"]["tapmin"] = pop!(branch, "MINIMUM TAP")
        sub_data["control_data"]["tapmax"] = pop!(branch, "MAXIMUM TAP")
        sub_data["control_data"]["circuit"] = branch["CIRCUIT"]

        sub_data["control_data"]["controlled_bus"] = abs(branch["CONTROLLED BUS"])

        dict_dctr = haskey(pwf_data, "DCTR") ? _create_dict_dctr(pwf_data["DCTR"]) : nothing
        constraint_type = "fix"
        if isa(dict_dctr, Dict) && haskey(dict_dctr, (sub_data["f_bus"], sub_data["t_bus"]))
            branch_dctr = dict_dctr[(sub_data["f_bus"], sub_data["t_bus"])]
            circuit = branch["CIRCUIT"]
            if haskey(branch_dctr, circuit)
                constraint_type = branch_dctr[circuit]["TYPE OF CONTROL"]
            end
        end

        tap_bounds = sub_data["control_data"]["tapmin"] !== nothing && sub_data["control_data"]["tapmax"] !== nothing
        if tap_bounds # tap control
            sub_data["control_data"]["control_type"] = "tap_control"
            sub_data["control_data"]["shift_control_variable"] = nothing
            sub_data["control_data"]["shiftmin"] = nothing
            sub_data["control_data"]["shiftmax"] = nothing
            if constraint_type == "VOLTAGE CONTROL"
                sub_data["control_data"]["constraint_type"] = "bounds"
                sub_data["control_data"]["valsp"] = branch_dctr[circuit]["SPECIFIED VALUE"]
                sub_data["control_data"]["controlled_bus"] = branch_dctr[circuit]["MEASUREMENT EXTREMITY"]

            else
                sub_data["control_data"]["constraint_type"] = "setpoint"
                sub_data["control_data"]["valsp"] = nothing
            end

        elseif constraint_type == "PHASE CONTROL" # phase control
            sub_data["control_data"]["control_type"] = "shift_control"
            sub_data["control_data"]["constraint_type"] = "setpoint"
            shift_type = branch_dctr[circuit]["CONTROL TYPE"]
            sub_data["control_data"]["shift_control_variable"] = shift_type == 'C' ? "current" : shift_type == 'P' ? "power" : "fixed"
            sub_data["control_data"]["shiftmin"] = branch_dctr[circuit]["MINIMUM PHASE"] / (180/pi)
            sub_data["control_data"]["shiftmax"] = branch_dctr[circuit]["MAXIMUM PHASE"] / (180/pi)
            sub_data["control_data"]["valsp"] = branch_dctr[circuit]["SPECIFIED VALUE"] / 100

            sub_data["control_data"]["controlled_bus"] = branch_dctr[circuit]["MEASUREMENT EXTREMITY"]

        else # fix
            sub_data["control_data"]["control_type"] = "fix"
            sub_data["control_data"]["constraint_type"] = nothing
            sub_data["control_data"]["shift_control_variable"] = nothing
            sub_data["control_data"]["shiftmin"] = nothing
            sub_data["control_data"]["shiftmax"] = nothing
            sub_data["control_data"]["valsp"] = nothing
        end

        ctrl_bus = pm_data["bus"]["$(sub_data["control_data"]["controlled_bus"])"]
        sub_data["control_data"]["vmsp"] = ctrl_bus["vm"]
        sub_data["control_data"]["vmmin"] = ctrl_bus["vmin"]
        sub_data["control_data"]["vmmax"] = ctrl_bus["vmax"]
        sub_data["control_data"]["control"] =  true
        if haskey(pwf_data, "DTPF CIRC")
            for (k,v) in pwf_data["DTPF CIRC"]
                for i in 1:5
                    if v["FROM BUS $i"] == sub_data["f_bus"] && v["TO BUS $i"] == sub_data["t_bus"] && v["CIRCUIT $i"] == branch["CIRCUIT"]
                        sub_data["control_data"]["control"] = false
                    end
                end
            end
        end
    end
    idx = string(sub_data["index"])
    pm_data["branch"][idx] = sub_data
end

function _pwf2pm_transformer!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false) # Two-winding transformer
    if !haskey(pm_data, "branch")
        pm_data["branch"] = Dict{String, Any}()
    end

    if haskey(pwf_data, "DLIN")
        for (i,branch) in pwf_data["DLIN"]
            if branch["TRANSFORMER"]
                _pwf2pm_transformer!(pm_data, pwf_data, branch, add_control_data = add_control_data)
            end
        end
    end
end


text:
using Unityper, Test, Setfield, ConstructionBase

@compactify begin
    @abstract struct AT
        common_field::Int = 0
    end
    struct A <: AT
        a::Bool = true
        b::Int = 10
    end
    struct B <: AT
        a::Int = 1
        b::Float64 = 1.0
        d::Complex = 1 + 1.0im # not isbits
    end
    struct C <: AT
        b::Float64 = 2.0
        d::Bool = false
        e::Float64 = 3.0
        k::Complex{Real} = 1 + 2im # not isbits
    end
    struct D <: AT
        b::Any = "hi" # not isbits
    end
end

@test length(fieldnames(AT)) == 6
at = AT(1, 1, 1.0, true, 1 + 1im, reinterpret(var"###AT###1", Int32(0)))
# This is actually mapped to Float64! Surprise!
@test !at.a
@test at.b isa Int
@test_throws Any at.c

@test Unityper.subtypes_fun(Val(AT)) == (:A, :B, :C, :D)

a = A()
@test a.a
@test a.b === 10
a = A(; a=false, b=42)
@test !a.a
@test a.b === 42
io = IOBuffer()
show(io, MIME"text/plain"(), a)
@test String(take!(io)) == "A(a = false, b = 42)::AT"
@test Unityper.isa_type_fun(Val(AT), Val(:A), a)
@test !Unityper.isa_type_fun(Val(AT), Val(:D), a)

b = B()
@test b.a === 1
@test b.b === 1.0
@test b.d === 1 + 1.0im
b = B(a=12, b=2.0, d=2+1im)
@test b.a === 12
@test b.b === 2.0
@test b.d === 2 + 1im

c = C()
@test c.b === 2.0
@test !c.d
@test c.e === 3.0
@test c.k === 1 + 2im

c = C(b=8.0, d=true, e=5.0, k=10+10im)
@test c.b === 8.0
@test c.d
@test c.e === 5.0
@test c.k === 10 + 10im

d = D()
@test d.b == "hi"
d = D(b=100)
@test d.b === 100

@test_throws Any @macroexpand @compactified a::AT begin
    A => "A"
    B => "B"
    C => "C"
    A => "D"
end
@test_throws Any @macroexpand @compactified a::AT begin
    A => "A"
    B => "B"
    C => "C"
end
foo(a) = @compactified a::AT begin
    A => "A"
    B => "B"
    C => "C"
    D => "D"
end
@test foo(a) == "A"
@test foo(b) == "B"
@test foo(c) == "C"
@test foo(d) == "D"

@test_throws Any @macroexpand @compactified a::AT begin
    A => "A"
    _ => "B"
    C => "C"
end
goo(a) = @compactified a::AT begin
    A => "A"
    _ => "B"
end
@test goo(a) == "A"
@test goo(b) == "B"
@test goo(c) == "B"
@test goo(d) == "B"

abstract type AAT end
@compactify begin
    @abstract struct BT <: AAT
        common_field::Int = 0
    end
    struct A1 <: BT
        a::Bool = true
        b::Int = 10
    end
    struct B1 <: BT
        a::Int = 1
        b::Float64 = 1.0
        d::Complex = 1 + 1.0im # not isbits
    end
    struct C1 <: BT
        b::Float64 = 2.0
        d::Bool = false
        e::Float64 = 3.0
        k::Complex{Real} = 1 + 2im # not isbits
    end
    struct D1 <: BT
        b::Any = "hi" # not isbits
    end
end
@test BT <: AAT

@test length(fieldnames(BT)) == 6
at = BT(1, 1, 1.0, true, 1 + 1im, reinterpret(var"###BT###1", Int32(0)))
# This is actually mapped to Float64! Surprise!
@test !at.a
@test at.b isa Int
@test_throws Any at.c
@test Unityper.subtypes_fun(Val(BT)) == (:A1, :B1, :C1, :D1)

a = A1()
@test a.a
@test a.b === 10
a = A1(; a=false, b=42)
@test !a.a
@test a.b === 42
@test_throws Any Unityper.isa_type_fun(Val(AT), Val(:A1), a)
@test_throws Any Unityper.isa_type_fun(Val(BT), Val(:A), a)
@test Unityper.isa_type_fun(Val(BT), Val(:A1), a)
@test !Unityper.isa_type_fun(Val(BT), Val(:D1), a)

b = B1()
@test b.a === 1
@test b.b === 1.0
@test b.d === 1 + 1.0im
b = B1(a=12, b=2.0, d=2+1im)
@test b.a === 12
@test b.b === 2.0
@test b.d === 2 + 1im

c = C1()
@test c.b === 2.0
@test !c.d
@test c.e === 3.0
@test c.k === 1 + 2im

c = C1(b=8.0, d=true, e=5.0, k=10+10im)
@test c.b === 8.0
@test c.d
@test c.e === 5.0
@test c.k === 10 + 10im

d = D1()
@test d.b == "hi"
d = D1(b=100)
@test d.b === 100

@compactify begin
    @abstract struct AT′{T}
        common_field::Int = 0
    end
    struct A′{T} <: AT′{T}
        a::Bool = true
        b::Int = 10
    end
    struct B′{T} <: AT′{T}
        a::Int = 1
        b::Float64 = 1.0
        d::Complex = 1 + 1.0im # not isbits
    end
end

a = A′{Int}()
@compactified a::AT′ begin
    A′ => "A"
    B′ => "B"
    _ => error("unreachable")
end

@test typeof(A′{Real}()) == AT′{Real}
@test typeof(B′{Int}()) == AT′{Int}

function ConstructionBase.setproperties_object(obj::AT′{T}, patch) where T
    nt = getproperties(obj)
    nt_new = merge(nt, patch)
    Unityper.rt_constructor(obj){T}(;nt_new...)
end

@set! a.a = false
@set! a.b = 1234
@test !a.a
@test a.b === 1234


text:
using Plots, LaTeXStrings

function remap_to_domain(x::AbstractVector, xmin::Number, xmax::Number)
    X = copy(x)
    L = xmax - xmin

    for i in eachindex(X)
        while X[i] < xmin
            X[i] += L
        end
        while X[i] > xmax
            X[i] -= L
        end
    end

    return X
end

function plot_distribution(F::AbstractMatrix{T}, xGrid::AbstractVector{T}, vGrid::AbstractVector{T}) where {T}
    heatmap(xGrid, vGrid, F',
            xlims = (minimum(xGrid), maximum(xGrid)),
            ylims = (minimum(vGrid), maximum(vGrid)),
            title = "Vlasov-Poisson",
            xlabel = L"$x$",
            ylabel = L"$v$",
            legend = :none,
            grid   = :none)
end

function plot_distribution(f::Function, xGrid::AbstractVector{T}, vGrid::AbstractVector{T}) where {T}
    F = zeros(length(xGrid), length(vGrid))

    for i in eachindex(xGrid)
        for j in eachindex(vGrid)
            F[i,j] = f(xGrid[i], vGrid[j])
        end
    end

    plot_distribution(F, xGrid, vGrid)
end

function plot_distribution(f::Function, xmin::T, xmax::T, vmin::T, vmax::T; nx = 100, nv = 100) where {T}
    xGrid = collect(LinRange(xmin, xmax, nx))
    vGrid = collect(LinRange(vmin, vmax, nv))
    plot_distribution(f, xGrid, vGrid)
end

function plot_distribution(X::AbstractVector, V::AbstractVector, W::AbstractVector, xmin, xmax, vmin, vmax; nx = 100, nv = 100)
    xGrid = collect(LinRange(xmin, xmax, nx))
    vGrid = collect(LinRange(vmin, vmax, nv))

    X = remap_to_domain(X, xmin, xmax)
    F = zeros(nx, nv)

    for k in eachindex(X,V)
        if xmin < X[k] < xmax && vmin < V[k] < vmax
            i = Int(floor((X[k]-xmin)/(xmax - xmin)*nx))
            j = Int(floor((V[k]-vmin)/(vmax - vmin)*nv))
            F[i+1,j+1] += W[k]
        end
    end

    plot_distribution(F, xGrid, vGrid)
end

function plot_particles(X::AbstractVector, V::AbstractVector, W::AbstractVector, xmin, xmax, vmin, vmax)
    scatter(remap_to_domain(X, xmin, xmax), V,
            xlims = (xmin, xmax),
            ylims = (vmin, vmax),
            title = "Vlasov-Poisson",
            xlabel = L"$x$",
            ylabel = L"$v$",
            legend = :none,
            grid   = :none)
end

function plot_particles(X::AbstractVector, V::AbstractVector, W::AbstractVector, xGrid::AbstractVector, vGrid::AbstractVector)
    plot_particles(X, V, W, minimum(xGrid), maximum(xGrid), minimum(vGrid), maximum(vGrid))
end


text:
immutable Extremity <: Loss
end

@inline function loss_value(loss::Extremity, p::Real, y::Real)
    return abs(p - 0.5)
end

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


text:
using Test
using NonPromoting

# Create and examine from floating types
for T in [Float16, Float32, Float64, BigFloat]
    for x in T[1//3, 0//1, -3//2, pi]
        y1 = NP{T}(x)
        r1 = convert(T, y1)
        @test isequal(r1, x)

        y2 = NP(x)
        r2 = convert(T, y2)
        @test isequal(r2, x)
    end
end

# Create and examine from integer types
for T in [Float16, Float32, Float64, BigFloat]
    for L in [UInt8, UInt16, Int32, Int64, Int128, BigInt]
        makeL = L ∘ (L<:Unsigned ? abs : identity)
        for x in makeL.([1, 0, -3])
            y = NP{T}(x)
            r = convert(T, y)
            @test isequal(r, T(x))
        end
    end
end

# Create and examine from rational types
for T in [Float16, Float32, Float64, BigFloat]
    for L in [Rational{S} for S in [Int32, Int64, Int128, BigInt]]
        for x in L[1//3, 0//1, -3//2]
            y = NP{T}(x)
            r = convert(T, y)
            @test isequal(r, T(x))
        end
    end
end

# Ensure that operations on NP again return NP
for T in [Float16, Float32, Float64, BigFloat]
    for fun in [one, zero]
        @test typeof(fun(NP{T})) === NP{T}
    end

    for fun in [+, -, abs, sign]
        @test typeof(fun(NP{T}(1))) === NP{T}
    end

    for fun in [+, -, *, /, hypot]
        @test typeof(fun(NP{T}(1), NP{T}(2))) === NP{T}
    end
end

# Ensure binary functions don't promote
for fun in [:(-), :(/), :(\), :(^),
            :atan, :copysign, :hypot, :modf, :rem,
            :(==), :(!=), :(<), :(>), :(<=), :(>=),
            :(+), :(*), :max, :min]
    for T1 in [Float16, Float32, Float64, BigFloat]
        for T2 in [Float16, Float32, Float64, BigFloat]
            x1 = T1(1)
            x2 = T2(2)
            x3 = T1(3)
            y1 = NP(x1)
            y2 = NP(x2)
            y3 = NP(x3)

            @test isequal(x1 + x2, x3)
            if T1 === T2
                @test isequal(y1 + y2, y3)
            else
                @test_throws Union{ErrorException, MethodError} (y1 + y2)
            end
            @test_throws Union{ErrorException, MethodError} (x1 + y2)
            @test_throws Union{ErrorException, MethodError} (y1 + x2)
        end
    end
end

# Check output
for T in [Float16, Float32, Float64, BigFloat]
    for x in T[1//3, 0//1, -3//2, pi]
        buf1 = IOBuffer()
        buf2 = IOBuffer()
        show(buf1, NP(x))
        show(buf2, x)
        str1 = String(take!(buf1))
        str2 = String(take!(buf2))
        @test str1 == str2
    end
end

max_stars_repo_name: eschnett/NonPromoting

max_issues_repo_name: eschnett/NonPromoting

max_forks_repo_name: eschnett/NonPromoting


text:
# License for this file: MIT (expat)
# Copyright 2017-2020, DLR Institute of System Dynamics and Control

module Modia3D

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.5.0"
const Date = "2021-07-27"

# println("\nImporting Modia3D Version $Version ($Date)")

# Abstract types
abstract type AbstractObject3DFeature end                                 # Data associated with one Object3D
abstract type AbstractVisualElement  <: AbstractObject3DFeature     end   # Visual element associated with one Object3D
abstract type AbstractShape       <: AbstractVisualElement end   # Geometry type
abstract type AbstractGeometry    <: AbstractShape end   # Immutable shape type that  has a volume and optionally mass. Can be used in collisions.

abstract type AbstractContactMaterial   end   # Contact properties of a solid (e.g. spring constant)
abstract type AbstractContactDetection  end   # Contact detection type
abstract type AbstractGravityField      end   # Gravity field type
abstract type AbstractRenderer          end   # Renderer type
abstract type AbstractDLR_VisualizationRenderer <: AbstractRenderer end   # Community or Professional edition of DLR_Visualization renderer

abstract type AbstractMassPropertiesInterface end
abstract type AbstractMassProperties    end

abstract type AbstractContactPairMaterial end # Constants needed to compute the contact response between two objects

abstract type AbstractObject3D end
abstract type AbstractTwoObject3DObject <: AbstractObject3D end  # Object related to two Object3Ds
abstract type AbstractJoint             <: AbstractTwoObject3DObject end  # Constraint between two Object3Ds

abstract type AbstractSceneOptions end

using StaticArrays
@inline cross(x::SVector{3,Float64}, y::SVector{3,Float64}) = @inbounds SVector{3,Float64}(x[2]*y[3]-x[3]*y[2],
                                                                                x[3]*y[1]-x[1]*y[3],
                                                                                x[1]*y[2]-x[2]*y[1])

# Enumerations
@enum Ternary      True False Inherited

"""
    @enum AnalysisType KinematicAnalysis QuasiStaticAnalysis DynamicAnalysis

Type of analyis that is actually carried out. The `AnalysisType` is set by the user
of the simulation model. Variables are declared in the model with [`VariableAnalysisType`](@ref)`.
Variables with [`VariableAnalysisType`](@ref)` <= AnalysisType` are removed from the analysis and do
not show up in the result. For example, an *acceleration* would be declared as `OnlyDynamicAnalysis`
and then this variable would not show up in the result, if `AnalysisType = KinematicAnalysis` or
`QuasiStaticAnalysis`.

Currently, only DynamicAnalysis is supported and used.
"""
@enum AnalysisType KinematicAnalysis QuasiStaticAnalysis DynamicAnalysis
@enum VariableAnalysisType AllAnalysis QuasiStaticAndDynamicAnalysis OnlyDynamicAnalysis NotUsedInAnalysis

# Used renderer (actual value is defined with __init__() below)
const renderer = Vector{AbstractRenderer}(undef,1)

import Unitful

numberType(value) = ModiaLang.baseType(eltype(value))

"""
    convertAndStripUnit(TargetType, requiredUnit, value)

Return the `value` of a variable converted to the `requiredUnit` and the required `TargetType`.
The `value` can be a scalar, array or collection of a primitive type, such as Float64,
or of a primitive type with units, or of a Measurement type.

- If `value` has no unit, it is converted to `TargetType`.
- If `value` has a unit, it is converted to the `requiredUnit`, then the unit is stripped,
  and finally it is converted to `TargetType`.

# Example

```julia
using Unitful

# :L has unit u"m"
convertAndStripUnit(Float32, 0.01u"km", u"m")  # = 10.0f0
convertAndStripUnit(Float32, 10.0)             # = 10.0f0
```
"""
convertAndStripUnit(TargetType, requiredUnit, value) =
    numberType(value) <: Unitful.AbstractQuantity && unit.(value) != Unitful.NoUnits ?
            convert(TargetType, ustrip.( uconvert.(requiredUnit, value))) : convert(TargetType, value)

# Include sub-modules
include(joinpath("Frames"           , "_module.jl"))
include(joinpath("Basics"           , "_module.jl"))
include(joinpath("Shapes"         , "_module.jl"))
include(joinpath("Composition"      , "_module.jl"))
include(joinpath("AnimationExport"  , "_module.jl"))
include(joinpath("PathPlanning"     , "_module.jl"))
include(joinpath("renderer"         , "DLR_Visualization"  , "_module.jl"))
include(joinpath("renderer"         , "NoRenderer"         , "_module.jl"))
include(joinpath("contactDetection" , "ContactDetectionMPR", "_module.jl"))
include(joinpath("Interface"        , "_module.jl"))
include(joinpath("ModiaInterface"   , "_module.jl"))

# Make symbols available that have been exported in sub-modules
using  .Frames
using  .Basics
using  .Shapes
using  .AnimationExport
using  .Composition
import .DLR_Visualization
import .NoRenderer
using .PathPlanning
using .Interface
#const connect = Composition.connect  # connect cannot be directly exported, due to a conflict with Base.connect
const run     = Interface.run        # run cannot be directly exported, due to a conflict with Base.run

# Called implicitely at the first import/using of Modia3D (when loading Modia3D to the current Julia session)
function __init__()
    info = DLR_Visualization.getSimVisInfo()
    (directory, dll_name, isProfessionalEdition, isNoRenderer) = info

    if isNoRenderer
        renderer[1] = NoRenderer.DummyRenderer(info)
    elseif isProfessionalEdition
        renderer[1] = DLR_Visualization.ProfessionalEdition(info)
    else
        renderer[1] = DLR_Visualization.CommunityEdition(info)
end; end

export Object3D

export Sphere, Ellipsoid, Box, Cylinder, Capsule, Beam, Cone
export Spring, GearWheel, CoordinateSystem, Grid, FileMesh

export Solid, Visual
export MassProperties
export Fix
export Revolute, Prismatic

export UniformGravityField, PointGravityField, NoGravityField
export VisualMaterial
export SceneOptions, SimulationModel
export print_ModelVariables

export PTP_path, pathEndTime, getPosition!, getPosition, getIndex, plotPath

export calculateRobotMovement

# Add import clauses used in examples and test
import StaticArrays
import Unitful
import LinearAlgebra
import Test
import ModiaLang

end # module


text:
function max_ones(bitstream ::String)
    max_count = 0
    current = 0
    for i = 1:length(bitstream)
        if bitstream[i] == '1'
            current += 1
        else
            if current > max_count
                max_count = current
            end
            current = 0
        end
    end
    return max(max_count,current)
end

n = parse(Int64, readline())
# Get binary interpretation of n
binary_n = bitstring(n)
# println(binary_n)
println(max_ones(binary_n))

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 10.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
using Distributions

const T_1df = [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0]
const T_2df = [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 2.0]
const xv    = (-10.0, -1.0, -0.1, 0.0, 0.1, 1.0, 10.0)

@testset "PDF" begin
  for x in xv
    @test isapprox(pdf(BulkAndTailsDist(T_1df), x), pdf(TDist(1.0), x))
    @test isapprox(pdf(BulkAndTailsDist(T_2df), x), pdf(TDist(2.0), x))
  end
end

@testset "LOGPDF" begin
  for x in xv
    @test isapprox(logpdf(BulkAndTailsDist(T_1df), x), logpdf(TDist(1.0), x))
    @test isapprox(logpdf(BulkAndTailsDist(T_2df), x), logpdf(TDist(2.0), x))
  end
end

@testset "CDF" begin
  for x in xv
    @test isapprox(cdf(BulkAndTailsDist(T_1df), x), cdf(TDist(1.0), x))
    @test isapprox(cdf(BulkAndTailsDist(T_2df), x), cdf(TDist(2.0), x))
  end
end


text:
function duftau!(
    champu :: Array{ComplexF64,2},
    champv :: Array{ComplexF64,2},
    m      :: MicMac,
    t      :: Float64,
    fft_u  :: Vector{ComplexF64},
    fft_v  :: Vector{ComplexF64},
    fft_du :: Vector{ComplexF64},
    fft_dv :: Vector{ComplexF64} )

    sigma = 1

    m.u .= exp.(1im .* t .* m.A1)
    m.v .= -1im .* m.llambda .* m.A2 .* conj.(m.u)

    m.ut .= (fft_u .* m.u) .* m.matr
    ifft!(m.ut, 1)

    m.vt .= (fft_v .* m.u) .* m.matr
    ifft!(m.vt, 1)

    m.z  .= (m.ut  .+ conj.(m.vt)) ./ 2

    m.ut .= (m.u .* fft_du) .* m.matr
    ifft!(m.ut, 1)

    m.vt .= (m.u .* fft_dv) .* m.matr
    ifft!(m.vt, 1)

    m.ut .= (m.ut .+ conj.(m.vt)) ./ 2

    m.vt .= 2 .* abs.(m.z).^2 .* m.ut .+ m.z.^2 .* conj.(m.ut)

    m.ut .= m.conjmatr .* m.vt
    m.vt .= m.conjmatr .* conj.(m.vt)

    fft!(m.ut,1)
    fft!(m.vt,1)

    m.ut .*= m.v
    m.vt .*= m.v

    transpose!(champu, m.ut)
    transpose!(champv, m.vt)

end


text:
using Base.Test
using AutomotiveDrivingModels
using ConditionalLinearGaussianDrivers


text:
module Inference

using ..Core
using ..Core: logZ
using ..Utilities
using DynamicPPL: Metadata, _tail, VarInfo, TypedVarInfo,
    islinked, invlink!, getlogp, tonamedtuple, VarName, getsym, vectorize,
    settrans!, _getvns, getdist, CACHERESET, AbstractSampler,
    Model, Sampler, SampleFromPrior, SampleFromUniform,
    Selector, AbstractSamplerState, DefaultContext, PriorContext,
    LikelihoodContext, MiniBatchContext, set_flag!, unset_flag!, NamedDist, NoDist,
    getspace, inspace
using Distributions, Libtask, Bijectors
using DistributionsAD: VectorOfMultivariate
using LinearAlgebra
using ..Turing: PROGRESS, Turing
using StatsFuns: logsumexp
using Random: GLOBAL_RNG, AbstractRNG, randexp
using DynamicPPL
using AbstractMCMC: AbstractModel, AbstractSampler
using DocStringExtensions: TYPEDEF, TYPEDFIELDS

import AbstractMCMC
import AdvancedHMC; const AHMC = AdvancedHMC
import AdvancedMH; const AMH = AdvancedMH
import ..Core: getchunksize, getADbackend
import DynamicPPL: get_matching_type,
    VarName, _getranges, _getindex, getval, _getvns
import EllipticalSliceSampling
import Random
import MCMCChains
import StatsBase: predict

export  InferenceAlgorithm,
        Hamiltonian,
        GibbsComponent,
        StaticHamiltonian,
        AdaptiveHamiltonian,
        SampleFromUniform,
        SampleFromPrior,
        MH,
        ESS,
        Emcee,
        Gibbs,      # classic sampling
        HMC,
        SGLD,
        SGHMC,
        HMCDA,
        NUTS,       # Hamiltonian-like sampling
        DynamicNUTS,
        IS,
        SMC,
        CSMC,
        PG,
        Prior,
        assume,
        dot_assume,
        observe,
        dot_observe,
        resume,
        predict,
        isgibbscomponent

#######################
# Sampler abstraction #
#######################
abstract type AbstractAdapter end
abstract type InferenceAlgorithm end
abstract type ParticleInference <: InferenceAlgorithm end
abstract type Hamiltonian{AD} <: InferenceAlgorithm end
abstract type StaticHamiltonian{AD} <: Hamiltonian{AD} end
abstract type AdaptiveHamiltonian{AD} <: Hamiltonian{AD} end

getchunksize(::Type{<:Hamiltonian{AD}}) where AD = getchunksize(AD)
getADbackend(::Hamiltonian{AD}) where AD = AD()

# Algorithm for sampling from the prior
struct Prior <: InferenceAlgorithm end

"""
    mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)

Decide if a proposal ``x'`` with log probability ``\\log p(x') = logp_proposal`` and
log proposal ratio ``\\log k(x', x) - \\log k(x, x') = log_proposal_ratio`` in a
Metropolis-Hastings algorithm with Markov kernel ``k(x_t, x_{t+1})`` and current state
``x`` with log probability ``\\log p(x) = logp_current`` is accepted by evaluating the
Metropolis-Hastings acceptance criterion
```math
\\log U \\leq \\log p(x') - \\log p(x) + \\log k(x', x) - \\log k(x, x')
```
for a uniform random number ``U \\in [0, 1)``.
"""
function mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)
    # replacing log(rand()) with -randexp() yields test errors
    return log(rand()) + logp_current ≤ logp_proposal + log_proposal_ratio
end

######################
# Default Transition #
######################

struct Transition{T, F<:AbstractFloat}
    θ  :: T
    lp :: F
end

function Transition(spl::Sampler, nt::NamedTuple=NamedTuple())
    theta = merge(tonamedtuple(spl.state.vi), nt)
    lp = getlogp(spl.state.vi)
    return Transition{typeof(theta), typeof(lp)}(theta, lp)
end

function additional_parameters(::Type{<:Transition})
    return [:lp]
end

DynamicPPL.getlogp(t::Transition) = t.lp

##########################################
# Internal variable names for MCMCChains #
##########################################

const TURING_INTERNAL_VARS = (internals = [
    "elapsed",
    "eval_num",
    "lf_eps",
    "lp",
    "weight",
    "le",
    "acceptance_rate",
    "hamiltonian_energy",
    "hamiltonian_energy_error",
    "max_hamiltonian_energy_error",
    "is_accept",
    "log_density",
    "n_steps",
    "numerical_error",
    "step_size",
    "nom_step_size",
    "tree_depth",
    "is_adapt"
],)

#########################################
# Default definitions for the interface #
#########################################

function AbstractMCMC.sample(
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer;
    kwargs...
)
    return AbstractMCMC.sample(Random.GLOBAL_RNG, model, alg, N; kwargs...)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::InferenceAlgorithm,
    N::Integer;
    kwargs...
)
    return AbstractMCMC.sample(rng, model, Sampler(alg, model), N; kwargs...)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    sampler::Sampler{<:InferenceAlgorithm},
    N::Integer;
    chain_type=MCMCChains.Chains,
    resume_from=nothing,
    progress=PROGRESS[],
    kwargs...
)
    if resume_from === nothing
        return AbstractMCMC.mcmcsample(rng, model, sampler, N;
                                       chain_type=chain_type, progress=progress, kwargs...)
    else
        return resume(resume_from, N; chain_type=chain_type, progress=progress, kwargs...)
    end
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::Prior,
    N::Integer;
    chain_type=MCMCChains.Chains,
    resume_from=nothing,
    progress=PROGRESS[],
    kwargs...
)
    if resume_from === nothing
        return AbstractMCMC.mcmcsample(rng, model, SampleFromPrior(), N;
                                       chain_type=chain_type, progress=progress, kwargs...)
    else
        return resume(resume_from, N; chain_type=chain_type, progress=progress, kwargs...)
    end
end

function AbstractMCMC.sample(
    model::AbstractModel,
    alg::InferenceAlgorithm,
    parallel::AbstractMCMC.AbstractMCMCParallel,
    N::Integer,
    n_chains::Integer;
    kwargs...
)
    return AbstractMCMC.sample(Random.GLOBAL_RNG, model, alg, parallel, N, n_chains;
                               kwargs...)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::InferenceAlgorithm,
    parallel::AbstractMCMC.AbstractMCMCParallel,
    N::Integer,
    n_chains::Integer;
    kwargs...
)
    return AbstractMCMC.sample(rng, model, Sampler(alg, model), parallel, N, n_chains;
                               kwargs...)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    sampler::Sampler{<:InferenceAlgorithm},
    parallel::AbstractMCMC.AbstractMCMCParallel,
    N::Integer,
    n_chains::Integer;
    chain_type=MCMCChains.Chains,
    progress=PROGRESS[],
    kwargs...
)
    return AbstractMCMC.mcmcsample(rng, model, sampler, parallel, N, n_chains;
                                   chain_type=chain_type, progress=progress, kwargs...)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    alg::Prior,
    parallel::AbstractMCMC.AbstractMCMCParallel,
    N::Integer,
    n_chains::Integer;
    chain_type=MCMCChains.Chains,
    progress=PROGRESS[],
    kwargs...
)
    return AbstractMCMC.sample(rng, model, SampleFromPrior(), parallel, N, n_chains;
                               chain_type=chain_type, progress=progress, kwargs...)
end

function AbstractMCMC.sample_init!(
    ::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:InferenceAlgorithm},
    N::Integer;
    kwargs...
)
    # Resume the sampler.
    set_resume!(spl; kwargs...)

    # Set the parameters to a starting value.
    initialize_parameters!(spl; kwargs...)
end

function initialize_parameters!(
    spl::Sampler;
    init_theta::Union{Nothing,Vector}=nothing,
    verbose::Bool=false,
    kwargs...
)
    islinked(spl.state.vi, spl) && invlink!(spl.state.vi, spl)
    # Get `init_theta`
    if init_theta !== nothing
        verbose && @info "Using passed-in initial variable values" init_theta
        # Convert individual numbers to length 1 vector; `ismissing(v)` is needed as `size(missing)` is undefined`
        init_theta = [ismissing(v) || size(v) == () ? [v] : v for v in init_theta]
        # Flatten `init_theta`
        init_theta_flat = foldl(vcat, map(vec, init_theta))
        # Create a mask to inidicate which values are not missing
        theta_mask = map(x -> !ismissing(x), init_theta_flat)
        # Get all values
        theta = spl.state.vi[spl]
        @assert length(theta) == length(init_theta_flat) "Provided initial value doesn't match the dimension of the model"
        # Update those which are provided (i.e. not missing)
        theta[theta_mask] .= init_theta_flat[theta_mask]
        # Update in `vi`
        spl.state.vi[spl] = theta
    end
end

##########################
# Chain making utilities #
##########################

"""
    getparams(t)

Return a named tuple of parameters.
"""
getparams(t) = t.θ
getparams(t::VarInfo) = tonamedtuple(TypedVarInfo(t))

function _params_to_array(ts::Vector)
    names = Vector{Symbol}()
    # Extract the parameter names and values from each transition.
    dicts = map(ts) do t
        nms, vs = flatten_namedtuple(getparams(t))
        for nm in nms
            if !(nm in names)
                push!(names, nm)
            end
        end
        # Convert the names and values to a single dictionary.
        return Dict(nms[j] => vs[j] for j in 1:length(vs))
    end
    # names = collect(names_set)
    vals = [get(dicts[i], key, missing) for i in eachindex(dicts),
        (j, key) in enumerate(names)]

    return names, vals
end

function flatten_namedtuple(nt::NamedTuple)
    names_vals = mapreduce(vcat, keys(nt)) do k
        v = nt[k]
        if length(v) == 1
            return [(Symbol(k), v)]
        else
            return mapreduce(vcat, zip(v[1], v[2])) do (vnval, vn)
                return collect(FlattenIterator(vn, vnval))
            end
        end
    end
    return [vn[1] for vn in names_vals], [vn[2] for vn in names_vals]
end

function get_transition_extras(ts::AbstractVector{<:VarInfo})
    valmat = reshape([getlogp(t) for t in ts], :, 1)
    return [:lp], valmat
end

function get_transition_extras(ts::AbstractVector)
    # Get the extra field names from the sampler state type.
    # This handles things like :lp or :weight.
    extra_params = additional_parameters(eltype(ts))

    # Get the values of the extra parameters.
    local extra_names
    all_vals = []

    # Iterate through each transition.
    for t in ts
        extra_names = Symbol[]
        vals = []

        # Iterate through each of the additional field names
        # in the struct.
        for p in extra_params
            # Check whether the field contains a NamedTuple,
            # in which case we need to iterate through each
            # key/value pair.
            prop = getproperty(t, p)
            if prop isa NamedTuple
                for (k, v) in pairs(prop)
                    push!(extra_names, Symbol(k))
                    push!(vals, v)
                end
            else
                push!(extra_names, Symbol(p))
                push!(vals, prop)
            end
        end
        push!(all_vals, vals)
    end

    # Convert the vector-of-vectors to a matrix.
    valmat = [all_vals[i][j] for i in 1:length(ts), j in 1:length(all_vals[1])]

    return extra_names, valmat
end

getlogevidence(sampler) = missing
function getlogevidence(sampler::Sampler)
    if isdefined(sampler.state, :average_logevidence)
        return sampler.state.average_logevidence
    elseif isdefined(sampler.state, :final_logevidence)
        return sampler.state.final_logevidence
    else
        return missing
    end
end

# Default MCMCChains.Chains constructor.
# This is type piracy (at least for SampleFromPrior).
function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Union{Sampler{<:InferenceAlgorithm},SampleFromPrior},
    N::Integer,
    ts::Vector,
    chain_type::Type{MCMCChains.Chains};
    save_state = false,
    kwargs...
)
    # Convert transitions to array format.
    # Also retrieve the variable names.
    nms, vals = _params_to_array(ts)

    # Get the values of the extra parameters in each transition.
    extra_params, extra_values = get_transition_extras(ts)

    # Extract names & construct param array.
    nms = [nms; extra_params]
    parray = hcat(vals, extra_values)

    # Get the average or final log evidence, if it exists.
    le = getlogevidence(spl)

    # Set up the info tuple.
    if save_state
        info = (range = rng, model = model, spl = spl)
    else
        info = NamedTuple()
    end

    # Conretize the array before giving it to MCMCChains.
    parray = MCMCChains.concretize(parray)

    # Chain construction.
    return MCMCChains.Chains(
        parray,
        nms,
        deepcopy(TURING_INTERNAL_VARS);
        evidence=le,
        info=info,
    ) |> sort
end

# This is type piracy (for SampleFromPrior).
function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Union{Sampler{<:InferenceAlgorithm},SampleFromPrior},
    N::Integer,
    ts::Vector,
    chain_type::Type{Vector{NamedTuple}};
    discard_adapt::Bool=true,
    save_state=false,
    kwargs...
)
    nts = Vector{NamedTuple}(undef, N)

    for (i, t) in enumerate(ts)
        params = getparams(t)

        k = collect(keys(params))
        vs = []
        for v in values(params)
            push!(vs, v[1])
        end

        push!(k, :lp)

        nts[i] = NamedTuple{tuple(k...)}(tuple(vs..., getlogp(t)))
    end

    return map(identity, nts)
end

function save(c::MCMCChains.Chains, spl::Sampler, model, vi, samples)
    nt = NamedTuple{(:spl, :model, :vi, :samples)}((spl, model, deepcopy(vi), samples))
    return setinfo(c, merge(nt, c.info))
end

function resume(
    c::MCMCChains.Chains,
    n_iter::Int;
    chain_type=MCMCChains.Chains,
    progress=PROGRESS[],
    kwargs...
)
    @assert !isempty(c.info) "[Turing] cannot resume from a chain without state info"

    # Sample a new chain.
    newchain = AbstractMCMC.mcmcsample(
        c.info[:range],
        c.info[:model],
        c.info[:spl],
        n_iter;
        resume_from=c,
        reuse_spl_n=n_iter,
        chain_type=MCMCChains.Chains,
        progress=progress,
        kwargs...
    )

    # Stick the new samples at the end of the old chain.
    return vcat(c, newchain)
end

function set_resume!(
    s::Sampler;
    resume_from::Union{MCMCChains.Chains, Nothing}=nothing,
    kwargs...
)
    # If we're resuming, grab the sampler info.
    if resume_from !== nothing
        s = resume_from.info[:spl]
    end
end

#########################
# Default sampler state #
#########################

"""
A blank `AbstractSamplerState` that contains only `VarInfo` information.
"""
mutable struct SamplerState{VIType<:VarInfo} <: AbstractSamplerState
    vi :: VIType
end

#######################################
# Concrete algorithm implementations. #
#######################################

include("ess.jl")
include("hmc.jl")
include("mh.jl")
include("is.jl")
include("AdvancedSMC.jl")
include("gibbs.jl")
include("../contrib/inference/sghmc.jl")
include("emcee.jl")

################
# Typing tools #
################

for alg in (:SMC, :PG, :MH, :IS, :ESS, :Gibbs, :Emcee)
    @eval DynamicPPL.getspace(::$alg{space}) where {space} = space
end
for alg in (:HMC, :HMCDA, :NUTS, :SGLD, :SGHMC)
    @eval DynamicPPL.getspace(::$alg{<:Any, space}) where {space} = space
end

floatof(::Type{T}) where {T <: Real} = typeof(one(T)/one(T))
floatof(::Type) = Real # fallback if type inference failed

function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{T},
) where {T}
    return T
end
function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{<:Union{Missing, AbstractFloat}},
)
    return Union{Missing, floatof(eltype(vi, spl))}
end
function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T,N}}) where {T,N}
    return Array{get_matching_type(spl, vi, T), N}
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T}}) where T
    return Array{get_matching_type(spl, vi, T)}
end
function get_matching_type(
    spl::Sampler{<:Union{PG, SMC}},
    vi,
    ::Type{TV},
) where {T, N, TV <: Array{T, N}}
    return TArray{T, N}
end

##############
# Utilities  #
##############

DynamicPPL.getspace(spl::Sampler) = getspace(spl.alg)
DynamicPPL.inspace(vn::VarName, spl::Sampler) = inspace(vn, getspace(spl.alg))

"""

    predict([rng::AbstractRNG,] model::Model, chain::MCMCChains.Chains; include_all=false)

Execute `model` conditioned on each sample in `chain`, and return the resulting `Chains`.

If `include_all` is `false`, the returned `Chains` will contain only those variables
sampled/not present in `chain`.

# Details
Internally calls `Turing.Inference.transitions_from_chain` to obtained the samples
and then converts these into a `Chains` object using `AbstractMCMC.bundle_samples`.

# Example
```jldoctest
julia> using Turing; setprogress!(false);
[ Info: [Turing]: progress logging is disabled globally

julia> @model function linear_reg(x, y, σ = 0.1)
           β ~ Normal(0, 1)

           for i ∈ eachindex(y)
               y[i] ~ Normal(β * x[i], σ)
           end
       end;

julia> σ = 0.1; f(x) = 2 * x + 0.1 * randn();

julia> Δ = 0.1; xs_train = 0:Δ:10; ys_train = f.(xs_train);

julia> xs_test = [10 + Δ, 10 + 2 * Δ]; ys_test = f.(xs_test);

julia> m_train = linear_reg(xs_train, ys_train, σ);

julia> chain_lin_reg = sample(m_train, NUTS(100, 0.65), 200);
┌ Info: Found initial step size
└   ϵ = 0.003125

julia> m_test = linear_reg(xs_test, Vector{Union{Missing, Float64}}(undef, length(ys_test)), σ);

julia> predictions = predict(m_test, chain_lin_reg)
Object of type Chains, with data of type 100×2×1 Array{Float64,3}

Iterations        = 1:100
Thinning interval = 1
Chains            = 1
Samples per chain = 100
parameters        = y[1], y[2]

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters     mean     std  naive_se     mcse       ess   r_hat
  ──────────  ───────  ──────  ────────  ───────  ────────  ──────
        y[1]  20.1974  0.1007    0.0101  missing  101.0711  0.9922
        y[2]  20.3867  0.1062    0.0106  missing  101.4889  0.9903

Quantiles
  parameters     2.5%    25.0%    50.0%    75.0%    97.5%
  ──────────  ───────  ───────  ───────  ───────  ───────
        y[1]  20.0342  20.1188  20.2135  20.2588  20.4188
        y[2]  20.1870  20.3178  20.3839  20.4466  20.5895

julia> ys_pred = vec(mean(Array(group(predictions, :y)); dims = 1));

julia> sum(abs2, ys_test - ys_pred) ≤ 0.1
true
```
"""
function predict(model::Model, chain::MCMCChains.Chains; kwargs...)
    return predict(Random.GLOBAL_RNG, model, chain; kwargs...)
end
function predict(rng::AbstractRNG, model::Model, chain::MCMCChains.Chains; include_all = false)
    spl = DynamicPPL.SampleFromPrior()

    # Sample transitions using `spl` conditioned on values in `chain`
    transitions = [
        transitions_from_chain(rng, model, chain[:, :, chn_idx]; sampler = spl)
        for chn_idx = 1:size(chain, 3)
    ]

    # Let the Turing internals handle everything else for you
    chain_result = reduce(
        MCMCChains.chainscat, [
            AbstractMCMC.bundle_samples(
                rng,
                model,
                spl,
                length(chain),
                transitions[chn_idx],
                MCMCChains.Chains
            ) for chn_idx = 1:size(chain, 3)
        ]
    )

    parameter_names = if include_all
        names(chain_result, :parameters)
    else
        filter(k -> ∉(k, names(chain, :parameters)), names(chain_result, :parameters))
    end

    return chain_result[parameter_names]
end

"""

    transitions_from_chain(
        [rng::AbstractRNG,]
        model::Model,
        chain::MCMCChains.Chains;
        sampler = DynamicPPL.SampleFromPrior()
    )

Execute `model` conditioned on each sample in `chain`, and return resulting transitions.

The returned transitions are represented in a `Vector{<:Turing.Inference.Transition}`.

# Details

In a bit more detail, the process is as follows:
1. For every `sample` in `chain`
   1. For every `variable` in `sample`
      1. Set `variable` in `model` to its value in `sample`
   2. Execute `model` with variables fixed as above, sampling variables NOT present
      in `chain` using `SampleFromPrior`
   3. Return sampled variables and log-joint

# Example
```julia-repl
julia> using Turing

julia> @model function demo()
           m ~ Normal(0, 1)
           x ~ Normal(m, 1)
       end;

julia> m = demo();

julia> chain = Chains(randn(2, 1, 1), ["m"]); # 2 samples of `m`

julia> transitions = Turing.Inference.transitions_from_chain(m, chain);

julia> [Turing.Inference.getlogp(t) for t in transitions] # extract the logjoints
2-element Array{Float64,1}:
 -3.6294991938628374
 -2.5697948166987845

julia> [first(t.θ.x) for t in transitions] # extract samples for `x`
2-element Array{Array{Float64,1},1}:
 [-2.0844148956440796]
 [-1.704630494695469]
```
"""
function transitions_from_chain(
    model::Turing.Model,
    chain::MCMCChains.Chains;
    kwargs...
)
    return transitions_from_chain(Random.GLOBAL_RNG, model, chain; kwargs...)
end
function transitions_from_chain(
    rng::AbstractRNG,
    model::Turing.Model,
    chain::MCMCChains.Chains;
    sampler = DynamicPPL.SampleFromPrior()
)
    vi = Turing.VarInfo(model)

    transitions = map(1:length(chain)) do i
        c = chain[i]
        md = vi.metadata
        for v in keys(md)
            for vn in md[v].vns
                vn_sym = Symbol(vn)

                # Cannot use `vn_sym` to index in the chain
                # so we have to extract the corresponding "linear"
                # indices and use those.
                # `ks` is empty if `vn_sym` not in `c`.
                ks = MCMCChains.namesingroup(c, vn_sym)

                if !isempty(ks)
                    # 1st dimension is of size 1 since `c`
                    # only contains a single sample, and the
                    # last dimension is of size 1 since
                    # we're assuming we're working with a single chain.
                    val = copy(vec(c[ks].value))
                    DynamicPPL.setval!(vi, val, vn)
                    DynamicPPL.settrans!(vi, false, vn)
                else
                    DynamicPPL.set_flag!(vi, vn, "del")
                end
            end
        end
        # Execute `model` on the parameters set in `vi` and sample those with `"del"` flag using `sampler`
        model(rng, vi, sampler)

        # Convert `VarInfo` into `NamedTuple` and save
        theta = DynamicPPL.tonamedtuple(vi)
        lp = Turing.getlogp(vi)
        Transition(theta, lp)
    end

    return transitions
end

end # module


text:
struct YarnManager <: ClusterManager
    ugi::UserGroupInformation
    clnt::YarnClient
    am::YarnAppMaster
    app::YarnApp
    stdout_ios::Vector{IO}
    launch_timeout::Integer
    keep_connected::Bool

    function YarnManager(; kwargs...)
        params = Dict()
        for (n,v) in kwargs
            params[n] = v
        end
        paramkeys = keys(params)
        @debug("YarnManager constructor: params: $params")

        user            = (:user            in paramkeys) ? params[:user]           : ""
        rmport          = (:rmport          in paramkeys) ? params[:rmport]         : 8032
        yarnhost        = (:yarnhost        in paramkeys) ? params[:yarnhost]       : "localhost"
        schedport       = (:schedport       in paramkeys) ? params[:schedport]      : 8030
        launch_timeout  = (:launch_timeout  in paramkeys) ? params[:launch_timeout] : 60
        keep_connected  = (:keep_connected  in paramkeys) ? params[:keep_connected] : true

        ugi     = UserGroupInformation(user)
        clnt    = YarnClient(yarnhost, rmport, ugi)
        am      = YarnAppMaster(yarnhost, schedport, ugi)
        app     = submit(clnt, am)

        new(ugi, clnt, am, app, IO[], launch_timeout, keep_connected)
    end
end

function disconnect(yarncm::YarnManager)
    ret = unregister(yarncm.am, true)
end

function show(io::IO, yarncm::YarnManager)
    print(io, "YarnManager for ")
    show(io, yarncm.clnt)
end

function setup_worker(host, port, cookie)
    @debug("YarnManager setup_worker: host:$host port:$port cookie:$cookie for container $(ENV[CONTAINER_ID])")
    c = connect(IPv4(host), port)
    Base.wait_connected(c)
    redirect_stdout(c)
    redirect_stderr(c)
    # identify container id so that rmprocs can clean things up nicely if required
    serialize(c, ENV[CONTAINER_ID])
    if cookie == nothing
        start_worker(c)
    else
        if isa(cookie, Symbol)
            cookie = string(cookie)[8:end] # strip the leading "cookie_"
        end
        start_worker(c, cookie)
    end
end

# used to reconstruct container id object from the environment string
function _container(cid::String)
    parts = split(cid, '_')
    (parts[1] == "container") || throw(YarnException("Invalid container id $cid"))
    cluster_timestamp = parse(Int64, parts[2])
    app_id = parse(Int32, parts[3])
    attempt_id = parse(Int32, parts[4])
    container_id = parse(Int64, parts[5])

    appid_proto = protobuild(ApplicationIdProto, Dict(:id => app_id, :cluster_timestamp => cluster_timestamp))
    app_attempt_id_proto = protobuild(ApplicationAttemptIdProto, Dict(:application_id => appid_proto, :attemptId => attempt_id))
    protobuild(ContainerIdProto, Dict(:app_id => appid_proto, :app_attempt_id => app_attempt_id_proto, :id => container_id))
end

_envdict(envdict::Dict) = envdict
function _envdict(envhash::Base.EnvDict)
    envdict = Dict{String,String}()
    for (n,v) in envhash
        envdict[n] = v
    end
    envdict
end

function _currprocname()
    p = joinpath(Sys.BINDIR, Sys.get_process_title())
    exists(p) && (return p)

    ("_" in keys(ENV)) && contains(ENV["_"], "julia") && (return ENV["_"])

    "julia"
end

function launch(manager::YarnManager, params::Dict, instances_arr::Array, c::Condition)
    @debug("YarnManager launch: params: $params")

    paramkeys   = keys(params)
    np          = (:np       in paramkeys)  ? params[:np]               : 1
    cmd         = (:exename  in paramkeys)  ? params[:exename]          : _currprocname()
    appenv      = (:env      in paramkeys)  ? _envdict(params[:env])    : Dict{String,String}()
    mem         = (:mem      in paramkeys)  ? params[:mem]              : YARN_CONTAINER_MEM_DEFAULT
    cpu         = (:cpu      in paramkeys)  ? params[:cpu]              : YARN_CONTAINER_CPU_DEFAULT
    loc         = (:loc      in paramkeys)  ? params[:loc]              : YARN_CONTAINER_LOCATION_DEFAULT
    priority    = (:priority in paramkeys)  ? params[:priority]         : YARN_CONTAINER_PRIORITY_DEFAULT

    stdout_ios = manager.stdout_ios
    ipaddr = getipaddr()
    (port, server) = listenany(ipaddr, 11000)
    rcv_stdouts = @async begin
        while length(stdout_ios) < np
            sock = accept(server)
            push!(stdout_ios, sock)
        end
    end

    cookie = ":cookie_" * cluster_cookie()
    initargs = "using Elly; Elly.setup_worker($(ipaddr.host), $(port), $(cookie))"
    clc = launchcontext(cmd="$cmd -e '$initargs'", env=appenv)

    @debug("YarnManager launch: initargs: $initargs")
    @debug("YarnManager launch: context: $clc")
    on_alloc = (cid) -> container_start(manager.am, cid, clc)
    callback(manager.am, on_alloc, nothing)

    container_allocate(manager.am, convert(Int, np); mem=mem, cpu=cpu, loc=loc, priority=priority)

    tend = time() + manager.launch_timeout
    while (time() < tend) && !istaskdone(rcv_stdouts)
        sleep(1)
    end

    (length(stdout_ios) == np) || throw(YarnException("Error starting all workers. $(length(stdout_ios)) of $(np) started."))

    configs = WorkerConfig[]
    @sync begin
        for io in manager.stdout_ios
            @async let io=io
                config = WorkerConfig()
                config.io = io

                # keep the container id in userdata to be used with rmprocs if required
                cidstr = deserialize(io)
                userdata = Dict(:container_id => _container(cidstr))
                config.userdata = userdata

                push!(configs, config)
            end
        end
        empty!(manager.stdout_ios)
    end

    append!(instances_arr, configs)
    notify(c)
    nothing
end

function manage(manager::YarnManager, id::Integer, config::WorkerConfig, op::Symbol)
    # This function needs to exist, but so far we don't do anything
    if op == :deregister
        @debug("YarnManager manage: id:$id, op:$op, nprocs:$(nprocs())")
        !manager.keep_connected && (1 == nprocs()) && (@async disconnect(manager))
    end
    nothing
end


text:
using HiQGA.transD_GP, PyPlot, Random
## model fixed parts, i.e., air, but only to be compatible with AEM ...
zfixed   = [-1e5]
ρfixed   = [1e12]
nmax = 200
# z grid spec starts, the first z and first ρ will be unused in MT
zstart = 0.
extendfrac1, dz1, n1 = 1.06, 1.15, 50
# top grid for AEM
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac1, dz=dz1, n=n1, showplot=false, forextension=true)
# deeper grid for MT
extendfrac2, dz2, n2 = 1.5, nothing, 15
zall, znall, zboundaries = transD_GP.zcontinue(;zall=zall, znall=znall, zboundaries=zboundaries,
                                            extendfrac=extendfrac2, dz=dz2, n=n2, showplot=true)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<20)] .= 800.
ρ[(z.>=20) .&(z.<40)] .= 100
ρ[(z.>=40) .&(z.<50)] .= 10
ρ[(z.>=50) .&(z.<70)] .= 1
ρ[(z.>=70) .&(z.<100)] .= 20
ρ[(z.>=100) .&(z.<120)] .= 80
ρ[(z.>=120) .&(z.<140)] .= 2
ρ[(z.>=140) .&(z.<250)] .= 100.
ρ[(z.>=250) .&(z.<270)] .= 100.
ρ[(z.>=270)]            .= 150
# add jitter to model in log10 domain
Random.seed!(11)
ρ[1:50] = 10 .^(0.1*randn(length(ρ[1:50])) + log10.(ρ[1:50]))
##
T = 10 .^range(-3, 0, length=10)[1:5]
F = transD_GP.MT1DInversion.create_synthetic(ρ =  ρ[2:end], zboundaries = zboundaries,
                                            freqs = 1 ./T, rseed=125, showplot=true, noisefrac=0.05, logscaledepth=false)
## now apply a stretch prior
# remember there must be a ρlow and ρhigh at every zall
ρlow, ρhigh = -0.5*ones(size(zall)), zall*0.00005 .+ 3.15
Δ = ρhigh - ρlow
F = transD_GP.MT1DInversion.makestretchop(F, ρlow=ρlow, Δ=Δ)
ax = gcf().axes[1]
transD_GP.MT1DInversion.plotpriorenv(F, ax=ax, lc = "r")
F.stretch = false # if wanting to toggle stretch priors

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
using Arrow, Tables, Test

include(joinpath(dirname(pathof(Arrow)), "../test/arrowjson.jl"))
# using .ArrowJSON

function runcommand(jsonname, arrowname, mode, verbose)
    if jsonname == ""
        error("must provide json file name")
    end
    if arrowname == ""
        error("must provide arrow file name")
    end

    if mode == "ARROW_TO_JSON"
        tbl = Arrow.Table(arrowname)
        df = ArrowJSON.DataFile(tbl)
        open(jsonname, "w") do io
            JSON3.write(io, df)
        end
    elseif mode == "JSON_TO_ARROW"
        df = ArrowJSON.parsefile(jsonname)
        open(arrowname, "w") do io
            Arrow.write(io, df)
        end
    elseif mode == "VALIDATE"
        df = ArrowJSON.parsefile(jsonname)
        tbl = Arrow.Table(arrowname)
        @test isequal(df, tbl)
    else
        error("unknown integration test mode: $mode")
    end
    return
end

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
export
	benchmark_solve!,
	shift_fill!

function interp_rows(N::Int,tf::Float64,X::AbstractMatrix)::Matrix
    n,N1 = size(X)
    t1 = range(0,stop=tf,length=N1)
    t2 = collect(range(0,stop=tf,length=N))
    X2 = zeros(n,N)
    for i = 1:n
        interp_cubic = CubicSplineInterpolation(t1, X[i,:])
        X2[i,:] = interp_cubic(t2)
    end
    return X2
end

function convertInf!(A::VecOrMat{Float64},infbnd=1.1e20)
    infs = isinf.(A)
    A[infs] = sign.(A[infs])*infbnd
    return nothing
end

function set_logger()
    if !(global_logger() isa SolverLogger)
        global_logger(SolverLogging.default_logger(true))
    end
end

function benchmark_solve!(solver; samples=10, evals=10)
    Z0 = deepcopy(get_trajectory(solver))
    if is_constrained(solver)
        λ0 = deepcopy(get_duals(solver))
    else
        λ0 = nothing
    end
    v0 = solver.opts.verbose
    s0 = solver.opts.show_summary
    solver.opts.verbose = 0
    solver.opts.show_summary = false
    b = @benchmark begin
        TO.initial_trajectory!($solver,$Z0)
        # set_duals!($solver, $λ0)
        solve!($solver)
    end samples=samples evals=evals
    solver.opts.verbose = v0
    solver.opts.show_summary = s0
    return b
end

function benchmark_solve!(solver, data::Dict; samples=10, evals=10)
	b = benchmark_solve!(solver, samples=samples, evals=evals)

   # Run stats
   push!(data[:time], time(median(b))*1e-6)  # ms
   push!(data[:iterations], iterations(solver))
   push!(data[:cost], cost(solver))
   return b
end

function shift_fill!(A::Vector, n=1, mode=:copy)
    N = length(A)
    N <= 1 && return nothing
	@inbounds for k = n+1:N
		A[k-n] = A[k]
    end
    if mode == :copy
        a_last = A[N-n]
        @inbounds for k = N-n:N
            A[k] = copy(a_last)
        end
    else mode == :zero
        @inbounds for k = N-n:N
            A[k] .= zero(A[k])
        end
    end
	return nothing
end


text:
# "Learning from Time-Changing Data with Adaptive Windowing"
#
# - v.0.3 (28 Jul 2017): Preparing the journal version of Pccf
#
# - v.0.2 (14 Dec 2016) : While doing experiments for change detection in multi-sensor settings ./201612-pccf3-code
#
#    Fixed the bug while checking the main condition
#
# - v.0.1 (12 Dec 2016) : While doing experiments for change detection in multi-sensor settings ./201612-pccf3-code
#
#    First implementation of `adwin' to use in multi-sensor settings
#    with Kalman Filter.
#
# https://github.com/abifet/adwin/blob/master/Java/ADWIN.java
#
using Distributions
#using Plots
#plotlyjs()

function scale01(sig)
    mx = maximum(sig)
    mn = minimum(sig)
    (sig - mn) / (mx-mn)
end

"Calc. ecut parameter for the partition of W = W0 dot W1"
function calcEcut(n0, n1, delta)
    @assert n0 >= 1
    @assert n1 >= 1
    m = 1.0/(1.0/n0 + 1.0/n1)
    deltaPrime = delta/(n0+n1)
    ecut = sqrt((1.0/(2.0*m)) * log(4/deltaPrime))
end

"Eq. 3.1"
function calcEcutEq31(w, n0, n1, delta)
    @assert n0 >= 1
    @assert n1 >= 1
    @assert n0 + n1 == length(w)
    sigma = std(w)
    m = 1.0/(1.0/n0 + 1.0/n1)
    deltaPrime = delta/(n0+n1)
    A = sqrt( (2/m) * sigma^2 * log(2/deltaPrime) )
    B = (2.0/(3.0 * m)) * log(2/deltaPrime)
    ecut = A + B
end

function recReduce!(x::Array{Float64,1},
                     ecut_custom::Float64,
                     min_len::Int=3,
                     rec_calls::Int = 1)
    n = length(x)
    if n <= min_len
        return x
    end
    cond_flag = true
    for k = 2:n-1
        cond = abs(mean(x[1:k])-mean(x[(k+1):end])) < ecut_custom
        if !cond
            cond_flag = false
            break
        end
    end
    if cond_flag
        return x
    else
        recReduce!(x[1:(end-1)], ecut_custom, min_len, rec_calls+1)
    end
end

function adwin_naive(sig::Array{Float64,1}; ecut_custom=0.5, min_len=10)
    ys = scale01(sig)
    n = length(ys)
    w::Array{Float64,1} = [ys[1]]
    r::Array{Float64,1}= [ys[1]]
    ws::Array{Int64,1} = [1] # widths
    for t=2:n
        insert!(w, 1, ys[t])
        w = recReduce!(w, ecut_custom, min_len)
        push!(ws, length(w))
    end
    ws
end

function genSigWithSlope(N, k=0.7)
    slope_part1 = randn(N)  .+ [1:N;] * k
    slope_part2 = randn(N)  .- [1:N;] * k
    sig = [rand(N);
           slope_part1;
           randn(N) .+ maximum(slope_part1);
           randn(floor(Int, 1.5*N))
           #slope_part2;
           #randn(N) .+ minimum(slope_part2)
           ]
    sig
end

"read output in inotebook and plot"
function demo_adwin_naive()
    srand(1234)
    N = 15
    settings = Dict("set1" => Dict("step" => 5.0, "ecut" => 0.5),
                    "set2" => Dict("step" => 2.5, "ecut" => 0.5/1.7),
                    "set3" => Dict("step" => 1.0, "ecut" => 0.5/1.7),
                    "set4_slope_sig" => Dict("step" => 1.0, "ecut" => 0.5/2.5)
                    )

    usesettings = settings["set4_slope_sig"]

    insig = scale01([randn(N); randn(N) + usesettings["step"]; randn(N)])
    insig2 = genSigWithSlope(N, 0.3)

    dat = insig2
    r = adwin_naive(dat, min_len=2, ecut_custom = usesettings["ecut"])

    # scale01(r), dat
    r, dat
end
# r = demo_adwin_naive()

function extractChps(x)
    chps::Array{Int,1}=[]
    n=length(x)
    for t=3:n
        if x[t] < x[t-1] && x[t-1] >= x[t-2]
            push!(chps, t)
        end
    end
    chps
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
using Documenter
using MCPhylo

makedocs(
    modules = MCPhylo,
    sitename="MCPhylo.jl",
    format = Documenter.HTML(),
    pages = [
        "Index" => "index.md",
        "Likelihood" => "likelihood.md",
        "Distributions" => "distributions.md",
        "Model" => "model.md",
        "Output" => "output.md",
        "Samplers" => "samplers.md",
        "Utils" => "utils.md",
        "Links" => "links.md",
        "For Developers" => "developers.md"],
   )

deploydocs(
    repo = "github.com/erathorn/MCPhylo.jl.git",
    devbranch="main"
     )

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
simpleloop = """
        x = 1/2*x
        y = 2y
    end
"""

# macro petter(n)
#     loop = """while true
#         x = x + y^$n
#         y = y + 1
#     end
#     """
#     eval(Expr(:(=), Symbol("petter$n"), loop))
# end

petter(n) = """
    while true
        x = x + y^$n
        y = y + 1
    end
"""

cohencu = """
    while n<=N
        n = n+1
        x = x+y
        y = y+z
        z = z+6
    end
"""

freire1 = """
    while x>r
        x = x-r
        r = r+1
    end
"""

freire2 = """
    while x-s > 0
        x = x-s
        s = s+6*r+3
        r = r+1
    end
"""

paperloop = """
    while true
        if true
            r = r - v
            v = v + 2
        else
            r = r + u
            u = u + 2
        end
    end
"""


text:
using AutoDiff
using MNIST

type SimpleParameter <: Parameter
	data
	timestamp
end

function SimpleParameter(data)
	return SimpleParameter(data,now())
end

type SimpleGradient <: Gradient
	data
end

function update(p::SimpleParameter, g::SimpleGradient)
	# update p with parameter g
	p.data-=g.data
	p.timestamp=now()
end

function zero_gradient(p::SimpleParameter)
	data = zeros( size(p.data, 1), size(p.data, 2) );
	return SimpleGradient(data);
end

#Replace error with an arbitrary function to minimize
#datum should be a tuple (digit,label) where label is length 10 array with a single 1
function loss(params::Array,datum)
	weights1=params[1]
	biases1=params[2]

	layer1=weights1*datum[1].+biases1
	tmp=exp(layer1)
	prediction=tmp./sum(tmp)

	error=sum((prediction-datum[2]).^2)
end

function loss(params::SimpleParameter, datum)
	loss(params.data,datum)
end
srand(1)
#Dummy inputs used for initializing variables in the gradient computation
dummy_weights1=0.001*randn((10,784))
dummy_biases1=0.001*randn((10,))
dummy_input=zeros(Float64,(784,))
dummy_output=zeros(Float64,(10,))

#assume dataset is a list of indices
function compute_gradient(state::WorkerState,dataset)
	println("[WORKER] Computing gradients")
	state.update_params(state.current_params.data)
	g=0
	for i in dataset
		example=trainfeatures(i)
		label=map(x->if x==trainlabel(i); 1.0; else 0.0; end, 0:9)
		g+=state.compute_gradient((example,label))
	end
	sleep(0.1)
	return SimpleGradient(state.learning_rate*g/length(dataset))
end


text:
using Flux, DiffEqFlux, DifferentialEquations, Plots

## Setup ODE to optimize
function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0] #True Params
prob = ODEProblem(lotka_volterra,u0,tspan,p)

# Verify ODE solution
sol = solve(prob,Tsit5())
plot(sol)

# Generate data from the ODE
data_sol = solve(prob,Tsit5(),saveat=0.1)
A1 = data_sol[1,:] # length 101 vector
A2 = data_sol[2,:] # length 101 vector
t = 0:0.1:10.0
scatter!(t,A1,color=[1],label = "rabbits")
scatter!(t,A2,color=[2],label = "wolves")

# Build a neural network that sets the cost as the difference from the
# generated data and true data

p = param([4., 1.0, 2.0, 0.4]) # Initial Parameter Vector
function predict_rd() # Our 1-layer neural network
  diffeq_rd(p,prob,Tsit5(),saveat=0.1)
end
loss_rd() = sum(abs2,predict_rd()-data_sol) # loss function

# Optimize the parameters so the ODE's solution stays near 1

data = Iterators.repeated((), 1000)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  #= display(loss_rd()) =#
  # using `remake` to re-create our `prob` with current parameters `p`
  scatter(t,A1,color=[1],label = "rabbit data")
  scatter!(t,A2,color=[2],label = "wolves data")
  display(plot!(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6),labels = ["rabbit model","wolf model"],color=[1 2]))
end
# Display the ODE with the initial parameter values.
cb()
Flux.train!(loss_rd, [p], data, opt, cb = cb)

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


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

println("--> Dataset 4")
config = "test_dataset4_config.txt"
p4 = run_NEAT_epochs(config, 2000, false)
print_best_solution_perfomance(p4)

#=
plot_neat_fitness(p4, imagename="fitness-dt4.png")
plot_neat_IP(p4, imagename="IP-dt4.png")
plot_neat_species(p4, imagename="species-dt4.png")
plot_output_distribution(p4, imagename="distribution-dt4.png")
plot_complexity(p4, imagename="solution-complexity-dt4.png")
plot_hit_rate(p4, imagename="hit-rate-dt4.png")
=#
#pop = getPopulation("checkpoint-dt4.jld2")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text: using SymbolicUtils


text:
using Pythia
cxx"""
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
"""

jetdef = icxx"""
  double Rparam = 0.4;
  fastjet::Strategy               strategy = fastjet::Best;
  fastjet::RecombinationScheme    recombScheme = fastjet::E_scheme;
  new fastjet::JetDefinition(fastjet::kt_algorithm, Rparam,
                                      recombScheme, strategy);
"""
fjinputs = icxx"new std::vector <fastjet::PseudoJet>";

pythia = newPythia()

icxx"""$pythia->readString("Beams:idA = 2212");"""
icxx"""$pythia->readString("Beams:idB = 2212");"""
icxx"""$pythia->readString("Beams:eCM = 8000.");"""
icxx"""$pythia->readString("HardQCD:gg2qqbar = on");"""
icxx"""$pythia->readString("PhaseSpace:pTHatMin = 50.0");"""
icxx"""$pythia->init();"""

function cluster(jetDef,e)
    icxx"$fjInputs->resize(0);"
    for p in particles(e)
        # Final State only
        if !icxx"$p.isFinal();"
            continue
        end

        # No neutrinos
        idabs = abs(icxx"$p.id();")
        if idabs == 12 || idabs == 14 || idabs == 16
            continue
        end

        if abs(icxx"$p.eta();" > 3.6)
            continue
        end

        icxx"$fjInputs->push_back( fastjet::PseudoJet( $p.px(),
        $p.py(), $p.pz(), $p.e() ) );"
    end
    icxx"new fastjet::ClusterSequence(*$fjinputs, *$jetDef);"
end

function jets(jetDef,e)
    clusterSeq = cluster(jetDef,e)
    sortedJets = icxx"""
        vector <fastjet::PseudoJet> inclusiveJets;
        inclusiveJets = $clusterSeq->inclusive_jets(20.0);
        new vector <fastjet::PseudoJet> (sorted_by_pt(inclusiveJets));
    """
end

for e in take(events(pythia),10)
    icxx"$pythia->process->list()"
    sortedJets = jets(jetDef,e)
    icxx"""
        for(unsigned int j = 0; j < $sortedJets->size(); j++) {
        cout << "Jet #" << j << " (px,py,pz,E) = (" << (*$sortedJets)[j].px() << ", " << (*$sortedJets)[j].py() << ", " << (*$sortedJets)[j].pz() << ", "
        << (*$sortedJets)[j].e() << ")\n";
        }
    """
end


text:
abstract type AbstractAggregation <: AbstractFunction end
getdim(::AbstractAggregation) = 1

@params struct NonNegSumOfSquares <: AbstractAggregation
    c::Real
end
(f::NonNegSumOfSquares)(x::AbstractVector) = f.c * sum(max.(0, x).^2)

@params struct NonNegSumOfRoots <: AbstractAggregation
    c::Real
end
(f::NonNegSumOfRoots)(x::AbstractVector) = f.c * sum(sqrt.(max.(0, x)))

@params struct NonNegSum <: AbstractAggregation
    c::Real
end
(f::NonNegSum)(x::AbstractVector) = f.c * sum(max.(0, x))

@params struct WeightedSum <: AbstractAggregation
    weights::AbstractVector
end
(f::WeightedSum)(x::AbstractVector) = dot(f.weights, x)


text:
julia> reverse("hey")
"yeh"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestCompositeParameters

using Mimi
using Test

@defcomp A begin
    p1 = Parameter(unit = "\$", default=3)
    p2 = Parameter()

    v1 = Variable()
    v2 = Variable(index=[time])
end

@defcomp B begin
    p1 = Parameter(unit = "thous \$", default=3)
    p3 = Parameter(description="B p3")
    p4 = Parameter(index=[time])

end

#------------------------------------------------------------------------------
# Test a failure to resolve namespace collisions

fail_expr1 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)
        # Will fail because you need to resolve the namespace collision of the p1's
    end
)

err1 = try eval(fail_expr1) catch err err end
@test occursin("unresolved parameter name collisions from subcomponents", sprint(showerror, err1))

#------------------------------------------------------------------------------
# Test a failure to resolve conflicting "description" fields

fail_expr2 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p1 = Parameter(A.p1, B.p1, unit="")
        p2 = Parameter(A.p2, B.p3)  # Will fail because the conflicting descriptions aren't resolved
    end
)

err2 = try eval(fail_expr2) catch err err end
@test occursin("subcomponents have conflicting values for the \"description\" field", sprint(showerror, err2))

#------------------------------------------------------------------------------
# Test trying to join parameters with different dimensions

fail_expr3 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p1 = Parameter(A.p1, B.p1, unit="")
        p2 = Parameter(A.p2, B.p4)  # Will fail because different dimensions
    end
)

err3 = try eval(fail_expr3) catch err err end
@test occursin("subcomponents have conflicting values for the \"dim_names\" field", sprint(showerror, err3))

#------------------------------------------------------------------------------
# Test a failure to auto-import a paramter because it's name has already been used

fail_expr4 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p3 = Parameter(A.p1, B.p1, unit="")    # should fail on auto import of other p3 parameter
    end
)

err4 = try eval(fail_expr4) catch err err end
@test occursin("this name has already been defined in the composite component's namespace", sprint(showerror, err4))

#------------------------------------------------------------------------------
# Test an attempt to import a parameter twice

fail_expr5 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p1 = Parameter(A.p1, B.p1, unit="")
        p1_repeat = Parameter(B.p1) # should not allow a second import of B.p1 (already connected)
    end
)

err5 = try eval(fail_expr5) catch err err end
@test occursin("Duplicate import", sprint(showerror, err5))

#------------------------------------------------------------------------------
# Test set_param! with unit collision

function get_model()
    m = Model()
    set_dimension!(m, :time, 10)
    add_comp!(m, A)
    add_comp!(m, B)
    return m
end

m1 = get_model()
err6 = try set_param!(m1, :p1, 5) catch err err end
@test occursin("components have conflicting values for the :unit field of this parameter", sprint(showerror, err6))

# use ignoreunits flag
set_param!(m1, :p1, 5, ignoreunits=true)

err7 = try run(m1) catch err err end
@test occursin("Cannot build model; the following parameters are not set", sprint(showerror, err7))

# Set separate values for p1 in A and B
m2 = get_model()
set_param!(m2, :A, :p1, 1)  # Set the value only for component A
@test length(m2.md.external_param_conns) == 1 # test that only one connection has been made
@test Mimi.UnnamedReference(:B, :p1) in Mimi.unconnected_params(m2.md)  # and that B.p1 is still unconnected

err8 = try set_param!(m2, :B, :p1, 2) catch err err end
@test occursin("the model already has an external parameter with this name", sprint(showerror, err8))

set_param!(m2, :B, :p1, :B_p1, 2)   # Use a unique name to set B.p1
@test length(m2.md.external_param_conns) == 2
@test Set(keys(m2.md.external_params)) == Set([:p1, :B_p1])

end


text:
using Documenter
using DSPG_Manual

# @testset "Documentation" begin
    using Documenter, DSPG_Manual

    DocMeta.setdocmeta!(DSPG_Manual,
                       :DocTestSetup,
                       :(using DSPG_Manual, RCall),
                       recursive = true)
    makedocs(sitename = "DSPG Manual",
             modules = [DSPG_Manual],
             pages = [
                 "Home" => "index.md",
                 "Onboarding" => "Onboarding.md",
                 "Data Science Framework" => "DataScienceFramework.md",
                 "Infrastructure" => [
                     "Zotero" => joinpath("Infrastructure", "Zotero.md"),
                     "GitHub" => joinpath("Infrastructure", "GitHub.md"),
                     "Rivanna" => [
                        "Dashboard" => joinpath("Infrastructure", "Rivanna", "Dashboard.md"),
                        "File Explorer" => joinpath("Infrastructure", "Rivanna", "FileExplorer.md"),
                        "Configuration" => joinpath("Infrastructure", "Rivanna", "Config.md"),
                        "Interactive Apps" => joinpath("Infrastructure", "Rivanna", "InteractiveApps.md"),
                        "Database" => joinpath("Infrastructure", "Rivanna", "Database.md"),
                     ]
                 ],
             ],
             source = joinpath(@__DIR__, "..", "docs", "src"),
             build = joinpath(@__DIR__, "..", "docs", "build"),
            #  format = Documenter.Writers.HTMLWriter.HTML(prettyurls = parse(Bool, get(ENV, "CI", "false")),)
             )
# end

max_stars_repo_name: uva-bi-sdad/DSPG-Manual

max_issues_repo_name: uva-bi-sdad/DSPG-Manual

max_forks_repo_name: uva-bi-sdad/DSPG-Manual


text:
export @mkdeepconvert, @mkdeepconvert1, @mkdeepconvert2, @bi_str, @bf_str, @bigint, @int128
export deepbigint2, deepbigfloat
export @bigfloat

@mkdeepconvert1(deepbigint, BigInt)
macro bi_str(s) deepbigint(s) end
_deepbigfloat(x::Real) = parse(BigFloat, string(x))
@mkdeepconvert(deepbigfloat,_deepbigfloat)
macro bf_str(s) deepbigfloat(s) end
macro bigint(s)
    return quote
        $(esc(deepbigint2(s)))
    end
end
macro bigfloat(s)
    return quote
        $(esc(deepbigfloat(s)))
    end
end

@mkdeepconvert3(deepbigint2, BigInt, Int)
# hmm, this does not work
@mkdeepconvert3(deepbigfloat2, BigFloat, Real)

macro int128(s)
    return quote
        $(esc(deepint128(s)))
    end
end

@mkdeepconvert3(deepint128, Int128, Int)


text:
# getindex.jl
# test

using LinearMapsAA: LinearMapAA, LinearMapAX
using Test: @test, @testset

"""
    LinearMapAA_test_getindex(A::LinearMapAX)
tests for `getindex`
"""
function LinearMapAA_test_getindex(A::LinearMapAX)
    B = Matrix(A)
    @test all(size(A) .>= (4,4)) # required by tests

    tf1 = [false; trues(size(A,1)-1)]
    tf2 = [false; trues(size(A,2)-2); false]
    ii1 = (3, 2:4, [2,4], :, tf1)
    ii2 = (2, 3:4, [1,4], :, tf2)
    for i2 in ii2
        for i1 in ii1
            @test B[i1,i2] == A[i1,i2]
        end
    end

    L = A._lmap
    test_adj = !((:fc in propertynames(L)) && isnothing(L.fc))
    if test_adj
        for i1 in ii2
            for i2 in ii1
                @test B'[i1,i2] == A'[i1,i2]
            end
        end
    end

    # "end"
    @test B[3,end-1] == A[3,end-1]
    @test B[end-2,3] == A[end-2,3]
    if test_adj
        @test B'[3,end-1] == A'[3,end-1]
    end

    # [?]
    @test B[1] == A[1]
    @test B[7] == A[7]
    if test_adj
        @test B'[3] == A'[3]
    end
    @test B[end] == A[end] # lastindex

    kk = [k in [3,5] for k = 1:length(A)] # Bool
    @test B[kk] == A[kk]

    # Some tests could rely on the fact that LinearMapAM <:i AbstractMatrix,
    # by inheriting from general Base.getindex, but all are provided here.
    @test B[[1, 3, 4]] == A[[1, 3, 4]]
    @test B[4:7] == A[4:7]

    true
end

N = 6; M = 8 # non-square to stress test
forw = x -> [cumsum(x); 0; 0]
back = y -> reverse(cumsum(reverse(y[1:N])))
A = LinearMapAA(forw, back, (M, N))
@test LinearMapAA_test_getindex(A)


text:
# used in shortest path calculations
# has_distances{T}(distmx::AbstractArray{T,2}) =
#     issparse(distmx)? (nnz(distmx) > 0) : !isempty(distmx)

type DefaultDistance<:AbstractArray{Int, 2}
    nv::Int
    DefaultDistance(nv::Int=typemax(Int)) = new(nv)
end

getindex(::DefaultDistance, s::Int, d::Int) = 1
getindex(::DefaultDistance, s::UnitRange, d::UnitRange) = DefaultDistance(length(s))
size(d::DefaultDistance) = (d.nv, d.nv)
transpose(d::DefaultDistance) = d
ctranspose(d::DefaultDistance) = d

"""
Calculates the eccentricity[ies] of a vertex `v`, vertex vector `vs`, or the
entire graph. An optional matrix of edge distances may be supplied.

The eccentricity of a vertex is the maximum shortest-path distance between it
and all other vertices in the graph.

Because this function must calculate shortest paths for all vertices supplied
in the argument list, it may take a long time.

The output is either a single float (when a single vertex is provided) or a
vector of floats corresponding to the vertex vector. If no vertex vector is
provided, the vector returned corresponds to each vertex in the graph.

Note: the eccentricity vector returned by `eccentricity()` may be used as input
for the rest of the distance measures below. If an eccentricity vector is
provided, it will be used. Otherwise, an eccentricity vector will be calculated
for each call to the function. It may therefore be more efficient to calculate,
store, and pass the eccentricities if multiple distance measures are desired.
"""
function eccentricity{T}(
    g::SimpleGraph,
    v::Int,
    distmx::AbstractArray{T, 2} = DefaultDistance()
)
    e = maximum(dijkstra_shortest_paths(g,v,distmx).dists)
    e == typemax(T) && error("Infinite path length detected")

    return e
end

eccentricity{T}(
    g::SimpleGraph,
    vs::AbstractArray{Int, 1}=vertices(g),
    distmx::AbstractArray{T, 2} = DefaultDistance()
) =
    [eccentricity(g,v,distmx) for v in vs]

eccentricity{T}(g::SimpleGraph, distmx::AbstractArray{T, 2}) =
    eccentricity(g, vertices(g), distmx)

"""Returns the maximum eccentricity of the graph."""
diameter{T}(all_e::Vector{T}) = maximum(all_e)
diameter{T}(g::SimpleGraph, distmx::AbstractArray{T, 2} = DefaultDistance()) =
    maximum(eccentricity(g, distmx))

"""Returns the set of all vertices whose eccentricity is equal to the graph's
diameter (that is, the set of vertices with the largest eccentricity).
"""
function periphery{T}(all_e::Vector{T})

    diam = maximum(all_e)
    return filter((x)->all_e[x] == diam, 1:length(all_e))
end

periphery{T}(g::SimpleGraph, distmx::AbstractArray{T, 2} = DefaultDistance()) =
    periphery(eccentricity(g, distmx))

"""Returns the minimum eccentricity of the graph."""
radius{T}(all_e::Vector{T}) = minimum(all_e)
radius{T}(g::SimpleGraph, distmx::AbstractArray{T, 2} = DefaultDistance()) =
    minimum(eccentricity(g, distmx))

"""Returns the set of all vertices whose eccentricity is equal to the graph's
radius (that is, the set of vertices with the smallest eccentricity).
"""
function center{T}(all_e::Vector{T})
    rad = radius(all_e)
    return filter((x)->all_e[x] == rad, 1:length(all_e))
end

center{T}(g::SimpleGraph, distmx::AbstractArray{T, 2} = DefaultDistance()) =
    center(eccentricity(g, distmx))


text:
function miyak00()
#
#  test polynomial suggested by Goedecker
#  square of Fibocacci polynomial
#
    p = poly([(1.1+1.1im)*ones(4);(3.2+2.3im)*[1;1];2.1+1.5im]);
    z = [[(1.1+1.1im);(3.2+2.3*im);2.1+1.5*im] [4;2;1]];
    p, PolyZeros(z)
end


text:
mutable struct FockProblem{F,O,C,M,KWS}
    fock::F
    P::O
    c::C
    H::M
    kws::KWS
end

function set!(f::FockProblem, v)
    copyto!(f.P, v)
    update!(f.fock)
end

function (f::FockProblem)()
    energy_matrix!(f.H, f.fock)
    real(f.c'f.H*f.c)
end

function (f::FockProblem)(v)
    set!(f, v)
    f()
end

function jac!(w, v, f::FockProblem)
    set!(f, v)
    # Threads.@threads
    for j = 1:length(f.fock.equations)
        mul!(view(w,:,j), f.kws[j], view(v,:,j))
    end
end


text:
 type _TreeState
   leaveIx::Int
   innerIx::Int
 end

 # a,b is the span of the subtree
 # terminates when a == b
 function construct_tree(parents::Vector{Int}, a::Int, b::Int)
   function _construct_tree(parents, a, b, state::_TreeState)
       #println("Split $a and $b with state $state")
       if a == b then
           state.leaveIx = state.leaveIx + 1
           return state.leaveIx
       end

       # Construct subtrees
       local split = rand(0:(b - a - 1))
       local leftIx = _construct_tree(parents, a, a + split, state)
       local rightIx = _construct_tree(parents, a + split + 1, b, state)

       state.innerIx = state.innerIx + 1
       # println("Node $(state.innerIx) ($a-$b), children $leftIx and $rightIx")
       parents[leftIx] = -state.innerIx
       parents[rightIx] = state.innerIx
       return state.innerIx
   end

   parents[size(parents, 1)] = 0
   return _construct_tree(parents, a, b, _TreeState(0, div(size(parents, 1) + 1, 2)))
 end

 #######

function generate_bhs(nleaves, input_size)
  parents = Array(Int, 2 * nleaves - 1)
  construct_tree(parents, 1, nleaves)
  hs = BinaryHierarchicalSoftmax{CPUDevice, Float64}(input_size, parents)
end

input_size = 100
hs = generate_bhs(200000, input_size)
init_gradient!(hs)

words = reshape([i for i=1:hs.nleaves], hs.nleaves, 1)
x = randn(1, input_size)
input = repmat(x, hs.nleaves)
output = forward!(hs, Any[input, words])
@assert abs(sum(exp(output)) - 1) / hs.nleaves <= 1e-16

# -- Timeit

ninput = 10000
input = randn(ninput, input_size)
targets = rand(1:hs.nleaves, ninput, 1)
println("Size of input = $(size(input)), # leaves = $(hs.nleaves)")
for i=1:10
  @time forward!(hs, Any[input, targets])
end

gradOutput = fill(-1., (ninput, 1)) # We want to increase sum(o) => decrease -sum(o)
sgd = StochasticGradient(1e-3)

inputs = Any[input, targets]
@time for i = 1:10
  init_gradient!(hs)
  forward!(hs, inputs)
  backward!(hs, inputs, gradOutput)
  optimize!(sgd, hs)
  println("Cost = $(sum(hs.output))")
end


text:
"""
    add_hs(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Adds Hydrogens to the molecule.

```julia
add_hs(mol)
```
"""
function add_hs(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:add_hs, librdkitcffi), Cshort, (Ref{Cstring}, Ref{Csize_t}), mol.mol, mol.mol_size)
end

"""
    remove_all_hs(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Removes all Hydrogens from the molecule.

```julia
remove_all_hs(mol)
```
"""
function remove_all_hs(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:remove_all_hs, librdkitcffi), Cshort, (Ref{Cstring}, Ref{Csize_t}), mol.mol, mol.mol_size)
end


text:
# ------------ Trace typing rules ------------ #

# Define an interpreter - allows multiple dispatch to handle multiple abstract interpretations.
struct TraceTypingInterpreter <: InterpretationContext end

# Fallback.
absint(ctx::TraceTypingInterpreter, args...) = Union{}
absint(ctx::TraceTypingInterpreter, trace, args...) = Missing

# Distributions.
absint(ctx::TraceTypingInterpreter, x::D, args...) where D <: Distributions.Distribution =  x
function absint(ctx::TraceTypingInterpreter, ::Type{D}, args...) where D <: Distributions.Distribution
    all(map(args) do a
            _lit(a)
        end) || return D
    return D(args...)
end

# Note - the usage of the @abstract macro below is strictly equivalent to the longform definitions for absint with ctx::T above.

# Learnables.
@abstract TraceTypingInterpreter learnable(::Symbol) = Any

# Calls to trace.
@abstract TraceTypingInterpreter trace(::Symbol, ::Type{DiscreteUniform}, args...) = Integers{1}
@abstract TraceTypingInterpreter trace(::Symbol, d::DiscreteUniform, args...) = DiscreteInterval{d.a, d.b}
@abstract TraceTypingInterpreter trace(::Symbol, ::Type{Uniform}, args...) = Reals{1}
@abstract TraceTypingInterpreter trace(::Symbol, d::Uniform, args...) = RealInterval(d.a, d.b)
@abstract TraceTypingInterpreter trace(::Symbol, ::Type{Gamma}, args...) = PositiveReals{1}
@abstract TraceTypingInterpreter trace(::Symbol, ::Gamma, args...) = PositiveReals{1}
@abstract TraceTypingInterpreter trace(::Symbol, ::Type{Normal}, args...) = Reals{1}
@abstract TraceTypingInterpreter trace(::Symbol, ::Normal, args...) = Reals{1}
@abstract TraceTypingInterpreter trace(::Symbol, ::Type{Bernoulli}, args...) = Discrete{2}
@abstract TraceTypingInterpreter trace(::Symbol, ::Bernoulli, args...) = Discrete{2}
@abstract TraceTypingInterpreter trace(::Symbol, jfn::JFunction, args...) = get_trace_type(jfn)

# Combinators.
@abstract TraceTypingInterpreter function trace(::Symbol, u::Gen.Unfold, args...)
    tt = get_trace_type(u.kernel)
    List{tt}
end
@abstract TraceTypingInterpreter function trace(::Symbol, u::Gen.Map, args...)
    List{get_trace_type(u.value.kernel)}
end
@abstract TraceTypingInterpreter function trace(::Symbol, u::Gen.Switch, args...)
    Sum{Tuple(map(u.value.branches) do br
                  get_trace_type(br)
              end)}
end


text:
using Seahawk
using Documenter

makedocs(;
    modules=[Seahawk],
    authors="Ole Kröger <o.kroeger@opensourc.es> and contributors",
    repo="https://github.com/Wikunia/Seahawk.jl/blob/{commit}{path}#L{line}",
    sitename="Seahawk.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Wikunia.github.io/Seahawk.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Wikunia/Seahawk.jl",
)


text:
using Documenter

push!(LOAD_PATH,  "../../src")

using Stipple, Stipple.Elements, Stipple.Layout, Stipple.Typography

makedocs(
    sitename = "Stipple - data dashboards and reactive UIs for Julia",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "Stipple API" => [
          "Elements" => "api/elements.md",
          "Layout" => "api/layout.md",
          "NamedTuples" => "api/namedtuples.md",
          "Stipple" => "api/stipple.md",
          "Typography" => "api/typography.md",
        ]
    ],
)

deploydocs(
  repo = "github.com/GenieFramework/Stipple.jl.git",
)


text:
using Robotlib
h = 0.004032;

dh = DHYuMi()
xi = DH2twistsPOE(dh)

pathopen = "/work/fredrikb/extRosetta/base2rob2xyz.csv"
pathsave = "/tmp/fredrikb/log.mat"

data    = orcalog2mat(pathopen, pathsave)
data    = readmat(pathsave)
ds      = 1
q       = getData("robot_1.*posRawAbs", data, ds)
q̇       = getData("robot_1.*velFlt", data, ds)
τ       = getData("robot_1.*trqRaw", data, ds)
f       = getData("force", data, ds)
# f[:,[2,5]] *= -1

abb2logical!(q)
abb2logical!(q̇)
abb2logical!(τ)
q = q*dh.GR'
q̇ = q̇*dh.GR'
τ = τ*inv(dh.GR')

q̈ = filtfilt(ones(50),[50.],smartDiff(q̇))

baseAnglesLeft  = [-0.63 , 0.95 , -0.18]
Rbase           = rpy2R(baseAnglesLeft,"xyz")

N   = size(q,1)
T  = cat(3,[fkinePOE(xi,q[i,:]') for i = 1:N]...);
v1 = T[1:3,4,1500] - T[1:3,4,500]
v2 = T[1:3,4,2200] - T[1:3,4,1500]
v3 = T[1:3,4,3000] - T[1:3,4,2200]

v1 /= norm(v1)
v2 /= norm(v2)
v3 /= norm(v3)
Rbase = [v1 v2 v3]'


text:
# This file is a part of UltraNest.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    UltraNest

Julia wrapper for Python nested sampling package
[UltraNest](https://github.com/JohannesBuchner/UltraNest).
"""
module UltraNest

using LinearAlgebra
using Random
using Statistics

using PyCall

"""
    const ultranest

The Python `ultranest` module.

Example:

```julia
using UltraNest
smplr = ultranest.ReactiveNestedSampler(paramnames, my_likelihood, kwargs...)
result = smplr.run()
```

See the the
[UltraNest Python documentation](https://johannesbuchner.github.io/UltraNest/)
regarding usage.

    !!! note

    Convention for matrices holding multiple parameter vectors (resp.
    multiple samples): In UltraNest.jl, using Julia's column-major array
    indexing, parameter vectors are stored as rows (not columns) in the
    matrices.
"""
const ultranest = PyNULL()
export ultranest

const scipy_stats = PyNULL()
const numpy = PyNULL()

include("ultra_nest.jl")

function __init__()
    copy!(ultranest, pyimport_conda("ultranest", "ultranest", "conda-forge"))
    copy!(scipy_stats, pyimport_conda("scipy.stats", "scipy"))
    copy!(numpy, pyimport("numpy"))
end

end # module


text:
module Cfdsof

# Write your package code here

include("readOpenFoamMesh.jl")
include("writeOpenFoamMesh.jl")
include("fields.jl")
include("boundaryConditions.jl")
include("solver.jl")
include("output.jl")

# Case general information
caseName = "elbow"

# Mesh
meshPath = "src/testCases/$caseName/constant/polyMesh"
mesh = readOpenFoamMesh(meshPath)

# Setup fields
k = setupScalarField("k", 1000)
T = setupScalarField("T")

# Boundary conditions
assignBC("velocity-inlet-5", T, "fixedValue", 200)
assignBC("velocity-inlet-6", T, "fixedValue", 100)
assignBC("pressure-outlet-7", T, "fixedValue", 500)
assignBC("wall-4", T, "zeroGradient")
assignBC("wall-8", T, "zeroGradient")
assignBC("frontAndBackPlanes", T, "zeroGradient")

# Solve the equations
solveLaplacian(k, T)

# Write outputs
writeVtu("src/testCases/vtu/$caseName", T, k)

end