text:
using CSV, DataFrames

testnum = 49

filename = string("/Users/crismoen/RunToSolve/Papers/Moen_2020_CFSRC_PurlinDesigner", "/section_export", "/Test", testnum, "nodes", ".csv")
nodes = CSV.read(filename, DataFrame, header=false)

path = "/Users/crismoen/.julia/dev/StructuresKit/docs/PurlinDesigner/papers/CFSRC2020/inputs/"

filename = string(path, "julia_sectionprops.csv")
sectprops = CSV.read(filename, DataFrame, header=false)

filename = string(path, "/Gao_Moen_2013_Table1", ".csv")
dims = CSV.read(filename, DataFrame, header=false)

I = sectprops[testnum-1, 1]

Ma = 1.0
Mb = 2.0

Yc = sectprops[testnum-1,7]   #from bottom fiber

#calculate stress
σa = Ma.*(nodes[:,3] .- Yc)./I
σb = Mb.*(nodes[:,3] .- Yc)./I

dσ = σb .- σa

plot(nodes[:,2], nodes[:,3], dσ)
plot!(nodes[:,2], nodes[:,3], zeros(length(q)))

dx = 1.0

q = zeros(Float64, numnodes)
s = zeros(Float64, numnodes)

for i = 1:numelem

    x1 = nodes[i, 2]
    y1 = nodes[i, 3]
    x2 = nodes[i+1, 2]
    y2 = nodes[i+1, 3]

    dS = norm(x2-x1, y2-y1)
    dA = dS*t

    dq = (dσ[i]*dA + (1/2*(dσ[i]+dσ[i+1]) * dA))/dx

    q[i+1] = q[i] + dq
    s[i+1] = s[i] + dS

end

plot(nodes[:,2], nodes[:,3], q)
plot!(nodes[:,2], nodes[:,3], zeros(length(q)))

plot(s, q)

t = dims[testnum-1,11]

numnodes = length(nodes[:,1])
numelem = numnodes - 1

q = zeros(Float64, numnodes)
s = zeros(Float64, numnodes)

using LinearAlgebra
using Statistics

for i = 1:numelem

    x1 = nodes[i, 2]
    y1 = nodes[i, 3]
    x2 = nodes[i+1, 2]
    y2 = nodes[i+1, 3]

    dS = norm(x2-x1, y2-y1)
    y = Yc - mean([y1, y2])
    dq = V/I*y*t*dS

    q[i+1] = q[i] + dq
    s[i+1] = s[i] + dS

end

using Plots
plot(s, q)

plot(nodes[:,2], nodes[:,3], q)
plot!(nodes[:,2], nodes[:,3], zeros(length(q)))


text:
#!/usr/bin/env julia
#
# @license Apache-2.0
#
# Copyright (c) 2018 The Stdlib Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import JSON

"""
    gen( x, y, name )

Generate fixture data and write to file.

# Arguments

* `x`: base domain
* `y`: exponent domain
* `name::AbstractString`: output filename

# Examples

``` julia
julia> x = range( -1000, stop = 1000, length = 2001 );
julia> y = range( 0, stop = 1000, length = 1001 );
julia> gen( x, y, \"data.json\" );
```
"""
function gen( x, y, name )
	z = Array{Float64}( undef, length(x) );
	for i in eachindex(x)
		z[ i ] = x[i]^y[i];
	end

	# Store data to be written to file as a collection:
	data = Dict([
		("x", x),
		("y", y),
		("expected", z)
	]);

	# Based on the script directory, create an output filepath:
	filepath = joinpath( dir, name );

	# Write the data to the output filepath as JSON:
	outfile = open( filepath, "w" );
	write( outfile, JSON.json(data) );
	close( outfile );
end

# Get the filename:
file = @__FILE__;

# Extract the directory in which this file resides:
dir = dirname( file );

# Squares (x small):
x = range( -1000, stop = 1000, length = 2001 );
y = zeros( 2001 );
y[:] .= 2.0;
gen( x, y, "squared_small.json" );

# Squares (x large):
x = range( 3.14e100, stop = 3.14e110, length = 2001 );
y = zeros( 2001 );
y[:] .= 2.0;
gen( x, y, "squared_large.json" );

# Cubes (x small):
x = range( -1000, stop = 1000, length = 2001 );
y = zeros( 2001 );
y[:] .= 3.0;
gen( x, y, "cubed_small.json" );

# Cubes (x large):
x = range( -9.99e99, stop = -4.641588833612779e102, length = 2001 );
y = zeros( 2001 );
y[:] .= 3.0;
gen( x, y, "cubed_large.json" );

# To the 4th power (x small):
x = range( -1000, stop = 1000, length = 2001 );
y = zeros( 2001 );
y[:] .= 4.0;
gen( x, y, "to_the_fourth_small.json" );

# To the 4th power (x large):
x = range( 1.234e61, stop = 3.14e69, length = 2001 );
y = zeros( 2001 );
y[:] .= 4.0;
gen( x, y, "to_the_fourth_large.json" );

# Base near unity (y small):
x = range( 0.9999990463256836, stop = 1.0000009536743164, length = 500 );
y = range( -10, stop = 10, length = 500 );
gen( x, y, "base_near_unity_small.json" );

# Base near unity (y large):
x = range( 0.9999990463256836, stop = 1.0000009536743164, length = 500 );
y = range( 100, stop = 709, length = 500 );
gen( x, y, "base_near_unity_large.json" );

# Base near unity (y huge):
x = range( 0.9999999999990463256836, stop = 1.0000000000009536743164, length = 500 );
y = range( 1.8446744073709552e13, stop = 4.294967296e14, length = 500 );
gen( x, y, "base_near_unity_huge.json" );

# Subnormal results:
z = range( 3.14e-309, stop = 3.14e-312, length = 2001 );
x = range( 3.14e-52, stop = 3.14e-54, length = 2001 );
y = log.(z) ./ log.(x);
gen( x, y, "subnormal_results.json" );

# Negative exponents (x small, y small):
x = range( 0, stop = 10, length = 5001 );
y = rand( 5001 ) .* -100;
gen( x, y, "negative_exp_small_small.json" );

# Negative exponents (x large, y small):
x = range( 3.14e10, stop = 1e15, length = 5001 );
y = rand( 5001 ) .* -100;
gen( x, y, "negative_exp_large_small.json" );

# Negative exponents (x small, y large):
x = range( 0, stop = 10, length = 5001 );
y = rand( 5001 ) .* -(2.0^31);
gen( x, y, "negative_exp_small_large.json" );

# Negative exponents (x large, y large):
x = range( 3.14e20, stop = 1e22, length = 5001 );
y = rand( 5001 ) .* -(2.0^31);
gen( x, y, "negative_exp_large_large.json" );

# x small, y small:
x = rand( 5001 ) .* 100;
y = rand( 5001 ) .* 100;
gen( x, y, "small_small.json" );

# x small, y large:
x = rand( 5001 ) .* 1.9988548118735103;
y = rand( 5001 ) .* 2.0^10;
gen( x, y, "small_large.json" );

# x large, y small:
x = rand( 5001 ) .* 5e20;
y = rand( 5001 ) .* 10;
gen( x, y, "large_small.json" );

# Random (x decimal, y decimal):
x = rand( 5001 ) .* 500;
y = rand( 5001 ) .* 100;
gen( x, y, "decimal_decimal.json" );

# Random (x decimal, y integer):
x = rand( 5001 ) .* 500;
y = round.( rand(5001) .* 100 );
gen( x, y, "decimal_integer.json" );

# Random (x integer, y decimal):
x = round.( rand(5001) .* 500 );
y = rand( 5001 ) .* 100;
gen( x, y, "integer_decimal.json" );

# Random (x integer, y integer):
x = round.( rand(5001) .* 500 );
y = round.( rand(5001) .* 100 );
gen( x, y, "integer_integer.json" );

ext: jl

max_stars_repo_name: ghalimi/stdlib

max_issues_repo_name: ghalimi/stdlib

max_forks_repo_name: sthagen/stdlib


text:
using Comonicon
using JLD2
using LinearAlgebra
using BetaVQE.VAN
using Yao
using Yao.EasyBuild
using BetaVQE
using Optimisers: ADAM
using StatsBase
using Random

# file storage
function build_key(nx, ny, Γ, β, depth, nsamples, nhiddens, lr; folder="data")
    mkpath(folder)
    key = "tfim"
    key *= "_nx$nx"
    key *= "_ny$ny"
    key *= "_Gamma$Γ"
    key *= "_beta$β"
    key *= "_d$depth"
    key *= "_batch$nsamples"
    key *= "_lr$lr"
    key *= "_nhiddens"

    for h in nhiddens
        key = key * "_$h"
    end
    return joinpath(folder, key)
end

@cast function exact_spectra(nx::Int, ny::Int, Γ::Float64; folder=joinpath("data", "exact"))
    mkpath(folder)
    h = hamiltonian(TFIM(nx, ny; Γ=Γ, periodic=false))
    H = mat(h)
    w, _ = eigen(Matrix(H))
    save(joinpath(folder, "nx$nx"*"ny$ny"*"Gamma$Γ"*".jld2"), "spectra", w)
end

@cast function scan_beta(nx::Int=2, ny::Int=2, Γ::Float64=1.0;
                         depth::Int=5, nsamples::Int=1000, nhiddens::Vector{Int}=[500], lr::Float64=0.01, niter::Int=500, cont::Bool=false)

    for β in collect(0.1:0.1:1.0)
        learn(nx, ny, Γ, β; depth=depth, nsamples=nsamples, nhiddens=nhiddens, lr=lr, niter=niter, cont=cont)
    end

end

@cast function scan_gamma(nx::Int=2, ny::Int=2, β::Float64=1.0;
                         depth::Int=5, nsamples::Int=1000, nhiddens::Vector{Int}=[500], lr::Float64=0.01, niter::Int=500, cont::Bool=false)

    for Γ in collect(0.0:1.0:4.0)
        learn(nx, ny, Γ, β; depth=depth, nsamples=nsamples, nhiddens=nhiddens, lr=lr, niter=niter, cont=cont)
    end

end

@cast function learn(nx::Int=2, ny::Int=2, Γ::Float64=1.0, β::Float64=1.0;
                  depth::Int=5, nsamples::Int=1000, nhiddens::Vector{Int}=[500], lr::Float64=0.01, niter::Int=500, cont::Bool=false)

    Random.seed!(42)
    key = build_key(nx, ny, Γ, β, depth, nsamples, nhiddens, lr; folder=joinpath("data", "tns3"))
    println(key)

    nbits = nx*ny

    if nhiddens[1] == 0
        network = PSAModel(nbits)
    else
        network = AutoRegressiveModel(nbits, nhiddens)
    end

    circuit = tns_circuit(nbits, depth, EasyBuild.pair_square(nx, ny; periodic=false); entangler=(n,i,j)->put(n,(i,j)=>general_U4()))

    h = hamiltonian(TFIM(nx, ny; Γ=Γ, periodic=false))
    F_exact, E_exact, S_exact, Cv_exact, γ_exact = BetaVQE.exact(β, h)

    chkp_file = key*".jld2"
    if cont && isfile(chkp_file)
        chkp = load(chkp_file)
        model_dispatch!(network, chkp["cparams"])
        dispatch!(circuit, chkp["qparams"])
        println("load chkp from $chkp_file")
    else
        message = "# $F_exact"
        println(message)

        logfile = open(key*".log", "w")
        write(logfile, message*"\n")
        close(logfile)
    end

    logfile = open(key*".log", "a")
    cparams, qparams = BetaVQE.train(β, h, network, circuit, logfile; optimizer=ADAM(lr), nbatch=nsamples, niter=niter)
    close(logfile)

    #compute observables
    samples = gen_samples(network, nsamples)
    E = BetaVQE.energy(h, network, circuit, samples)
    E2 = BetaVQE.energy2(h, network, circuit, samples)
    S = BetaVQE.entropy(h, network, circuit, samples)
    γ = BetaVQE.purity(h, network, circuit, samples)
    F = E - S/β

    Cv = β^2*(E2 - E^2)

    save(key*".jld2", "cparams", cparams, "qparams", qparams, "exact", (F_exact, E_exact, S_exact, Cv_exact, γ_exact), "result", (F, E, S, Cv, γ))
end

@cast function inference(nx::Int=2, ny::Int=2, Γ::Float64=1.0, β::Float64=1.0;
                  depth::Int=5, nsamples::Int=1000, nhiddens::Vector{Int}=[500], lr::Float64=0.01, niter::Int=500)

    Random.seed!(42)
    key = build_key(nx, ny, Γ, β, depth, nsamples, nhiddens, lr; folder=joinpath("data", "tns3"))
    println(key)

    nbits = nx*ny

    if nhiddens[1] == 0
        network = PSAModel(nbits)
    else
        network = AutoRegressiveModel(nbits, nhiddens)
    end

    circuit = tns_circuit(nbits, depth, EasyBuild.pair_square(nx, ny; periodic=false); entangler=(n,i,j)->put(n,(i,j)=>general_U4()))

    h = hamiltonian(TFIM(nx, ny; Γ=Γ, periodic=false))

    chkp_file = key*".jld2"
    if isfile(chkp_file)
        chkp = load(chkp_file)
        model_dispatch!(network, chkp["cparams"])
        dispatch!(circuit, chkp["qparams"])
        println("load chkp from $chkp_file")
    else
        throw(ArgumentError("file not exists: $chkp_file"))
    end

    samples = gen_samples(network, nsamples)
    samples = unique(samples, dims=2)
    s = BetaVQE.spectra(h, network, circuit, samples)
    jldopen(chkp_file, "a+") do file
        file["spectra"] = s
    end
end

@main


text:
struct IterationOutput{M, R}
    solution::M
    value::R
    gradient::M
end

to_iterator(x::Any) = x
to_iterator(n::T) where T <: Number = Iterators.repeated(n)

include("optimizers/gradient_descent.jl")
include("optimizers/polyak.jl")
include("optimizers/nesterov.jl")
include("optimizers/barzilai_borwein.jl")
include("optimizers/adagrad.jl")
include("optimizers/rmsprop.jl")
include("optimizers/adadelta.jl")
include("optimizers/adam.jl")


text:
################################################################################
# Displacement Jacobian
################################################################################
@testset "Displacement Jacobian" begin
    @testset "Rotational" begin
        mech = Dojo.get_mechanism(:pendulum)
        joint0 = mech.joints[1]
        rot0 = joint0.rotational
        rot0.axis_offset = rand(QuatRotation).q

        xa = rand(3)
        qa = rand(QuatRotation).q
        xb = rand(3)
        qb = rand(QuatRotation).q

        # displacement_jacobian_configuration
        X0, Q0 = Dojo.displacement_jacobian_configuration(:parent, rot0, xa, qa, xb, qb,
            attjac=true)
        J0 = [X0 Q0]
        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.displacement(rot0, z[1:3], Quaternion(z[4:7]...), xb, qb),
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # displacement_jacobian_configuration
        X0, Q0 = Dojo.displacement_jacobian_configuration(:child, rot0, xa, qa, xb, qb,
            attjac=true)
        J0 = [X0 Q0]
        attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.displacement(rot0, xa, qa, z[1:3], Quaternion(z[4:7]...)),
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8
    end

    @testset "Translational" begin
        mech = Dojo.get_mechanism(:slider)
        joint0 = mech.joints[1]
        tra0 = joint0.translational

        xa = rand(3)
        qa = rand(QuatRotation).q
        xb = rand(3)
        qb = rand(QuatRotation).q

        # displacement_jacobian_configuration
        X0, Q0 = Dojo.displacement_jacobian_configuration(:parent, tra0, xa, qa, xb, qb,
            attjac=true)
        J0 = [X0 Q0]
        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.displacement(tra0, z[1:3], Quaternion(z[4:7]...), xb, qb),
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # displacement_jacobian_configuration
        X0, Q0 = Dojo.displacement_jacobian_configuration(:child, tra0, xa, qa, xb, qb,
            attjac=true)
        J0 = [X0 Q0]
        attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.displacement(tra0, xa, qa, z[1:3], Quaternion(z[4:7]...)),
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8
    end
end

################################################################################
# Impulse transform Jacobian
################################################################################
@testset "Impulse transform Jacobian: rotational" begin
    @testset "Rotational" begin
        mech = Dojo.get_mechanism(:pendulum)
        joint0 = mech.joints[1]
        rot0 = joint0.rotational
        rot0.axis_offset = rand(QuatRotation).q

        xa = rand(3)
        qa = rand(QuatRotation).q
        xb = rand(3)
        qb = rand(QuatRotation).q

        p0 = rand(3)
        J0 = Dojo.impulse_transform_jacobian(:parent, :parent, rot0, xa, qa, xb, qb, p0)
        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:parent, rot0, z[1:3], Quaternion(z[4:7]...), xb, qb) * p0,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.impulse_transform_jacobian(:parent, :child, rot0, xa, qa, xb, qb, p0)
        attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:parent, rot0, xa, qa, z[1:3], Quaternion(z[4:7]...)) * p0,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.impulse_transform_jacobian(:child, :parent, rot0, xa, qa, xb, qb, p0)
        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:child, rot0, z[1:3], Quaternion(z[4:7]...), xb, qb) * p0,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.impulse_transform_jacobian(:child, :child, rot0, xa, qa, xb, qb, p0)
        attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:child, rot0, xa, qa, z[1:3], Quaternion(z[4:7]...)) * p0,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8
    end

    @testset "Translational" begin
        mech = Dojo.get_mechanism(:slider)
        joint0 = mech.joints[1]
        tra0 = joint0.translational

        xa = rand(3)
        qa = rand(QuatRotation).q
        xb = rand(3)
        qb = rand(QuatRotation).q

        p0 = rand(3)
        J0 = Dojo.impulse_transform_jacobian(:parent, :parent, tra0, xa, qa, xb, qb, p0)
        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:parent, tra0, z[1:3], Quaternion(z[4:7]...), xb, qb) * p0,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.impulse_transform_jacobian(:parent, :child, tra0, xa, qa, xb, qb, p0)
        attjac = cat(I(3), Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:parent, tra0, xa, qa, z[1:3], Quaternion(z[4:7]...)) * p0,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.impulse_transform_jacobian(:child, :parent, tra0, xa, qa, xb, qb, p0)
        attjac = cat(I(3), Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:child, tra0, z[1:3], Quaternion(z[4:7]...), xb, qb) * p0,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.impulse_transform_jacobian(:child, :child, tra0, xa, qa, xb, qb, p0)
        attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_transform(:child, tra0, xa, qa, z[1:3], Quaternion(z[4:7]...)) * p0,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8
    end
end

################################################################################
# Impulse map Jacobian
################################################################################
@testset "Impulse map" begin
    @testset "impulse map: rotational" begin
        mech = Dojo.get_mechanism(:pendulum)
        joint0 = mech.joints[1]
        rot0 = joint0.rotational
        rot0.axis_offset = rand(QuatRotation).q

        xa = rand(3)
        qa = rand(QuatRotation).q
        xb = rand(3)
        qb = rand(QuatRotation).q
        λ = rand(2)
        η = rand(0)
        Dojo.impulse_map(:parent, rot0, xa, qa, xb, qb, η)
        Dojo.impulse_map(:child, rot0, xa, qa, xb, qb, η)

        m = 1.0
        J = I(3)
        pbody = Body(m, J)
        cbody = Body(m, J)
        pbody.state.x2 = xa
        pbody.state.q2 = qa
        cbody.state.x2 = xb
        cbody.state.q2 = qb
        # impulse_map_parent_jacobian_parent
        J0 = Dojo.impulse_map_jacobian(:parent, :parent, rot0, pbody, cbody, λ)
        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:parent, rot0, z[1:3], Quaternion(z[4:7]...), xb, qb, 0) * λ,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # impulse_map_parent_jacobian_child
        J0 = Dojo.impulse_map_jacobian(:parent, :child, rot0, pbody, cbody, λ)
        attjac = cat(I(3), Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:parent, rot0, xa, qa, z[1:3], Quaternion(z[4:7]...), 0) * λ,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # impulse_map_child_jacobian_parent
        J0 = Dojo.impulse_map_jacobian(:child, :parent, rot0, pbody, cbody, λ)
        attjac = cat(I(3), Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:child, rot0, z[1:3], Quaternion(z[4:7]...), xb, qb, 0) * λ,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # impulse_map_child_jacobian_child
        J0 = Dojo.impulse_map_jacobian(:child, :child, rot0, pbody, cbody, λ)
        attjac = cat(I(3), Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:child, rot0, xa, qa, z[1:3], Quaternion(z[4:7]...), 0) * λ,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8
    end

    @testset "Translational" begin
        mech = get_mechanism(:pendulum)
        joint0 = mech.joints[1]
        tra0 = joint0.translational

        xa = rand(3)
        qa = rand(QuatRotation).q
        xb = rand(3)
        qb = rand(QuatRotation).q
        λ = rand(3)
        η = rand(0)
        Dojo.impulse_map(:parent, tra0, xa, qa, xb, qb, η)
        Dojo.impulse_map(:child, tra0, xa, qa, xb, qb, η)

        m = 1.0
        J = I(3)
        pbody = Body(m, J)
        cbody = Body(m, J)
        pbody.state.x2 = xa
        pbody.state.q2 = qa
        cbody.state.x2 = xb
        cbody.state.q2 = qb
        # impulse_map_parent_jacobian_parent
        J0 = Dojo.impulse_map_jacobian(:parent, :parent, tra0, pbody, cbody, λ)

        attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:parent, tra0, z[1:3], Quaternion(z[4:7]...), xb, qb, 0) * λ,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # impulse_map_parent_jacobian_child
        J0 = Dojo.impulse_map_jacobian(:parent, :child, tra0, pbody, cbody, λ)
        attjac = cat(I(3), Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:parent, tra0, xa, qa, z[1:3], Quaternion(z[4:7]...), 0) * λ,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # impulse_map_child_jacobian_parent
        J0 = Dojo.impulse_map_jacobian(:child, :parent, tra0, pbody, cbody, λ)
        attjac = cat(I(3), Dojo.LVᵀmat(qa), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:child, tra0, z[1:3], Quaternion(z[4:7]...), xb, qb, 0) * λ,
            [xa; Dojo.vector(qa)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8

        # impulse_map_child_jacobian_child
        J0 = Dojo.impulse_map_jacobian(:child, :child, tra0, pbody, cbody, λ)
        attjac = cat(I(3), Dojo.LVᵀmat(qb), dims=(1,2))
        J1 = ForwardDiff.jacobian(
            z -> Dojo.impulse_map(:child, tra0, xa, qa, z[1:3], Quaternion(z[4:7]...), 0) * λ,
            [xb; Dojo.vector(qb)]
            ) * attjac
        @test norm(J0 - J1, Inf) < 1.0e-8
    end
end


text:
################################################################################
# Correlation Functions
################################################################################
"""
	function correlate(u::Vector{T},v::Vector{T},offset::Integer,window::Integer = -1) where {T<:Number}

Calculates correlation between vectors u and v with given offset. Specify averaging window to limit range of correlation. If the window extends beyond the end of one vector, it treats out-of-bounds indices as zero.
"""
function correlate(u::Vector{T},v::Vector{T},offset::Integer,window::Integer = -1) where {T<:Number}

	@assert offset <= length(u) "Offset out of bounds"
	@assert window <= length(u) && window <= length(v) "Window must be smaller than input vector lengths"

	if window == -1
		window = length(u)
	end

	v1 = view(u,1:window)
	v2 = view(v,1+offset:min(window+offset,length(v)))
	if window+offset > length(v)
		v2 = vcat(v2,zeros(window+offset-length(v)))
	end

	return dot(v1,v2)/window
end

export correlate

"""
	function autocorrelate(u::Vector{T},offset::Integer, window::Integer = -1) where {T<:Number}

Calculates correlation of vector u with itself.
"""
function autocorrelate(u::Vector{T},offset::Integer, window::Integer = -1) where {T<:Number}
	correlate(u,u,offset,window)
end

export autocorrelate

"""
	ncorrelate(u::Vector{T},v::Vector{T},offset::Integer,window::Integer = -1) where {T<:Number}

Calculates normalized correlation between vectors u and v with given offset. Specify averaging window to limit range of correlation. If the window extends beyond the end of one vector, it treats out-of-bounds indices as zero.
"""
function ncorrelate(u::Vector{T},v::Vector{T},offset::Integer,window::Integer = -1) where {T<:Number}

	@assert offset <= length(u) "Offset out of bounds"
	@assert window <= length(u) && window <= length(v) "Window must be smaller than input vector lengths"

	if window == -1
		window = length(u)
	end

	v1 = view(u,1:window)
	v2 = view(v,1+offset:min(window+offset,length(v)))
	if window+offset > length(v)
		v2 = vcat(v2,zeros(window+offset-length(v)))
	end

	return dot(v1,v2)/(window*mean(v1)*mean(v2))
end

export ncorrelate
"""
    function corrTimes(τ::Vector,γCounts::Vector)

Returns an array of τ values for which the γCounts autocorrelation is non-zero.
"""
function corrTimes(τ::Vector,γCounts1::Vector,γCounts2::Vector)
    @assert length(γCounts1) == length(γCounts2) "Count vectors must have the same length"
    return τ[map(i->correlate(γCounts1,γCounts2,i,length(τ)) > 0 ? true : false,collect(0:length(τ)-1))]
end

export corrTimes

"""
    function autocorrTimes(τ::Vector,γCounts::Vector)

Returns an array of τ values for which the γCounts autocorrelation is non-zero.
"""
function autocorrTimes(τ::Vector,γCounts::Vector)
    return τ[map(i->autocorrelate(γCounts,i,length(τ)) > 0 ? true : false,collect(0:length(τ)-1))]
end

export autocorrTimes

"""
    γCoincidentτ(beam1::Vector{T}, beam2::Vector{T}) where {T<:Number}

Returns the separation index between the first coincident counts of beam 1 and 2.
Returns -1 if no coincident counts are calculated.
"""
function γCoincidentIndex(beam1::Vector{T}, beam2::Vector{T}) where {T<:Number}
    @assert length(beam1) == length(beam2) "Beam vectors must have equal length"
    # iterate through beam 1
    for (i,nbar) in enumerate(beam1)
        if poissonCount(nbar) > 0
            for (j,mbar) in enumerate(beam2[i:end])
                if poissonCount(mbar) > 0
                    # return the separation index between coincident counts
                    return j-1 # subtract 1 since that's where Julia indexing starts
                end
            end
            return -1
        end
    end
    return -1
end

export γCoincidentIndex

"""
    firstNonzero(a::Vector)

Returns the index of the first non-zero element in a.
Returns -1 if all values are zero.
"""
function firstNonzero(a::Vector)
    for (i,val) in enumerate(a)
        if val != 0
            return i
        end
    end
    return -1
end

export firstNonzero
################################################################################
# Counting Related Functions
################################################################################
"""
    function countDeltaTimes(τ::Vector,γCounts::Vector)

Returns an array of τ values for which the γCounts autocorrelation is non-zero.
"""
function countTimes(times::Vector,γCounts::Vector)
    out = Vector{Real}(undef,0)
    for (i,counts) in enumerate(γCounts)
        if counts != 0
            countTimes = times[i]*ones(counts)
            out = vcat(out,countTimes)
        end
    end
    return out
end

export countTimes

################################################################################
# Functions from MGST2021
################################################################################

"""
    stauAvg(τ::Number,params::eFieldParams,n::Integer)

Returns the average of the Doppler noise term
"""
function stauAvg(τ::Number,source::LightSource)
    return stauAvg(τ,source.σ,source.n)
end

"""
    stauAvg(τ::Number,σ::Number,n::Integer)

Returns the average of the Doppler noise term
"""
function stauAvg(τ::Number,σ::Number,n::Integer)
    term1 = n
    term2 = n*(n-1)
    term2 *= exp(-σ^2*τ^2)
    return term1 + term2
end

export stauAvg

"""
    stauVar(τ::Number,source::LightSource)

Returns the variance of the Doppler noise term
"""
function stauVar(τ::Number,source::LightSource)
    return stauVar(τ,source.σ,source.n)
end

"""
    stauVar(τ::Number,σ::Number,n::Integer)

Returns the variance of the Doppler noise term
"""
function stauVar(τ::Number,σ::Number,n::Integer)
    στ2 = σ^2*τ^2

    prod1 = 8*n*(n-1)
    prod1 *= exp(-2*στ2)

    prod2 = n-1+cosh(στ2)

    prod3 = sinh(στ2/2)^2

    return prod1*prod2*prod3
end

export stauVar

"""
    stau(τ::Number,field::eField)

Calculates the value of the Doppler noise term for the given eField
"""
function stau(τ::Number,field::eField)
    return stau(τ,field.νn)
end

"""
    stau(τ::Number,νn::Vector)

Calculates the value of the Doppler noise term for the given τ and frequencies
"""
function stau(τ::Number,νn::Vector)
    terms = exp.(-2.0*π*im*τ*νn)
    sumterms = sum(terms)
    return real(sumterms*conj(sumterms))
end

export stau

"""
    g2Calc(τ::Number,n::Integer,source::LightSource)

Returns the average calculated value of g2(τ) from MGST2021.
"""
function g2Calc(τ::Number,source::LightSource)
    n = source.n
    em2 = real.(source.Em .* conj.(source.Em))
    em4 = em2 .* em2
    sumEm2 = sum(em2)
    sumEm4 = sum(em4)

    g2τ = 1.0

    term2 = -sumEm4/(n*sumEm2^2)

    g2τ += term2

    term3 = sum(em2 .* exp.(-im*τ*Δm(source)))/sumEm2
    term3 *= conj(term3)
    term3 = real(term3)
    term3 *= stauAvg(τ,source)/n^2

    return g2τ+term3
end

export g2Calc


text:
module TestServices

using Bio.Services.EUtils
using EzXML
using Base.Test

@testset "EUtils" begin
    @testset "einfo" begin
        res = einfo(db="pubmed")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)
    end

    @testset "esearch" begin
        res = esearch(db="pubmed", term="asthma")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)
    end

    @testset "epost" begin
        ctx = Dict()
        res = epost(ctx, db="protein", id="NP_005537")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)
        @test haskey(ctx, :WebEnv)
        @test haskey(ctx, :query_key)
    end

    @testset "esummary" begin
        # esummary doesn't seem to support accession numbers
        res = esummary(db="protein", id="15718680,157427902,119703751")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)

        res = esummary(db="protein", id=["15718680", "157427902", "119703751"])
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)

        # esearch then esummary
        query = "asthma[mesh] AND leukotrienes[mesh] AND 2009[pdat]"
        ctx = Dict()
        res = esearch(ctx, db="pubmed", term=query, usehistory=true, retmode="xml")
        @test res.status == 200
        res = esummary(ctx, db="pubmed")
        @test res.status == 200

        ctx = Dict()
        res = esearch(ctx, db="pubmed", term=query, usehistory=true, retmode="json")
        @test res.status == 200
        res = esummary(ctx, db="pubmed")
        @test res.status == 200
    end

    @testset "efetch" begin
        res = efetch(db="nuccore", id="NM_001178.5", retmode="xml", idtype="acc")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)

        # epost then efetch
        ctx = Dict()
        res = epost(ctx, db="protein", id="NP_005537")
        @test res.status == 200
        res = efetch(ctx, db="protein", retmode="xml")
        @test res.status == 200
    end

    @testset "elink" begin
        res = elink(dbfrom="protein", db="gene", id="NM_001178.5")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)
    end

    @testset "egquery" begin
        res = egquery(term="asthma")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)
    end

    @testset "espell" begin
        res = espell(db="pubmed", term="athma")
        @test res.status == 200
        @test startswith(res.headers["Content-Type"], "text/xml")
        @test isa(parsexml(res.data), EzXML.Document)
    end

    @testset "ecitmatch" begin
        res = ecitmatch(
            db="pubmed",
            retmode="xml",
            bdata="proc natl acad sci u s a|1991|88|3248|mann bj|Art1|")
        @test res.status == 200
        # ECitMatch is mysterious: it returns a plain text data even though it
        # requires retmode="xml".
        # @show res.headers["Content-Type"]
    end
end

end # TestServices


text:
"""
    get_image(x)
Get an image
"""
function get_image(x)

    y = load(x);
    return y

end


text:
using Turing
using Distributions
using Test

@model gdemo(x) = begin
  s ~ InverseGamma(2, 3)
  # println(s)
  m ~ Normal(0, sqrt(s))
  # println(m)
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  # println(vi)
  return s, m
end

alg = HMC(3000, 0.15, 10)

res = sample(gdemo([1.5, 2.0]), alg)

println(mean(res[:s])," ≈ ", 49/24, "?")
println(mean(res[:m])," ≈ ", 7/6, "?")
@test mean(res[:s]) ≈ 49/24 atol=0.2
@test mean(res[:m]) ≈ 7/6 atol=0.2


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

import Zygote
using Zygote: @adjoint

using Flux, DiffEqFlux
using OrdinaryDiffEq
using DiffEqCallbacks

using DiffEqFlux: ODEFunction, basic_tgrad, ODEProblem, ZygoteVJP, InterpolatingAdjoint, solve

import ForwardDiff
import Optim

import SciMLBase: RightRootFind

using FMIImport: fmi2ComponentStateModelInitialized, fmi2ComponentStateModelSetableFMUstate, fmi2ComponentStateModelUnderEvaluation
using FMIImport: fmi2Struct, FMU2, FMU2Component

include("FMI2_neural.jl")

"""
The mutable struct representing an abstract (simulation mode unknown) NeuralFMU.
"""
abstract type NeuralFMU end

"""
Structure definition for a NeuralFMU, that runs in mode `Model Exchange` (ME).
"""
mutable struct ME_NeuralFMU <: NeuralFMU
    neuralODE::NeuralODE
    solution::ODESolution
    fmu::FMU

    tspan
    saveat
    saved_values
    recordValues

    valueStack

    handleEvents::Bool

    ME_NeuralFMU() = new()
end

"""
Structure definition for a NeuralFMU, that runs in mode `Co-Simulation` (CS).
"""
mutable struct CS_NeuralFMU{T} <: NeuralFMU
    model
    fmu::T

    tspan
    saveat
    valueStack

    CS_NeuralFMU{T}() where {T} = new{T}()
end

##### EVENT HANDLING START

# Read next time event from fmu and provide it to the integrator
function time_choice(c::FMU2Component, integrator)
    eventInfo = fmi2NewDiscreteStates(c)
    fmi2EnterContinuousTimeMode(c)
    if eventInfo.nextEventTimeDefined == fmi2True
        Zygote.@ignore @debug "time_choice(_, _): Next event defined at $(eventInfo.nextEventTime)s"
        return eventInfo.nextEventTime
    end

    return Inf # ToDo: better `nothing`
end

# Handles events and returns the values and nominals of the changed continuous states.
function handleEvents(c::FMU2Component, enterEventMode::Bool, exitInContinuousMode::Bool)

    if enterEventMode
        fmi2EnterEventMode(c)
    end

    eventInfo = fmi2NewDiscreteStates(c)

    valuesOfContinuousStatesChanged = eventInfo.valuesOfContinuousStatesChanged
    nominalsOfContinuousStatesChanged = eventInfo.nominalsOfContinuousStatesChanged

    #set inputs here
    #fmiSetReal(myFMU, InputRef, Value)

    while eventInfo.newDiscreteStatesNeeded == fmi2True
        # update discrete states
        eventInfo = fmi2NewDiscreteStates(c)
        valuesOfContinuousStatesChanged = eventInfo.valuesOfContinuousStatesChanged
        nominalsOfContinuousStatesChanged = eventInfo.nominalsOfContinuousStatesChanged

        if eventInfo.terminateSimulation == fmi2True
            Zygote.@ignore @error "Event info returned error!"
        end
    end

    if exitInContinuousMode
        fmi2EnterContinuousTimeMode(c)
    end

    return valuesOfContinuousStatesChanged, nominalsOfContinuousStatesChanged
end

# Returns the event indicators for an FMU.
function condition(nfmu::ME_NeuralFMU, out, x, t, integrator) # Event when event_f(u,t) == 0

    # ToDo: set inputs here
    #fmiSetReal(myFMU, InputRef, Value)

    if isa(t, ForwardDiff.Dual)
        t = ForwardDiff.value(t)
    end

    if all(isa.(x, ForwardDiff.Dual))
        x = collect(ForwardDiff.value(e) for e in x)
    end

    fmi2SetTime(nfmu.fmu, t)
    # ToDo: Evaluate on light-weight model (sub-model) without fmi2GetXXX or similar and the bottom ANN
    nfmu.neuralODE.model(x) # evaluate NeuralFMU (set new states)

    indicators = fmi2GetEventIndicators(nfmu.fmu)
    Zygote.@ignore @debug "condition(_, _, $(x), $(t), _): Event indicators=$(indicators)"

    copy!(out, indicators)
end

function f_optim(x, nfmu, right_x_fmu) # , idx, direction::Real)
    # propagete the new state-guess `x` through the NeuralFMU
    nfmu.neuralODE.model(x)
    indicators = fmi2GetEventIndicators(nfmu.fmu)
    return Flux.Losses.mse(right_x_fmu, fmi2GetContinuousStates(nfmu.fmu)) # - min(-direction*indicators[idx], 0.0)
end

# Handles the upcoming events.
function affectFMU!(nfmu::ME_NeuralFMU, integrator, idx)
    # Event found - handle it

    t = integrator.t
    if isa(t, ForwardDiff.Dual)
        t = ForwardDiff.value(t)
    end

    x = integrator.u
    if all(isa.(x, ForwardDiff.Dual))
        x = collect(ForwardDiff.value(e) for e in x)
    end

    left_x_fmu = fmi2GetContinuousStates(nfmu.fmu)

    continuousStatesChanged, nominalsChanged = handleEvents(nfmu.fmu.components[end], true, Bool(sign(idx)))

    Zygote.@ignore @debug "affectFMU!(_, _, $idx): Event [$idx] detected at t=$(t)s (statesChanged=$(continuousStatesChanged))"
    #indicators = fmi2GetEventIndicators(nfmu.fmu)

    if continuousStatesChanged == fmi2True

        left_x = x

        right_x_fmu = fmi2GetContinuousStates(nfmu.fmu) # the new FMU state after handled events

        Zygote.@ignore @debug "affectFMU!(_, _, $idx): NeuralFMU state event from $(left_x) (fmu: $(left_x_fmu)). Indicator [$idx]: $(indicators[idx]). Optimizing new state ..."

        # ToDo: Problem-related parameterization of optimize-call
        #result = optimize(x_seek -> f_optim(x_seek, nfmu, right_x_fmu), left_x, LBFGS(); autodiff = :forward)
        #result = Optim.optimize(x_seek -> f_optim(x_seek, nfmu, right_x_fmu, idx, sign(indicators[idx])), left_x, NelderMead())
        result = Optim.optimize(x_seek -> f_optim(x_seek, nfmu, right_x_fmu), left_x, NelderMead())

        #display(result)

        right_x = Optim.minimizer(result)
        integrator.u = right_x

        Zygote.@ignore @debug "affectFMU!(_, _, $idx): NeuralFMU state event to   $(right_x) (fmu: $(right_x_fmu)). Indicator [$idx]: $(fmi2GetEventIndicators(nfmu.fmu)[idx]). Minimum: $(Optim.minimum(result))."
    end

    if nominalsChanged == fmi2True
        x_nom = fmi2GetNominalsOfContinuousStates(nfmu.fmu)
    end
end

# Does one step in the simulation.
function stepCompleted(nfmu::ME_NeuralFMU, x, t, integrator)
    (status, enterEventMode, terminateSimulation) = fmi2CompletedIntegratorStep(nfmu.fmu, fmi2True)
    if enterEventMode == fmi2True
        affectFMU!(nfmu, integrator, 0)
    else
        # ToDo: set inputs here
        #fmiSetReal(myFMU, InputRef, Value)
    end
end

# save FMU values
function saveValues(c::FMU2Component, recordValues, u, t, integrator)

    if isa(t, ForwardDiff.Dual)
        t = ForwardDiff.value(t)
    end

    if all(isa.(x, ForwardDiff.Dual))
        x = collect(ForwardDiff.value(e) for e in x)
    end

    fmi2SetTime(c, t)
    # ToDo: Evaluate on light-weight model (sub-model) without fmi2GetXXX or similar and the bottom ANN
    nfmu.neuralODE.model(x) # evaluate NeuralFMU (set new states)

    (fmiGetReal(c, recordValues)...,)
end

##### EVENT HANDLING END

"""
Constructs a ME-NeuralFMU where the FMU is at an arbitrary location inside of the NN.

# Arguents
    - `fmu` the considered FMU inside the NN
    - `model` the NN topology (e.g. Flux.chain)
    - `tspan` simulation time span
    - `alg` a numerical ODE solver

# Keyword arguments
    - `saveat` time points to save the NeuralFMU output, if empty, solver step size is used (may be non-equidistant)
    - `fixstep` forces fixed step integration
    - `recordFMUValues` additionally records internal FMU variables (currently not supported because of open issues)
"""
function ME_NeuralFMU(fmu::FMU2,
                      model,
                      tspan,
                      alg=nothing;
                      saveat=[],
                      recordFMUValues = nothing,
                      kwargs...)

    nfmu = ME_NeuralFMU()
    nfmu.fmu = fmu

    ext_model = nothing

    ext_model = Chain(model.layers...)

    nfmu.handleEvents = false
    nfmu.saved_values = nothing

    if (nfmu.fmu.modelDescription.numberOfEventIndicators > 0)
        @info "This ME-NeuralFMU has event indicators. Event-handling is in BETA-testing, so it is disabled by default. If you want to try event-handling during training for this discontinuous FMU, use the attribute `myNeuralFMU.handleEvents=true`."
    end

    nfmu.recordValues = prepareValueReference(fmu, recordFMUValues)

    nfmu.neuralODE = NeuralODE(ext_model, tspan, alg; saveat=saveat, kwargs...)

    nfmu.tspan = tspan
    nfmu.saveat = saveat

    nfmu
end

"""
Constructs a CS-NeuralFMU where the FMU is at an arbitrary location inside of the NN.

# Arguents
    - `fmu` the considered FMU inside the NN
    - `model` the NN topology (e.g. Flux.chain)
    - `tspan` simulation time span

# Keyword arguments
    - `saveat` time points to save the NeuralFMU output, if empty, solver step size is used (may be non-equidistant)
"""
function CS_NeuralFMU(fmu,
                      model,
                      tspan;
                      saveat=[],
                      recordValues = [])

    nfmu = CS_NeuralFMU{typeof(fmu)}()

    nfmu.fmu = fmu

    nfmu.model = Chain(model.layers...)

    nfmu.tspan = tspan
    nfmu.saveat = saveat

    nfmu
end

"""
Evaluates the ME_NeuralFMU in the timespan given during construction or in a custom timespan from `t_start` to `t_stop` for a given start state `x_start`.

# Keyword arguments
    - `reset`, the FMU is reset every time evaluation is started (default=`true`).
    - `setup`, the FMU is set up every time evaluation is started (default=`true`).
    - `rootSearchInterpolationPoints`, the number of interpolation points during root search before event-handling (default=`100`). Try to increase this value, if event-handling fails with exception.
"""
function (nfmu::ME_NeuralFMU)(x_start::Union{Array{<:Real}, Nothing} = nothing,
                              t_start::Real = nfmu.tspan[1],
                              t_stop = nfmu.tspan[end];
                              reset::Bool=true,
                              setup::Bool=true,
                              rootSearchInterpolationPoints::Integer=100,
                              kwargs...)

    c = nfmu.fmu.components[end]
    callbacks = []
    sense = nothing
    savedValues = nothing
    x0 = nothing

    if reset
        if c.state == fmi2ComponentStateModelInitialized
            fmi2Terminate(c)
        end
        if c.state == fmi2ComponentStateModelSetableFMUstate
            fmi2Reset(c)
        end
    end

    if setup
        fmi2SetupExperiment(nfmu.fmu, t_start, t_stop)
        fmi2EnterInitializationMode(nfmu.fmu)

        Zygote.ignore() do
            # First evaluation of the FMU
            if x_start == nothing
                x0 = fmi2GetContinuousStates(c)
            else
                x0 = x_start
            end
            x0_nom = fmi2GetNominalsOfContinuousStates(c)

            fmi2SetContinuousStates(c, x0)
        end

        fmi2ExitInitializationMode(nfmu.fmu)
    end

    tspan = getfield(nfmu.neuralODE,:tspan)
    t_start = tspan[1]
    t_stop = tspan[end]

    #################

    saving = (length(nfmu.recordValues) > 0)

    eventHandling = (nfmu.fmu.modelDescription.numberOfEventIndicators > 0)
    handleTimeEvents = false

    Zygote.ignore() do
        if nfmu.handleEvents && eventHandling
            eventInfo = fmi2NewDiscreteStates(c)
            handleTimeEvents = (eventInfo.nextEventTimeDefined == fmi2True)
        end

        handleEvents(c, false, false)

        # Get states of handling initial Events
        x0 = fmi2GetContinuousStates(c)
        x0_nom = fmi2GetNominalsOfContinuousStates(c)

        fmi2EnterContinuousTimeMode(c)
    end

    Zygote.@ignore @debug "NeuralFMU experimental event handling: $(nfmu.handleEvents) stateEvents: $(eventHandling) timeEvents: $(handleTimeEvents)"

    Zygote.ignore() do
        # Event handling only if there are event-indicators
        if nfmu.handleEvents && eventHandling

            # Only ForwardDiffSensitivity supports differentiation and callbacks
            sense = ForwardDiffSensitivity(;chunk_size=0,convert_tspan=true)

            eventCb = VectorContinuousCallback((out, x, t, integrator) -> condition(nfmu, out, x, t, integrator),
                                            (integrator, idx) -> affectFMU!(nfmu, integrator, idx),
                                            Int64(nfmu.fmu.modelDescription.numberOfEventIndicators);
                                            rootfind=RightRootFind,
                                            save_positions=(false, false),
                                            interp_points=rootSearchInterpolationPoints)
            push!(callbacks, eventCb)

            stepCb = FunctionCallingCallback((x, t, integrator) -> stepCompleted(nfmu, x, t, integrator);
                                            func_everystep=true,
                                            func_start=true)
            push!(callbacks, stepCb)

            if handleTimeEvents
                timeEventCb = IterativeCallback((integrator) -> time_choice(c, integrator),
                                                (integrator) -> affectFMU!(nfmu, integrator, 0),
                                                Float64;
                                                initial_affect = false,
                                                save_positions=(false,false))

                push!(callbacks, timeEventCb)
            end
        else
            sense = InterpolatingAdjoint(autojacvec=ZygoteVJP())
        end

        if saving
            savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(nfmu.recordValues))...})

            savingCB = SavingCallback((u,t,integrator) -> saveValues(c, nfmu.recordValues, u, t, integrator),
                              savedValues,
                              saveat=nfmu.saveat)
            push!(callbacks, savingCB)
        end
    end

    p = nfmu.neuralODE.p
    dudt_(u,p,t) = nfmu.neuralODE.re(p)(u)
    ff = ODEFunction{false}(dudt_,tgrad=basic_tgrad)
    prob = ODEProblem{false}(ff,x0,tspan,p)

    #sense = ForwardDiffSensitivity() # (;chunk_size=0,convert_tspan=true)
    #callbacks = []

    # Zygote.ignore() do
    #     @info "Tspan: $(tspan)"
    #     @info "Sensealg: $(sense)"
    #     @info "Callbacks: $(length(callbacks))"
    # end

    if length(callbacks) > 0
        nfmu.solution = solve(prob, nfmu.neuralODE.args...; sensealg=sense, saveat=nfmu.saveat, callback=CallbackSet(callbacks...), nfmu.neuralODE.kwargs...)
    else
        nfmu.solution = solve(prob, nfmu.neuralODE.args...; sensealg=sense, saveat=nfmu.saveat, nfmu.neuralODE.kwargs...)
    end

    if saving
        return nfmu.solution, savedValues
    else
        return nfmu.solution
    end
end

"""
Evaluates the CS_NeuralFMU in the timespan given during construction or in a custum timespan from `t_start` to `t_stop` with a given time step size `t_step`.

Via optional argument `reset`, the FMU is reset every time evaluation is started (default=`true`).
"""
function (nfmu::CS_NeuralFMU{T})(inputFct,
                                 t_step::Real,
                                 t_start::Real = nfmu.tspan[1],
                                 t_stop::Real = nfmu.tspan[end];
                                 reset::Bool = true,
                                 setup::Bool = true) where {T}

    c = nfmu.fmu.components[end]

    if reset
        if c.state == fmi2ComponentStateModelInitialized
            fmi2Terminate(c)
        end
        if c.state == fmi2ComponentStateModelSetableFMUstate
            fmi2Reset(c)
        end
    end

    if setup
        while fmi2SetupExperiment(nfmu.fmu, t_start, t_stop) != 0
        end
        while fmi2EnterInitializationMode(nfmu.fmu) != 0
        end
        while fmi2ExitInitializationMode(nfmu.fmu) != 0
        end
    end

    ts = t_start:t_step:t_stop
    nfmu.fmu.components[end].skipNextDoStep = true # skip first fim2DoStep-call
    model_input = inputFct.(ts)
    valueStack = nfmu.model.(model_input)
    return valueStack

    # savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(nfmu.model(inputFct(t_start))))...})

    # ts = t_start:t_step:t_stop
    # i = 1
    # for t in ts
    #     DiffEqCallbacks.copyat_or_push!(savedValues.t, i, t)
    #     DiffEqCallbacks.copyat_or_push!(savedValues.saveval, i, (nfmu.model(inputFct(t))...,), Val{false})
    #     i += 1
    # end

    # return true, savedValues
end
function (nfmu::CS_NeuralFMU{Vector{T}})(inputFct,
                                         t_step::Real,
                                         t_start::Real = nfmu.tspan[1],
                                         t_stop::Real = nfmu.tspan[end];
                                         reset::Bool = true,
                                         setup::Bool = true) where {T}

    if reset
        for fmu in nfmu.fmu
            c = fmu.components[end]
            if c.state == fmi2ComponentStateModelInitialized
                fmi2Terminate(c)
            end
            if c.state == fmi2ComponentStateModelSetableFMUstate
                fmi2Reset(c)
            end
        end
    end

    if setup
        for fmu in nfmu.fmu
            while fmi2SetupExperiment(fmu, t_start, t_stop) != 0
            end
            while fmi2EnterInitializationMode(fmu) != 0
            end
            while fmi2ExitInitializationMode(fmu) != 0
            end

            fmu.components[end].skipNextDoStep = true # skip first fim2DoStep-call
        end
    end

    ts = t_start:t_step:t_stop
    model_input = inputFct.(ts)
    valueStack = nfmu.model.(model_input)
    return valueStack

    # savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...})

    # ts = t_start:t_step:t_stop
    # i = 1
    # for t in ts
    #     DiffEqCallbacks.copyat_or_push!(savedValues.t, i, t)
    #     DiffEqCallbacks.copyat_or_push!(savedValues.saveval, i, (valueStack[i]...,), Val{false})
    #     i += 1
    # end

    # return true, savedValues
end

# adapting the Flux functions
function Flux.params(neuralfmu::ME_NeuralFMU)
    Flux.params(neuralfmu.neuralODE)
end

function Flux.params(neuralfmu::CS_NeuralFMU)
    Flux.params(neuralfmu.model)
end

# FMI version independent dosteps

"""
Wrapper. Call ```fmi2EvaluateME``` for more information.
"""
function fmiEvaluateME(str::fmi2Struct,
                     x::Array{<:Real},
                     t::Real = -1.0,
                     setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                     setValues::Array{<:Real} = zeros(Real, 0),
                     getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))
    fmi2EvaluateME(str, x, t,
                setValueReferences,
                setValues,
                getValueReferences)
end

"""
Wrapper. Call ```fmi2DoStepCS``` for more information.
"""
function fmiDoStepCS(str::fmi2Struct,
                     dt::Real,
                     setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                     setValues::Array{<:Real} = zeros(Real, 0),
                     getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))
    fmi2DoStepCS(str, dt, setValueReferences, setValues, getValueReferences)
end

"""
Wrapper. Call ```fmi2InputDoStepCSOutput``` for more information.
"""
function fmiInputDoStepCSOutput(str::fmi2Struct,
                                dt::Real,
                                u::Array{<:Real})
    fmi2InputDoStepCSOutput(str, dt, u)
end

# define neutral gradients (=feed-trough) for ccall-functions (ToDo: Remove)
function neutralGradient(c̄)
    tuple(c̄,)
end

@adjoint fmi2SetupExperiment(fmu, startTime, stopTime) = fmi2SetupExperiment(fmu, startTime, stopTime), c̄ -> neutralGradient(c̄)
@adjoint fmi2EnterInitializationMode(fmu) = fmi2EnterInitializationMode(fmu), c̄ -> neutralGradient(c̄)
@adjoint fmi2ExitInitializationMode(fmu) = fmi2ExitInitializationMode(fmu), c̄ -> neutralGradient(c̄)
@adjoint fmi2Reset(fmu) = fmi2Reset(fmu), c̄ -> neutralGradient(c̄)
@adjoint fmi2Terminate(fmu) = fmi2Terminate(fmu), c̄ -> neutralGradient(c̄)


text:
function funk(A,B)
    m = []
    s = 0
    a,b = size(A)
    for i in 1:a
        for j in 1:b
            if j > i
                push!(m,A[i,j])
            elseif i > j
                push!(m,B[i,j])
            end
            if B[i,j]>0 && rem(B[i,j],3)==0
                s+=B[i,j]
            end
        end
    end
    return m, s
end


text:
# Advent of Code 2018 - Day 9
# https://adventofcode.com/2018/day/9
#
# --- Day 9: Marble Mania ---
#
# You talk to the Elves while you wait for your navigation system
# to initialize. To pass the time, they introduce you to their
# favorite marble game.  The Elves play this game by taking turns
# arranging the marbles in a circle according to very particular
# rules. The marbles are numbered starting with 0 and increasing
# by 1 until every marble has a number.  First, the marble
# numbered 0 is placed in the circle. At this point, while it
# contains only a single marble, it is still a circle: the marble
# is both clockwise from itself and counter-clockwise from
# itself. This marble is designated the current marble.  Then,
# each Elf takes a turn placing the lowest-numbered remaining
# marble into the circle between the marbles that are 1 and 2
# marbles clockwise of the current marble. (When the circle is
# large enough, this means that there is one marble between the
# marble that was just placed and the current marble.) The marble
# that was just placed then becomes the current marble.  However,
# if the marble that is about to be placed has a number which is
# a multiple of 23, something entirely different happens. First,
# the current player keeps the marble they would have placed,
# adding it to their score. In addition, the marble 7 marbles
# counter-clockwise from the current marble is removed from the
# circle and also added to the current player's score. The marble
# located immediately clockwise of the marble that was removed
# becomes the new current marble.
#
# For example, suppose there are 9 players. After the marble with
# value 0 is placed in the middle, each player (shown in square
# brackets) takes a turn. The result of each of those turns would
# produce circles of marbles like this, where clockwise is to the
# right and the resulting current marble is in parentheses
#
#     [-] (0)
#     [1]  0 (1)
#     [2]  0 (2) 1
#     [3]  0  2  1 (3)
#     [4]  0 (4) 2  1  3
#     [5]  0  4  2 (5) 1  3
#     [6]  0  4  2  5  1 (6) 3
#     [7]  0  4  2  5  1  6  3 (7)
#     [8]  0 (8) 4  2  5  1  6  3  7
#     [9]  0  8  4 (9) 2  5  1  6  3  7
#     [1]  0  8  4  9  2(10) 5  1  6  3  7
#     [2]  0  8  4  9  2 10  5(11) 1  6  3  7
#     [3]  0  8  4  9  2 10  5 11  1(12) 6  3  7
#     [4]  0  8  4  9  2 10  5 11  1 12  6(13) 3  7
#     [5]  0  8  4  9  2 10  5 11  1 12  6 13  3(14) 7
#     [6]  0  8  4  9  2 10  5 11  1 12  6 13  3 14  7(15)
#     [7]  0(16) 8  4  9  2 10  5 11  1 12  6 13  3 14  7 15
#     [8]  0 16  8(17) 4  9  2 10  5 11  1 12  6 13  3 14  7 15
#     [9]  0 16  8 17  4(18) 9  2 10  5 11  1 12  6 13  3 14  7 15
#     [1]  0 16  8 17  4 18  9(19) 2 10  5 11  1 12  6 13  3 14  7 15
#     [2]  0 16  8 17  4 18  9 19  2(20)10  5 11  1 12  6 13  3 14  7 15
#     [3]  0 16  8 17  4 18  9 19  2 20 10(21) 5 11  1 12  6 13  3 14  7 15
#     [4]  0 16  8 17  4 18  9 19  2 20 10 21  5(22)11  1 12  6 13  3 14  7 15
#     [5]  0 16  8 17  4 18(19) 2 20 10 21  5 22 11  1 12  6 13  3 14  7 15
#     [6]  0 16  8 17  4 18 19  2(24)20 10 21  5 22 11  1 12  6 13  3 14  7 15
#     [7]  0 16  8 17  4 18 19  2 24 20(25)10 21  5 22 11  1 12  6 13  3 14  7 15
#     The goal is to be the player with the highest score after the last marble is used up. Assuming the example above ends after the marble numbered 25, the winning score is 23+9=32 (because player 5 kept marble 23 and removed marble 9, while no other player got any points in this very short example game).
#
#     Here are a few more examples:
#
#     10 players; last marble is worth 1618 points: high score is 8317
#     13 players; last marble is worth 7999 points: high score is 146373
#     17 players; last marble is worth 1104 points: high score is 2764
#     21 players; last marble is worth 6111 points: high score is 54718
#     30 players; last marble is worth 5807 points: high score is 37305

using CircularList, Dates, Test

function play(num_players, last_marble)
    leaderboard = fill(0, num_players)
    game = circularlist(0, capacity = last_marble)
    marble = 0
    while marble < last_marble
        marble += 1
        if marble % 23 == 0
            player = (marble % num_players) + 1
            leaderboard[player] += marble   # win current marble
            shift!(game, 7, :backward)
            leaderboard[player] += game.current.data  # win the other marble
            delete!(game)
            forward!(game)
        else
            forward!(game)
            insert!(game, marble)
        end
    end
    return (scores = leaderboard,
        highest = maximum(leaderboard),
        num_players = num_players,
        last_marble = last_marble
        )
end

# @test play(9, 25).highest == 32
# @test play(10, 1618).highest == 8317
# @time play(431, 70950)
@time play(431, 7095000)
# play(9, 25)


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

using SLiDE

READ_DIR = joinpath("data", "readfiles");

# include("io.jl")
# include("share_cfs.jl")
# include("share_gsp.jl")
# include("share_pce.jl")
# include("share_sgf.jl")
# include("share_utd.jl")

# ******************************************************************************************
#   READ SUPPLY/USE DATA
# ******************************************************************************************
println("  Reading supply/use data...")
DATA_DIR = joinpath("data", "input")
io = Dict(k => read_file(joinpath(DATA_DIR, string(k, ".csv"))) for k in [:supply, :use])

# ******************************************************************************************
#   READ SHARING DATA.
# ******************************************************************************************
files_share = write_yaml(READ_DIR, XLSXInput("generate_yaml.xlsx", "share", "B1:G150", "share"))
y = [read_file(files_share[ii]) for ii in 1:length(files_share)]
files_share = run_yaml(ensurearray(files_share))
shr = Dict(Symbol(y[ii]["PathOut"][end][1:end-4]) =>
    read_file(joinpath(y[ii]["PathOut"]...)) for ii in 1:length(y))

# Filter data and extrapolate values as appropriate.
shr = Dict(k => sort(filter_with(df, set; extrapolate = true)) for (k, df) in shr)
shr[:cfs][!,:value] .*= 1E3

shr_comp = Dict()
d = Dict(k => copy(shr[k]) for k in keys(shr))
d = merge(d, io)

# ******************************************************************************************
#   Do the things!
# ******************************************************************************************

#
io_keys = [:a0,:fd0,:fs0,:id0,:m0,:md0,:ms0,:ta0,:tm0,:va0,:x0,:y0,:ys0]
partition!(io, set)

io_check = Dict()
[benchmark!(io_check, k, bio, io) for k in io_keys]

# shr = merge(Dict(:va0 => io[:va0]), shr)
# share_pce!(shr)
# share_sgf!(shr)
# share_utd!(shr, set)
# share_region!(shr, set)
# share_labor!(shr, set)
# share_rpc!(shr, set)

# benchmark_keys = sort(setdiff(intersect(collect(keys(bshr)),collect(keys(shr))), [:wg,:hw,:notrd,:ng]))
# [benchmark!(shr_check, k, bshr, shr) for k in benchmark_keys];

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
function split2(str, on)
    i = findfirst(isequal(on), str)
    i === nothing && return str, ""
    return (SubString(str, firstindex(str), prevind(str, first(i))),
            SubString(str, nextind(str, last(i))))
end

"""
`SnoopCompile.read("compiledata.csv")` reads the log file produced by the compiler and returns the functions as a pair of arrays. The first array is the amount of time required to compile each function, the second is the corresponding function + types. The functions are sorted in order of increasing compilation time. (The time does not include the cost of nested compiles.)
"""
function read(filename)
    times = Vector{UInt64}()
    data = Vector{String}()
    toplevel = false # workaround for Julia#22538
    for line in eachline(filename)
        if toplevel
            if endswith(line, '"')
                toplevel = false
            end
            continue
        end
        time, str = split2(line, '\t')
        length(str) < 2 && continue
        (str[1] == '"' && str[end] == '"') || continue
        if startswith(str, """"<toplevel thunk> -> """)
            # consume lines until we find the terminating " character
            toplevel = true
            continue
        end
        tm = tryparse(UInt64, time)
        tm === nothing && continue
        push!(times, tm)
        push!(data, str[2:prevind(str, lastindex(str))])
    end
    # Save the most costly for last
    p = sortperm(times)
    return (permute!(times, p), permute!(data, p))
end

# pattern match on the known output of jl_static_show
extract_topmod(e::QuoteNode) = extract_topmod(e.value)
function extract_topmod(e)
    Meta.isexpr(e, :.) &&
        return extract_topmod(e.args[1])
    Meta.isexpr(e, :call) && length(e.args) == 3 && e.args[1] == :getfield &&
        return extract_topmod(e.args[2])
    Meta.isexpr(e, :call) && length(e.args) == 2 && e.args[1] == :typeof &&
        return extract_topmod(e.args[2])
    #Meta.isexpr(e, :call) && length(e.args) == 2 && e.args[1] == :Symbol &&
    #    return Symbol(e.args[2])
    # parametrized anonymous functions
    Meta.isexpr(e, :call) && e.args[1].args[1] == :getfield &&
        return extract_topmod(e.args[1].args[2].args[2])
    isa(e, Symbol) &&
        return e
    return :unknown
end

function parse_call(line; subst=Vector{Pair{String, String}}(), exclusions=String[])
    match(anonrex, line) === nothing || return false, line, :unknown, ""
    for (k, v) in subst
        line = replace(line, k=>v)
    end
    if any(b -> occursin(b, line), exclusions)
        println(line, " contains an excluded substring")
        return false, line, :unknown, ""
    end

    curly = ex = Meta.parse(line, raise=false)
    while Meta.isexpr(curly, :where)
        curly = curly.args[1]
    end
    if !Meta.isexpr(curly, :curly)
        @warn("failed parse of line: ", line)
        return false, line, :unknown, ""
    end
    func = curly.args[2]
    topmod = (func isa Expr ? extract_topmod(func) : :Main)

    check = Meta.isexpr(func, :call) && length(func.args) == 3 && func.args[1] == :getfield
    name = (check ? func.args[3].args[2] : "")

    # make some substitutions to try to form a leaf types tuple
    changed = false
    for i in 3:length(curly.args)
        e = curly.args[i]
        if Meta.isexpr(e, :where) && length(e.args) == 2 && Meta.isexpr(e.args[1], :curly) && length(e.args[1].args) == 3 && e.args[1].args[1] === :Vararg
            # Unwrap a Varargs argument, make it a single instance of the allowed Varargs type.
            e = e.args[1].args[2]
            curly.args[i] = e
            changed = true
        end
        if e === :Function
            # `Function` is not an allowable type for `precompile()`
            e = :(typeof(identity))
        elseif e === :Any
            # `Any` is an abstract type, and can't be precompiled for. Replace it with something that can.
            e = :Int
        elseif e == :(Type{T} where T)
            # Type is an abstract type; replace it with a concrete type.
            e = :(Type{Int})
        else
            continue
        end
        curly.args[i] = e
        changed = true
    end
    # In Julia 0.6, functions with symbolic names like `Base.:(+)` don't output correctly.
    # So if we didn't change the arg list above, use the original string.
    if changed
        line = string(ex)
    end
    return true, line, topmod, name
end

"""
`pc = parcel(calls; subst=[], exclusions=[])` assigns each compile statement to the module that owns the function. Perform string substitution via `subst=["Module1"=>"Module2"]`, and omit functions in particular modules with `exclusions=["Module3"]`. On output, `pc[:Module2]` contains all the precompiles assigned to `Module2`.

Use `SnoopCompile.write(prefix, pc)` to generate a series of files in directory `prefix`, one file per module.
"""
function parcel(calls::AbstractVector{String};
    subst=Vector{Pair{String, String}}(),
    exclusions=String[],
    remove_exclusions::Bool = true,
    check_eval::Bool = false,
    blacklist=nothing,                  # deprecated keyword
    remove_blacklist=nothing)           # deprecated keyword

    if blacklist !== nothing
        Base.depwarn("`blacklist` is deprecated, please use `exclusions` to pass a list of excluded names", :parcel)
        append!(exclusions, blacklist)
    end
    if remove_blacklist !== nothing
        Base.depwarn("`remove_blacklist` is deprecated, please use `remove_exclusions` to pass a list of excluded names", :parcel)
        remove_exclusions = remove_blacklist
    end
    pc = Dict{Symbol, Vector{String}}()

    sym_module = Dict{Symbol, Module}() # 1-1 association between modules and module name

    for c in calls
        local keep, pcstring, topmod
        keep, pcstring, topmod, name = parse_call(c, subst=subst, exclusions=exclusions)
        keep || continue
        # Add to the appropriate dictionary
        if !haskey(pc, topmod)
            pc[topmod] = String[]
        end
        prefix = (isempty(name) ? "" : "isdefined($topmod, Symbol(\"$name\")) && ")
        push!(pc[topmod], prefix * "precompile($pcstring)")

        # code from parcel_snoopi
        # TODO moduleof a symobl!
        # eval(topmod) # throws error for Test
        # sym_module[topmod] = moduleof(topmod)
    end

    # TODO
    # for mod in keys(pc)
    #     # check_eval remover
    #     if check_eval
    #         pc[mod] = remove_if_not_eval!(pc[mod], sym_module[mod])
    #     end
    # end
    return pc
end

"""
`pc = format_userimg(calls; subst=[], exclusions=[])` generates precompile directives intended for your base/userimg.jl script. Use `SnoopCompile.write(filename, pc)` to create a file that you can `include` into `userimg.jl`.
"""
function format_userimg(calls; subst=Vector{Pair{String, String}}(), exclusions=String[])
    pc = Vector{String}()
    for c in calls
        keep, pcstring, topmod, name = parse_call(c, subst=subst, exclusions=exclusions)
        keep || continue
        prefix = (isempty(name) ? "" : "isdefined($topmod, Symbol(\"$name\")) && ")
        push!(pc, prefix * "precompile($pcstring)")
    end
    return pc
end

# TODO Make this work with snoopc:

"""
    remove_if_not_eval!(pcstatements, modul::Module)

Removes everything statement in `pcstatements` can't be `eval`ed in `modul`.

# Example

```jldoctest; setup=:(using SnoopCompile), filter=r":\\d\\d\\d"
julia> pcstatements = ["precompile(sum, (Vector{Int},))", "precompile(sum, (CustomVector{Int},))"];

julia> SnoopCompile.remove_if_not_eval!(pcstatements, Base)
┌ Warning: Faulty precompile statement: precompile(sum, (CustomVector{Int},))
│   exception = UndefVarError: CustomVector not defined
└ @ Base precompile_Base.jl:375
1-element Vector{String}:
 "precompile(sum, (Vector{Int},))"
```
"""
function remove_if_not_eval!(pcstatements, modul::Module)
    todelete = Set{eltype(pcstatements)}()
    for line in pcstatements
        try
            if modul === Core
                #https://github.com/timholy/SnoopCompile.jl/issues/76
                Core.eval(Main, Meta.parse(line))
            else
                Core.eval(modul, Meta.parse(line))
            end
        catch e
            @warn("Faulty precompile statement: $line", exception = e, _module = modul, _file = "precompile_$modul.jl")
            push!(todelete, line)
        end
    end
    return setdiff!(pcstatements, todelete)
end


text:
export BasicDQNLearner

using Flux: gradient, params
using Zygote: ignore
using Setfield: @set

import Functors

"""
    BasicDQNLearner(;kwargs...)

See paper: [Playing Atari with Deep Reinforcement Learning](https://www.cs.toronto.edu/~vmnih/docs/dqn.pdf)

This is the very basic implementation of DQN. Compared to the traditional Q
learning, the only difference is that, in the optimising step it uses a batch of
transitions sampled from an experience buffer instead of current transition. And
a neural network is used to estimate the Q-value.  You can start from this
implementation to understand how everything is organized and how to write your
own customized algorithm.

# Keyword Arguments

- `approximator`::[`Approximator`](@ref): used to get Q-values of a state.
- `loss_func=huber_loss`: the loss function to use.
- `γ::Float32=0.99f0`: discount rate.
"""
Base.@kwdef mutable struct BasicDQNLearner{Q} <: AbstractLearner
    approximator::Q
    loss_func::Any = huber_loss
    γ::Float32 = 0.99f0
    # for debugging
    loss::Float32 = 0.0f0
end

Functors.functor(x::BasicDQNLearner) = (Q=x.approximator,), y -> @set x.approximator = y.Q

(L::BasicDQNLearner)(s::AbstractArray) = L.approximator(s)

function RLCore.optimise!(
    learner::BasicDQNLearner,
    batch::NamedTuple{SS′ART},
)

    Q = learner.approximator
    γ = learner.γ
    loss_func = learner.loss_func

    s, s′, a, r, t = send_to_device(device(Q), batch)
    a = CartesianIndex.(a, 1:length(a))

    gs = gradient(params(Q)) do
        q = Q(s)[a]
        q′ = vec(maximum(Q(s′); dims=1))
        G = @. r + γ * (1 - t) * q′
        loss = loss_func(G, q)
        ignore() do
            learner.loss = loss
        end
        loss
    end

    optimise!(Q, gs)
end


text:
# ---
# title: 1339. Maximum Product of Splitted Binary Tree
# id: problem1339
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Tree, Depth-first Search
# link: <https://leetcode.com/problems/maximum-product-of-splitted-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree `root`. Split the binary tree into two subtrees by
# removing 1 edge such that the product of the sums of the subtrees are
# maximized.
#
# Since the answer may be too large, return it modulo 10^9 + 7.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/21/sample_1_1699.png)**
#
#
#
#     Input: root = [1,2,3,4,5,6]
#     Output: 110
#     Explanation: Remove the red edge and get 2 binary trees with sum 11 and 10. Their product is 110 (11*10)
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/21/sample_2_1699.png)
#
#
#
#     Input: root = [1,null,2,3,4,null,null,5,6]
#     Output: 90
#     Explanation:  Remove the red edge and get 2 binary trees with sum 15 and 6.Their product is 90 (15*6)
#
#
# **Example 3:**
#
#
#
#     Input: root = [2,3,9,10,7,8,6,5,4,11,1]
#     Output: 1025
#
#
# **Example 4:**
#
#
#
#     Input: root = [1,1]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * Each tree has at most `50000` nodes and at least `2` nodes.
#   * Each node's value is between `[1, 10000]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "X-2-5 scenario" begin

@time for X in 6:7
    filepath = "./data/quick_adjacency_decomposition/"
    filename = "$X-2-5"

    println("$filename scenario")

    PM = LocalSignaling(X,5,2)
    BG_seed = BellGame(cat([1 0 0 0 0;1 0 0 0 0;1 0 0 0 0;0 1 0 0 0;0 0 1 0 0],zeros(Int64,(5,X-5)), dims=2),2)

    vertices = LocalPolytope.vertices(PM)
    canonical_facets = LocalPolytope.adjacency_decomposition(vertices, BG_seed, PM, max_vertices = 100)
    games = collect(keys(canonical_facets))

    pretty_print_txt(games, filepath*filename)
    write_ieq(filepath*"ieq/"*filename, convert(IEQ, games))
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "IOAPI"
version = v"3.2.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/cjcoats/ioapi-3.2.git", "ef5d5f4e112c249b593b19426421f25d79ae094b")
]

# Bash recipe for building across all platforms
script = raw"""
apk add tcsh # Build script is in csh

cd $WORKSPACE/srcdir/ioapi-3.2/ioapi/
install_license ../LICENSE
export HOME=${WORKSPACE}/srcdir
export BIN=tmp
export BINDIR=${HOME}/ioapi-3.2/${BIN}
mkdir $BINDIR

sed -i 's/-ffast-math /-fPIC /g' Makeinclude* # Get rid of forbidden -ffast-math flag
sed -i 's/-DNEED_ARGS=1/-DNEED_ARGS=1 -DIOAPI_NCF4=1/g' Makeinclude* # Specify NetCDF v4

if [[ "${nbits}" == 32 ]]; then
    sed -i 's/-m64/-m32/g' Makeinclude* # Specify 32-bit build
elif [[ "${target}" == *aarch64* ]]; then
    sed -i 's/-m64//g' Makeinclude* # Get rid of x86 flag
fi

if [[ "${target}" == *mingw* ]]; then
    # Add missing stdint header.
    sed -i 's/<stdio.h>/<stdio.h>\n\#include <stdint.h>/g' bufint3.c

    # This sys/wait.h header is not available and doesn't seem to be used.
    sed -i 's/\#include <sys\/wait.h>//g' systemf.c
fi

cp Makefile.nocpl Makefile
cp Makeinclude.Linux2_x86_64gfort Makeinclude.$BIN

make # Parallel make (-j > 1) does not work

cp *.h ${includedir} # C header files
cp fixed_src/* ${includedir} # FORTRAN .EXT (include) files

# Convert static library to dynamic library
gfortran -shared -fPIC -fopenmp -o ${libdir}/libioapi.${dlext} -L${libdir} -Wl,$(flagon --whole-archive) ${BINDIR}/libioapi.a -Wl,$(flagon --no-whole-archive | cut -d' ' -f1) -lnetcdf -Wl,$(flagon --no-whole-archive | cut -d' ' -f1) -lnetcdff
rm ${BINDIR}/libioapi.a

cd ../m3tools/
cp Makefile.nocpl Makefile
make

cd $BINDIR

rm *.o
mv *.mod ${includedir} # Move FORTRAN mod files to include dir, they are used by some dependencies
mv * $bindir # Move executables to bindir
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("x86_64", "macos"; ),
    Platform("aarch64", "linux"; libc = "glibc"),
    #Platform("x86_64", "windows"),
    #Platform("i686", "windows")
]
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("dayagg", :dayagg),
    ExecutableProduct("wrftom3", :wrftom3),
    ExecutableProduct("kfxtract", :kfxtract),
    ExecutableProduct("m3combo", :m3combo),
    ExecutableProduct("m3wndw", :m3wndw),
    ExecutableProduct("wrfgriddesc", :wrfgriddesc),
    ExecutableProduct("m3probe", :m3probe),
    ExecutableProduct("m3agmax", :m3agmax),
    ExecutableProduct("m3totxt", :m3totxt),
    ExecutableProduct("timediff", :timediff),
    ExecutableProduct("factor", :factor),
    ExecutableProduct("presz", :presz),
    ExecutableProduct("presterp", :presterp),
    ExecutableProduct("gregdate", :gregdate),
    ExecutableProduct("m3stat", :m3stat),
    ExecutableProduct("m3pair", :m3pair),
    ExecutableProduct("projtool", :projtool),
    LibraryProduct("libioapi", :libioapi),
    ExecutableProduct("datshift", :datshift),
    ExecutableProduct("juldate", :juldate),
    ExecutableProduct("m3cple", :m3cple),
    ExecutableProduct("wndwdesc", :wndwdesc),
    ExecutableProduct("mpasstat", :mpasstat),
    ExecutableProduct("vertintegral", :vertintegral),
    ExecutableProduct("jul2greg", :jul2greg),
    ExecutableProduct("m3fake", :m3fake),
    ExecutableProduct("m3merge", :m3merge),
    ExecutableProduct("greg2jul", :greg2jul),
    ExecutableProduct("latlon", :latlon),
    ExecutableProduct("mpasdiff", :mpasdiff),
    ExecutableProduct("airs2m3", :airs2m3),
    ExecutableProduct("m3edhdr", :m3edhdr),
    ExecutableProduct("m3diff", :m3diff),
    ExecutableProduct("randomstat", :randomstat),
    ExecutableProduct("mtxbuild", :mtxbuild),
    ExecutableProduct("mtxcple", :mtxcple),
    ExecutableProduct("insertgrid", :insertgrid),
    ExecutableProduct("m3xtract", :m3xtract),
    ExecutableProduct("vertot", :vertot),
    ExecutableProduct("juldiff", :juldiff),
    ExecutableProduct("camxtom3", :camxtom3),
    ExecutableProduct("mpastom3", :mpastom3),
    ExecutableProduct("vertimeproc", :vertimeproc),
    ExecutableProduct("findwndw", :findwndw),
    ExecutableProduct("m3hdr", :m3hdr),
    ExecutableProduct("m3tproc", :m3tproc),
    ExecutableProduct("mtxblend", :mtxblend),
    ExecutableProduct("selmrg2d", :selmrg2d),
    ExecutableProduct("m3tshift", :m3tshift),
    ExecutableProduct("timeshift", :timeshift),
    ExecutableProduct("julshift", :julshift),
    ExecutableProduct("bcwndw", :bcwndw),
    ExecutableProduct("gridprobe", :gridprobe),
    ExecutableProduct("m3agmask", :m3agmask),
    ExecutableProduct("m3interp", :m3interp),
    ExecutableProduct("m3mask", :m3mask),
    ExecutableProduct("mtxcalc", :mtxcalc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="NetCDF_jll", uuid="7243133f-43d8-5620-bbf4-c2c921802cf3"))
    Dependency(PackageSpec(name="NetCDFF_jll", uuid="78e728a9-57fe-5d11-897c-5014b89e5f84"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using BinDeps
@BinDeps.setup

version = "3360"
url = "ftp://heasarc.gsfc.nasa.gov/software/fitsio/c/cfitsio$version.tar.gz"

libcfitsio = library_dependency("libcfitsio")
provides(Sources, URI(url), libcfitsio, unpacked_dir="cfitsio")
depsdir = BinDeps.depsdir(libcfitsio)
srcdir = joinpath(depsdir, "src", "cfitsio")
prefix = joinpath(depsdir, "usr")
@unix_only libfilename = "libcfitsio.so"
@osx_only libfilename = "libcfitsio.dylib"
provides(BuildProcess,
         (@build_steps begin
            GetSources(libcfitsio)
            @build_steps begin
                ChangeDirectory(srcdir)
                FileRule(joinpath(prefix,"lib",libfilename),
                         @build_steps begin
                         `./configure --prefix=$prefix`
                         `make shared install`
                         end)
            end
         end),
         libcfitsio)

@BinDeps.install [:libcfitsio => :libcfitsio]


text:
using LinearAlgebra
using Statistics
module LearningHorse

    export LinearRegression
    module LinearRegression
        module SGD #this is Slope Gradient Descent
            include("./LinearRegression/SGD.jl")
        end
        module MR #this is Multiple Regression
            include("./LinearRegression/MR.jl")
        end
        module RR #this is Ridge Regression
            include("./LinearRegression/RR.jl")
        end
        module LR #this is Lasso Regression
            include("./LinearRegression/LR.jl")
        end
        module BFM #this use Basis Function Model for regression
            include("./LinearRegression/BFM.jl")
        end
    end

    export Classification
    module Classification #this is Classification
        include("./Classification/Other.jl") #This file defined LE and CTOH
        module MS #This is Multiclass softmax
            include("./Classification/MS.jl")
        end
        module OVR #this is One-vs-Rest
            include("./Classification/OVR.jl")
        end
    end
    export Tree
    module Tree #This is Decision Tree
        module DT #This is Decision Tree
            include("./Tree/DT.jl")
        end
        include("./Tree/MV.jl")
    end

    export Ensemble
    module Ensemble
        module RF
            include("./Ensemble/RF.jl")
        end
    end

    export NeuralNetwork
    module NeuralNetwork
        include("./NeuralNetwork/NN.jl")
    end
    export Preprocessing
    module Preprocessing
        module SS() #this is Standard Scaler
            include("./Preprocessing/SS.jl")
        end
        module MM() #this is MinMax Scaler
            include("./Preprocessing/MM.jl")
        end
        module RS() #this is Robust Scaler
            include("./Preprocessing/RS.jl")
        end
        module Data() # This is the dataloader and datasplitter.
            include("./Preprocessing/Data.jl")
        end
    end
    export LossFunction
    module LossFunction #Loss Fuunctions
        include("./LossFunction/LossFunctions.jl")
    end

end


text:
# Module written by Ed Scheinerman, ers@jhu.edu
# distributed under terms of the MIT license

module SimpleGraphs
using DataStructures
using Iterators

"""
`AbstractSimpleGraph` is a parent class for `SimpleGraph` and `SimpleDigraph`.
"""
abstract AbstractSimpleGraph
export AbstractSimpleGraph

include("simple_core.jl")
include("simple_ops.jl")
include("simple_constructors.jl")
include("platonic.jl")
include("simple_connect.jl")
include("simple_matrices.jl")
include("disjoint_sets_helper.jl")
include("simple_converters.jl")
include("simple_coloring.jl")
include("simple_euler.jl")
include("simple_ham.jl")
include("bisect.jl")
include("trans_orient.jl")

include("d_simple_core.jl")
include("d_simple_constructors.jl")
include("d_simple_matrices.jl")

end # module SimpleGraphs


text:
#!/usr/bin/julia

using Pkg
pkg"activate ."

using Makie

scene = Scene()

x = 0:0.1:1
y = sin.(x)

scatter!(x, y, markersize=0.05)
lines!(x, y, color=:blue)

sleep(10)

exit()


text:
function load_data()
    obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "test_data.uvfits"))
    obs.add_scans()
    obsavg = obs.avg_coherent(0.0, scan_avg=true)

    m = ehtim.model.Model()
    m = m.add_gauss(1.0, μas2rad(40.0), μas2rad(20.0), π/3, 0.0, 0.0)
    m = m.add_gauss(0.5, μas2rad(20.0), μas2rad(10.0), π/6, μas2rad(30.0), μas2rad(30.0))
    obsm = m.observe_same_nonoise(obsavg)

    vis = extract_vis(obsm)
    amp = extract_amp(obsm)
    lcamp = extract_lcamp(obsm)
    cphase = extract_cphase(obsm, cut_trivial=true)

    return m, vis, amp, lcamp, cphase
end

function test_model(θ)
    m1 = θ.f1*rotated(stretched(Gaussian(), θ.σ1*θ.τ1, θ.σ1), θ.ξ1)
    m2 = θ.f2*rotated(stretched(Gaussian(), θ.σ2*θ.τ2, θ.σ2), θ.ξ2)
    return m1 + shifted(m2, θ.x, θ.y)
end

function test_prior()
    return (f1=Uniform(0.8, 1.2),
             σ1 = Uniform(μas2rad(1.0), μas2rad(40.0)),
             τ1 = Uniform(0.35, 0.65),
             ξ1 = Uniform(-π/2, π/2),
             f2 = Uniform(0.3, 0.7),
             σ2 = Uniform(μas2rad(1.0), μas2rad(40.0)),
             τ2 = Uniform(0.35, 0.65),
             ξ2 = Uniform(-π/2, π/2),
             x = Uniform(-μas2rad(40.0), μas2rad(40.0)),
             y = Uniform(-μas2rad(40.0), μas2rad(40.0))
            )
end


text:
using ScikitLearn
using PyCall

const naive_bayes = PyNULL()
const model_selection = PyNULL()
function __init__()
    copy!(naive_bayes, PyCall.pyimport_conda("sklearn.naive_bayes", "scikit-learn"))
    copy!(model_selection, PyCall.pyimport_conda("sklearn.model_selection", "scikit-learn"))
end

"""
function ecrelieff(data::Array{<:Real,2}, target::Array{<:Integer,1}, m, k,
                    dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                    mode::String="k_nearest", sig::Real=1.0, f_type::String="continuous")::Array{Int64,1}

Compute feature rankings using Evaporative Cooling ReliefF algorithm.

---
# Reference:
- B.A. McKinney, D.M. Reif, B.C. White, J.E. Crowe, Jr., J.H. Moore.
Evaporative cooling feature selection for genotypic data involving interactions.
"""
function ecrelieff(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1, k::Integer=10,
                   dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                   mode::String="k_nearest", sig::Real=1.0, f_type::String="continuous")::Array{Int64,1}

    # Compute mu values.
    mu_vals = get_mu_vals(data, target)

	# Set initial tinfo value.
	tinfo = 1

	# Initialize vector of feature ranks.
	# rank = Array{Int64}(undef, length(weights))
	rank = Array{Int64}(undef, size(data, 2))

	# Initialize vector of indices in original weights array.
	# index_vector = collect(1:length(weights))
	index_vector = collect(1:size(data, 2))

	# Initialize initial rank value.
	# rank_value_next = length(weights)
	rank_value_next = size(data, 2)

	# Initialize variable that holds current best tinfo value.
	best_tinfo = tinfo

	# While there are unranked features, perform evaporation.
	@inbounds while length(index_vector) > 1

        # Compute ReliefF weights.
        weights = relieff(data, target, m, k, dist_func, mode=mode, sig=sig, f_type=f_type)

        # Get maximal ReliefF weight and compute epsilon values.
        max_weight = maximum(weights)
        eps = (max_weight .- weights)./max_weight

		### Grid search for best temperature from previous step ###

		# Initialize variables that hold current maximal CV score
		# and index of removed feature.
		max_cv_val = 0.0
		idx_removed_feature = -1

		# Perform grid search for best value of tinfo.
		@inbounds for tinfo_nxt = tinfo-0.3:0.1:tinfo+0.3

			# Compute weights for next value of tinfo.
			weights_nxt = eps .- tinfo_nxt.*mu_vals

			# Rank features.
			enumerated_weights = [weights_nxt collect(1:length(weights_nxt))]'
			rank_weights = zeros(Int64, length(weights_nxt))
			s = enumerated_weights[:, sortperm(enumerated_weights[1,:], rev=false)]
			rank_weights[convert.(Int64, s[2, :])] = 1:length(weights_nxt)

			# Remove lowest ranked feature.
			msk_rm = rank_weights .!= length(weights_nxt)

			# Perform 5-fold cross validation.
			data_filt = data[:, msk_rm]

			# Compare to current maximal CV value.
			cv_val_nxt = Statistics.mean(model_selection.cross_val_score(naive_bayes.GaussianNB(), data_filt, target; cv=5))

			# If CV value greater than current maximal, save current
			# CV score, tinfo and feature weights.
			if cv_val_nxt > max_cv_val
				# Set current tinfo value as best.
				best_tinfo = tinfo_nxt
				# Get index of removed feature.
				idx_removed_feature = findall(.!msk_rm)[1]
			end
		end

		# Remove evaporated feature from data matrix,
		# and data at corresponding index from eps vector and
		# mu values vector.
		data = data[:, 1:end .!= idx_removed_feature]
		eps = eps[1:end .!= idx_removed_feature]
		mu_vals = mu_vals[1:end .!= idx_removed_feature]

		# Get index of evaporated feature in original data matrix.
		rank_index_next = index_vector[idx_removed_feature]

		# Delete removed feature from vector of indices.
		index_vector = deleteat!(index_vector, idx_removed_feature)

		# Assign rank to removed feature.
		rank[rank_index_next] = rank_value_next

		# Decrease next rank value.
		rank_value_next -= 1

	end

	# Assign final rank.
	rank[index_vector[1]] = rank_value_next

	# Return vector of ranks.
	return rank

end

"""
    function joint_entropy_pair(distribution1::Array{<:Real,2}, distribution2::Array{<:Real,2})::Float64

Compute joint entropy of two distributions (auxiliary function).
"""
function joint_entropy_pair(distribution1::Array{<:Real,1}, distribution2::Array{<:Real,1})::Float64

    stacked = hcat(distribution1, distribution2)
    unique_stacked = unique(stacked, dims=1)
    counts_pairs = Array{Int64}(undef, size(unique_stacked, 1))
    @inbounds for idx = 1:size(unique_stacked, 1)
        counts_pairs[idx] = sum(all(stacked' .== reshape(unique_stacked[idx, :], length(unique_stacked[idx, :]), 1), dims=1))
    end

    p_pairs = counts_pairs/length(distribution1)
    return sum(p_pairs.*log.(p_pairs))
end

"""
    scaled_mutual_information(distribution1::Array{<:Real,2}, distribution2::Array{<:Real,2})::Float64

Compute scaled mutual information of two distributions (auxiliary function).
"""
function scaled_mutual_information(distribution1::Array{<:Real,1}, distribution2::Array{<:Real,1})::Float64
    entropy_dist1 = StatsBase.entropy(distribution1)
    entropy_dist2 = StatsBase.entropy(distribution2)
    return (entropy_dist1 + entropy_dist2 - joint_entropy_pair(distribution1, distribution2))/entropy_dist1
end

"""
    mu_vals(data::Array{<:Real,2}, target::Array{<:Integer,1})

Compute mu values from data and target (auxiliary function).
"""
function get_mu_vals(data::Array{<:Real,2}, target::Array{<:Integer,1})
    mu_vals = Array{Float64}(undef, size(data, 2))
    @inbounds for idx = 1:size(data, 2)
        mu_vals[idx] = scaled_mutual_information(Array{Real,1}(data[:, idx]), target)
    end
    return mu_vals
end


text:
using MINLPLib

LIBS = ["bcp","global", "ibm", "inf", "minlp",
        "minlp2", "morg", "mpec", "mult3", "mult4",
        "PODLib", "poly", "prince", "qcqp", "qcqp2",
        "qcqp3"]
LIBS = ["PODLib"]
function collect_dimension(libname::AbstractString, dat::Dict)

    names = MINLPLib.fetch_names(libname)
    for i in names
        meta = MINLPLib.fetch_meta(libname, i)
        if !isempty(meta)
            haskey(meta, "LOAD") ? dat[i]["lt"] = meta["LOAD"] : dat[i]["lt"] = -1.0
            dat[i]["sense"] = meta["OBJSENSE"]
            dat[i]["nvars"] = meta["NVARS"]
            dat[i]["nbinvars"] = meta["NBINVARS"]
            dat[i]["nintvars"] = meta["NINTVARS"]
            dat[i]["ncons"] = meta["NCONS"]
            dat[i]["nlincons"] = meta["NLINCONS"]
            dat[i]["nnlcons"] = meta["NNLCONS"]
            dat[i]["nothercons"] = meta["NCONS"] - meta["NLINCONS"] - meta["NNLCONS"]
            nlo = ""
            if haskey(meta, "NLOPERATOR")
                for j in meta["NLOPERATOR"]
                    nlo = "$(nlo) $(j)"
                end
            end
            dat[i]["nloperators"] = nlo
        end
    end

    return dat
end

function generate_lib_page(libname::AbstractString)
    isdir("temp") || mkdir("temp")
    f = open(joinpath("temp",libname), "w")
    names = MINLPLib.fetch_names(libname)
    N = length(names)
    DAT = Dict(i => Dict() for i in names)
    DAT = collect_dimension(libname, DAT)
    write(f, "$(N)\n\n")
    write(f, "| NAME | LT | NLOPERATOR | SENSE | VAR | BINVAR | INTVAR | CON | LINCON | NLCONS | OTHERCONS |\n")
    write(f, "|------|----|------------|-------|-----|--------|--------|-----|--------|--------|-----------|\n")
    err = 0
    for i in names
        try
            write(f,"| $i | $(round(DAT[i]["lt"],2)) | $(DAT[i]["nloperators"]) | $(DAT[i]["sense"]) | $(DAT[i]["nvars"]) | $(DAT[i]["nbinvars"]) | $(DAT[i]["nintvars"]) | $(DAT[i]["ncons"]) | $(DAT[i]["nlincons"]) | $(DAT[i]["nnlcons"]) | $(DAT[i]["nothercons"]) |\n")
        catch e
            err += 1
            write(f,"| $i | - | - | - | - | - | - | - | - | - |\n")
        end
    end
    write(f, "\nERROR $(err)\n")
    close(f)
end

for l in LIBS
    generate_lib_page(l)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# compiler driver and main interface

export cufunction

# (::CompilerContext)
const compile_hook = Ref{Union{Nothing,Function}}(nothing)

# Main entry point for compiling a Julia function + argtypes to a callable CUDA function
function cufunction(dev::CuDevice, @nospecialize(f), @nospecialize(tt); kwargs...)
    CUDAnative.configured || error("CUDAnative.jl has not been configured; cannot JIT code.")
    isa(f, Core.Function) || throw(ArgumentError("Kernel argument to `cufunction` should be a function."))

    ctx = CompilerContext(f, tt, supported_capability(dev), true; kwargs...)

    if compile_hook[] != nothing
        global globalUnique
        previous_globalUnique = globalUnique
        compile_hook[](ctx)
        globalUnique = previous_globalUnique
    end

    (module_asm, module_entry) = compile_function(ctx)

    # enable debug options based on Julia's debug setting
    jit_options = Dict{CUDAdrv.CUjit_option,Any}()
    if Base.JLOptions().debug_level == 1
        jit_options[CUDAdrv.GENERATE_LINE_INFO] = true
    elseif Base.JLOptions().debug_level >= 2
        jit_options[CUDAdrv.GENERATE_DEBUG_INFO] = true
    end
    cuda_mod = CuModule(module_asm, jit_options)
    cuda_fun = CuFunction(cuda_mod, module_entry)

    @debug begin
        attr = attributes(cuda_fun)
        bin_ver = VersionNumber(divrem(attr[CUDAdrv.FUNC_ATTRIBUTE_BINARY_VERSION],10)...)
        ptx_ver = VersionNumber(divrem(attr[CUDAdrv.FUNC_ATTRIBUTE_PTX_VERSION],10)...)
        regs = attr[CUDAdrv.FUNC_ATTRIBUTE_NUM_REGS]
        local_mem = attr[CUDAdrv.FUNC_ATTRIBUTE_LOCAL_SIZE_BYTES]
        shared_mem = attr[CUDAdrv.FUNC_ATTRIBUTE_SHARED_SIZE_BYTES]
        constant_mem = attr[CUDAdrv.FUNC_ATTRIBUTE_CONST_SIZE_BYTES]
        """Compiled $f to PTX $ptx_ver for SM $bin_ver using $regs registers.
           Memory usage: $local_mem B local, $shared_mem B shared, $constant_mem B constant"""
    end

    return cuda_fun, cuda_mod
end

# Compile a function to PTX, returning the assembly and an entry point.
# Not to be used directly, see `cufunction` instead.
# FIXME: this pipeline should be partially reusable from eg. code_llvm
#        also, does the kernel argument belong in the compiler context?
function compile_function(ctx::CompilerContext; strip_ir_metadata::Bool=false)
    ## high-level code generation (Julia AST)

    @debug "(Re)compiling function" ctx

    check_method(ctx)

    ## low-level code generation (LLVM IR)

    mod, entry = irgen(ctx)

    # link libdevice, if it is necessary
    libdevice = load_libdevice(ctx)
    for f in functions(mod)
        if isdeclaration(f) && intrinsic_id(f)==0 && haskey(functions(libdevice), LLVM.name(f))
            libdevice_copy = LLVM.Module(libdevice)
            link_libdevice!(ctx, mod, libdevice_copy)
            break
        end
    end

    # optimize the IR
    entry = optimize!(ctx, mod, entry)

    check_invocation(ctx, entry)

    # check generated IR
    check_ir(ctx, mod)
    verify(mod)

    if strip_ir_metadata
        strip_debuginfo!(mod)
    end

    ## machine code generation (PTX assembly)

    module_asm = mcgen(ctx, mod, entry)

    return module_asm, LLVM.name(entry)
end