text:
module ConstructiveSolidGeometry

include("types.jl")
include("util.jl")
include("ray_trace.jl")
include("region_operators.jl")
include("cell_utils.jl")
include("plotting.jl")

export Coord
export Ray
export Surface
export Plane
export Cone
export Sphere
export InfCylinder
export Box
export Region
export Cell
export Geometry
export +,-,*,^,|,~
export reflect
export generate_random_ray
export raytrace
export find_intersection
export halfspace
export is_in_cell
export find_cell_id
export plot_geometry_2D
export plot_cell_2D
export dot
export magnitude
export unitize
export cross

using Plots

end # module

ext: jl


text:
# Decomposes a scene into smaller scenes made up of pairwise agents
function decompose_scene(scene::Scene, egoid::Int)
    scenes = Dict{Int, Scene}()
    ego = set_veh_id(get_by_id(scene, egoid), 2)
    for i=1:egoid-1
        if has_veh(i, scene)
            veh = set_veh_id(get_by_id(scene, i), 1)
            scenes[i] = Scene([veh, ego])
        end
    end
    scenes
end

# Gets the observation indices of the desired agents (by id)
function decompose_indices(mdp::AdversarialDrivingMDP; ids)
    i, indices = 1, Int[]
    for a in agents(mdp)
        id(a) in ids && push!(indices, (i:i+a.entity_dim - 1)...)
        i += a.entity_dim
    end
    indices
end

## T-Intersection helper functions
ez_Tint_vehicle(;id::Int64, s::Float64, v::Float64, lane::Int64) = BlinkerVehicle(roadway = Tint_roadway,
                                              lane=lane, s=s, v = v,
                                              id = id, goals = Tint_goals[lane],
                                              blinker = Tint_should_blink[lane])

function get_Tint_vehicle(;id::Int64, s::Float64, v::Float64, lane::Int64)
    (rng::AbstractRNG = Random.GLOBAL_RNG) -> ez_Tint_vehicle(id=id, s=s, v=v, lane=lane)
end

function get_rand_Tint_vehicle(;id::Int64, s_dist, v_dist, lane_dist)
    (rng::AbstractRNG = Random.GLOBAL_RNG) -> ez_Tint_vehicle(id=id, s=rand(rng, s_dist), v=rand(rng, v_dist), lane=rand(rng, lane_dist))
end

up_left(;id::Int64, s::Float64 = 40., v::Float64 = 20.) = get_Tint_vehicle(id = id, s = s, v = v, lane = 5)
left_straight(;id::Int64, s::Float64 = 15., v::Float64 = 19.) = get_Tint_vehicle(id = id, s = s, v = v, lane = 2)
left_turnright(;id::Int64, s::Float64 = 15., v::Float64 = 19.) = get_Tint_vehicle(id = id, s = s, v = v, lane = 1)
right_straight(;id::Int64, s::Float64 = 30., v::Float64 = 20.) = get_Tint_vehicle(id = id, s = s, v = v, lane = 3)
right_turnleft(;id::Int64, s::Float64 = 40., v::Float64 = 14.) = get_Tint_vehicle(id = id, s = s, v = v, lane = 4)

rand_up_left(;id::Int64, s_dist::Distribution, v_dist::Distribution) = get_rand_Tint_vehicle(id=id, s_dist=s_dist, v_dist=v_dist, lane_dist = [5])
rand_left(;id::Int64, s_dist::Distribution, v_dist::Distribution) = get_rand_Tint_vehicle(id=id, s_dist=s_dist, v_dist=v_dist, lane_dist = [1,2])
rand_right(;id::Int64, s_dist::Distribution, v_dist::Distribution) = get_rand_Tint_vehicle(id=id, s_dist=s_dist, v_dist=v_dist, lane_dist = [3,4])

## Pedestrian Crosswalk scenario helpers
ez_ped_vehicle(;id::Int64, s::Float64, v::Float64) = BlinkerVehicle(roadway = ped_roadway,
                                              lane=1, s=s, v = v,
                                              id = id, goals = ped_goals[1],
                                              blinker = ped_should_blink[1])
ez_pedestrian(;id::Int64, s::Float64, v::Float64) = NoisyPedestrian(roadway = ped_roadway, lane = 2, s=s, v=v, id=id)

get_ped_vehicle(;id::Int64, s::Float64, v::Float64) = (rng::AbstractRNG = Random.GLOBAL_RNG) -> ez_ped_vehicle(id=id, s=s, v=v)
get_rand_ped_vehicle(;id::Int64, s_dist, v_dist) = (rng::AbstractRNG = Random.GLOBAL_RNG) -> ez_ped_vehicle(id=id, s=rand(rng, s_dist), v=rand(rng, v_dist))
get_pedestrian(;id::Int64, s::Float64, v::Float64) = (rng::AbstractRNG = Random.GLOBAL_RNG) -> ez_pedestrian(id=id, s=s, v=v)
get_rand_pedestrian(;id::Int64, s_dist, v_dist) = (rng::AbstractRNG = Random.GLOBAL_RNG) -> ez_pedestrian(id=id, s=rand(rng, s_dist), v=rand(rng, v_dist))

## Create gif from rollout

function scenes_to_gif(scenes, roadway, filename; others = [], others_fn = [], fps = 10, camera = nothing, canvas = (1200, 800), repeat_last_frame = 1, verbose = false )
    push!(scenes, fill(scenes[end], repeat_last_frame)...)
    files = []
    for i=1:length(scenes)
        verbose && println("frame $i")
        frame = AutomotiveVisualization.render([roadway, others..., [f(scenes[i]) for f in others_fn]..., scenes[i]], canvas_width=canvas[1], canvas_height=canvas[2], camera = camera)
        file = string("/tmp/test_", lpad(i,4, "0"), ".png")
        push!(files, file)
        write(file, frame)
    end
    delay = 100 / fps
    run(`convert -delay $delay /tmp/test_"*".png  $filename`)
    run(`rm $files`)
end

# Create a random IntelligentDriverModel
function random_IDM(rng::AbstractRNG)
    headway_t = rand(rng, Uniform(0.5, 1.5)) # desired time headway [s]
    v_des = rand(rng, Uniform(10., 15.))
    s_min = rand(rng, Uniform(2.0, 5.0)) # minimum acceptable gap [m]
    d_cmf = rand(rng, Uniform(1.5, 2.5)) # comfortable deceleration [m/s²] (positive)
    d_max = rand(rng, Uniform(5.0, 9.0)) # maximum deceleration [m/s²] (positive)
    IntelligentDriverModel(T = headway_t, v_des = v_des, s_min = s_min, d_cmf = d_cmf, d_max = d_max)
end

## action conversion function

function continuous_policy(t, rng = Random.GLOBAL_RNG)
    function get_action(s)
        noise = Noise((rand(rng, t[:noise_s]), 0.), rand(rng, t[:noise_v]))
        Disturbance[BlinkerVehicleControl(0., rand(rng, t[:da]), rand(rng, t[:toggle_goal]), rand(rng, t[:toggle_blinker]), noise)]
    end
    FunctionPolicy(get_action)
end

# Converts a Multi-variate timeseries to the actions of BlinkerVehicle
function create_actions_1BV(d)
    actions = Vector{Disturbance}[]
    for i=1:length(first(d)[2])
        noise = Noise((d[:noise_s][i], 0.), d[:noise_v][i])
        a = BlinkerVehicleControl(0., d[:da][i], d[:toggle_goal][i], d[:toggle_blinker][i], noise)
        push!(actions, [a])
    end
    actions
end

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
struct SemiconductorPoint{T}
    Data::SemiconductorData
    p::T
end

SemiconductorPoint1D(Data::SemiconductorData,p::Real=0.0) = SemiconductorPoint(Data,p)

function SemiconductorPoint1D(Data::SemiconductorData,p::Real,ψ,φfₑ=0.0*Eunit,φfₚ=0.0*Eunit)
    D = SemiconductorData(Data,ψ,φfₑ,φfₚ)
    SemiconductorPoint1D(D,p)
end

function SemiconductorPoint1D(S1D::SemiconductorPoint,p)
    SemiconductorPoint1D(S1D.Data,p)
end

Eg(S::SemiconductorPoint) = Eg(S.Data)
Ec(S::SemiconductorPoint) = Ec(S.Data)
Ev(S::SemiconductorPoint) = Ev(S.Data)
Ef(S::SemiconductorPoint) = Ef(S.Data)
Ei(S::SemiconductorPoint) = Ei(S.Data)
#= n(S::SemiconductorPoint) = n(S.Data)
p(S::SemiconductorPoint) = p(S.Data) =#
n(S::SemiconductorPoint,ψ=0.0*Eunit,φfₑ=0.0*Eunit) = n(S.Data,ψ,φfₑ)
p(S::SemiconductorPoint,ψ=0.0*Eunit,φfₚ=0.0*Eunit) = p(S.Data,ψ,φfₚ)
Na(S::SemiconductorPoint) = Na(S.Data)
Nd(S::SemiconductorPoint) = Nd(S.Data)
N(S::SemiconductorPoint) = N(S.Data)
N(S::SemiconductorPoint,ψ) = N(S.Data,ψ)
P(S::SemiconductorPoint) = P(S.Data)
P(S::SemiconductorPoint,ψ) = P(S.Data,ψ)
Pt(S::SemiconductorPoint) = Pt(S.Data)
Nt(S::SemiconductorPoint) = Nt(S.Data)
ni(S::SemiconductorPoint) = ni(S.Data)
φf(S::SemiconductorPoint) = φf(S.Data)
r(S::SemiconductorPoint) = S.p
ψ(S::SemiconductorPoint) = S.Data.ψ
φfₑ(S::SemiconductorPoint) = S.Data.φfₑ
φfₚ(S::SemiconductorPoint) = S.Data.φfₚ
ϵᵣ(S::SemiconductorPoint) = ϵᵣ(S.Data)
ξ(S::SemiconductorPoint) = ξ(S.Data,ψ(S))
ξ(S::SemiconductorPoint,ψ) = ξ(S.Data,ψ)
ξ(S::SemiconductorPoint,ψ,ϕ) = ξ(S.Data,ψ,ϕ)
ρ(S::SemiconductorPoint) = ρ(S.Data)
ρ(S::SemiconductorPoint,ψ) = ρ(S.Data,ψ)
Δχ(S::SemiconductorPoint,δχ) = Δχ(S.Data,δχ)
χ(S::SemiconductorPoint) = χ(S.Data)
Temperature(S::SemiconductorPoint) = Temperature(S.Data)
Nv(S::SemiconductorPoint) = Nv(S.Data)
Nc(S::SemiconductorPoint) = Nc(S.Data)
DistFun(S::SemiconductorPoint,E) = DistFun(S.Data,E)
C(S::SemiconductorPoint) = C(S.Data)

function applyPsi(S::SemiconductorPoint,ψ)
    D = S.Data
    _p = S.p
    _φfₚ = φfₚ(S)
    _φfₑ = φfₑ(S)
    return SemiconductorPoint1D(D,_p,ψ,_φfₑ,_φfₚ)
end

function applyEfe(S::SemiconductorPoint,φfₑ)
    D = S.Data
    _p = S.p
    psi = ψ(S)
    _φfₚ = φfₚ(S)
    return SemiconductorPoint1D(D,_p,psi,φfₑ,_φfₚ)
end

function applyEfp(S::SemiconductorPoint,φfₚ)
    D = S.Data
    _p = S.p
    psi = ψ(S)
    _φfₑ = φfₑ(S)
    return SemiconductorPoint1D(D,_p,psi,_φfₑ,φfₚ)
end


text:
ENV["COMPILE_STATIC"] = true
include(joinpath(Pkg.dir("PackageCompiler"), "juliac.jl"))

args = ["-vcasj", "--cc-flags=-Wl,-rpath,\$ORIGIN", "--compile=all", "foo_AOT.jl"]

print_with_color(:cyan, "juliac.jl $(join(args, " "))\n")
julia_main(args)

includePath = joinpath(JULIA_HOME, Base.INCLUDEDIR, "julia")
command = `g++ -shared -o foo_driver.so ../driver/foo_driver.cpp -std=c++11 -I$includePath -DJULIA_ENABLE_THREADING=1 -fPIC -Wl,--export-dynamic -m64 foo_AOT.so '-Wl,-rpath,$ORIGIN'`
print_with_color(:cyan, command, "\n")
run(command)

command = `g++ -o app ../app/app.cpp -std=c++11 -fPIC -m64 foo_driver.so '-Wl,-rpath,$ORIGIN' -fuse-ld=gold`
print_with_color(:cyan, command, "\n")
run(command)

command = `./app`
print_with_color(:cyan, command, "\n")
run(command)


text:
function update_v(n_obs::Int, Z::Float64)
    return rand(Gamma(n_obs, 1 / Z))
end

function update_M!(M::Array, γ::Array, K::Int, N::Int)
    # Update the mass parameter
    prior = [2.0, 0.25]
    for k = 1:K
        @inbounds current_γ = γ[:, k]
        current_M = M[k]
        log_likelihood = - sum(logpdf.(Gamma(current_M / N, 1.0), current_γ))
        log_likelihood_0 = - sum(logpdf.(Gamma(prior[1], prior[2]), current_M))
        proposed_mass = (current_M + rand(Normal()) / 10)
        if proposed_mass <= 0.0
            alpha = 0.0
        else
            new_log_likelihood = - sum(logpdf.(Gamma(proposed_mass / N, 1), current_γ))
            new_log_likelihood_0 = - sum(logpdf(Gamma.(prior[1], prior[2]), proposed_mass))
            alpha = exp(-new_log_likelihood - new_log_likelihood_0 + log_likelihood + log_likelihood_0)
        end
        if rand() < alpha
            M[k] = proposed_mass
        end
    end
    return
end

@inline function update_Z(Φ::Array, Φ_index::Array, Γ::Array)
    # Update the normalising constant
    # Z = sum(exp.((Φ_index * (log.(Φ .+ 1)) + sum(Γ, dims = 2))))
    Φ_log = log.(Φ .+ 1)
    norm_temp = Φ_index * Φ_log + sum(Γ, dims = 2)
    Z = 0.0
    for i in eachindex(norm_temp)
        Z += exp(norm_temp[i])
    end
    return Z
end

function calculate_likelihood(s::Array, Φ::Array, γ::Array, Z::Float64)
    likelihood = zeros(Float64, (size(s, 1)))
    Φ_log = similar(Φ)
    for i in eachindex(Φ)
        Φ_log[i] = log(1 + Φ[i])
    end
    for i in 1:size(s, 1)
        for k in 1:size(s, 2)
            likelihood[i] += log(γ[s[i, k], k])
        end
        if size(s, 2) > 1
            ϕ = 1
            for k1 in 1:(size(s, 2) - 1)
                for k2 in (k1 + 1):(size(s, 2))
                    likelihood .+= Φ_log[ϕ] * s[i, k1] == s[i, k2]
                    ϕ += 1
                end
            end
        end
    end
    return sum(exp.(likelihood) ./ Z)
end

function update_γ!(γ::Array, Φ::Array, v::Float64, M, s::Array, Φ_index::Array, γ_combn::Array, Γ::Array, N::Int, K::Int)
    β_0 = 1.0
    Φ_log = log.(Φ .+ 1)
    α_star = Matrix{Float64}(undef, N, K)
    for k = 1:K
        for n = 1:N
            # @inbounds α_star[n, k] = α_0 + sum(s_k .== n)
            # @inbounds α_star[n, k] = M[k] / N + sum(s[:, k] .== n)
            @inbounds α_star[n, k] = M[k] / N + countn(s[:, k], n)
        end
    end
    norm_temp = Φ_index * Φ_log + sum(Γ, dims = 2)
    for i in eachindex(norm_temp)
        norm_temp[i] = exp(norm_temp[i])
    end
    @inbounds for k = 1:K
         pertinent_rows = findZindices(k, K, 1, N)
         for n = 1:N
            old_γ = γ[n, k] + 0.0
            β_star = β_0 + v * sum(view(norm_temp, pertinent_rows, :)) / γ[n, k]
            γ[n, k] = rand(Gamma(α_star[n, k], 1 / β_star)) + eps(Float64)
            for i in pertinent_rows
                norm_temp[i] *= γ[n, k] / old_γ
            end
            pertinent_rows .+= N ^ (k - 1)
        end
    end
    return
end

function update_Φ!(Φ, v::Float64, s, Φ_index, γ, K::Int, Γ)
    # Prior parameters
    α_0 = 1.0
    β_0 = 5.0
    Φ_lab = calculate_Φ_lab(K)
    Φ_log = log.(Φ .+ 1)
    norm_temp = Φ_index * Φ_log + sum(Γ, dims = 2)
    for i in eachindex(norm_temp)
        @inbounds norm_temp[i] = exp(norm_temp[i])
    end
    @inbounds for i in 1:length(Φ)
        # Get relevant allocations
        Φ_current = Φ[i] + 0.0
        n_agree = 0
        for j in 1:size(s, 1)
            # For some reason checking if the difference is zero
            # is quicker/fewer allocations than checking if the
            # two are equal
            # Perhaps it removes checks for types?
            n_agree += (s[j, Φ_lab[i, 1]] - s[j, Φ_lab[i, 2]]) == 0
        end
        # Get relevant terms in the normalisation constant
        pertinent_rows = findall(Φ_index[:, i])
        β_star = β_0 + (v * sum(norm_temp[pertinent_rows, :]) / (1 + Φ_current))
        weights = lgamma.((0:n_agree) .+ α_0)
        weights += logpdf.(Binomial(n_agree, 0.5), 0:n_agree)
        weights -= (0:(n_agree)) .* log(1 / β_star)
        α_star = α_0 + sample(0:n_agree, Weights(exp.(weights .- maximum(weights))))
        Φ[i] = rand(Gamma(α_star, 1 / β_star)) # + eps(Float64)
        # Update the normalising constant values to account for this update
        norm_temp[pertinent_rows, :] .*= (1 + Φ[i]) / (1 + Φ_current)
    end
    return
end


text:
#Contains report global properties
type Report <: Display
  cwd::AbstractString
  basename::AbstractString
  formatdict::Dict{Symbol,Any}
  pending_code::AbstractString
  cur_result::AbstractString
  rich_output::AbstractString
  fignum::Int
  figures::Array{AbstractString}
  term_state::Symbol
  cur_chunk
  mimetypes::Array{AbstractString}
end

function Report(cwd, basename, formatdict, mimetypes)
    Report(cwd, basename, formatdict, "", "", "", 1, AbstractString[], :text, nothing, mimetypes)
end

#Default mimetypes in order, can be overridden for some inside `run method` formats
const default_mime_types = ["image/svg+xml", "image/png", "text/html", "text/plain"]
#From IJulia as a reminder
#const supported_mime_types = [ "text/html", "text/latex", "image/svg+xml", "image/png", "image/jpeg", "text/plain", "text/markdown" ]

function Base.display(report::Report, data)
    #Set preferred mimetypes for report based on format
    for m in report.mimetypes
        if mimewritable(m, data)
            display(report, m, data)
            #Always show plain text as well for term mode
            if m ≠ "text/plain" && report.cur_chunk.options[:term]
              display(report, "text/plain", data)
            end
            break
        end
    end
end

function Base.display(report::Report, m::MIME"image/png", data)
    figname = add_figure(report, ".png")
    open(figname, "w") do io
      writemime(io, m, data)
    end
end

function Base.display(report::Report, m::MIME"image/svg+xml", data)
    figname = add_figure(report, ".svg")
    open(figname, "w") do io
      writemime(io, m, data)
    end
end

function Base.display(report::Report, m::MIME"application/pdf", data)
    figname = add_figure(report, ".pdf")
    open(figname, "w") do io
      writemime(io, m, data)
    end
end

#Text is written to stdout, called from "term" mode chunks
function Base.display(report::Report, m::MIME"text/plain", data)
    s = reprmime(m, data)
    print("\n" * s)
end

#Catch "rich_output"
function Base.display(report::Report, m::MIME"text/html", data)
    s = reprmime(m, data)
    report.rich_output *= "\n" * s
end

#Catch "rich_output"
function Base.display(report::Report, m::MIME"text/markdown", data)
    s = reprmime(m, data)
    report.rich_output *= "\n" * s
end

function Base.display(report::Report, m::MIME"text/latex", data)
    s = reprmime(m, data)
    report.rich_output *= "\n" * s
end

"""Add saved figure name to results and return the name"""
function add_figure(report::Report, ext)
  chunk = report.cur_chunk
  full_name, rel_name = get_figname(report, chunk, ext = ext)
  push!(report.figures, rel_name)
  report.fignum += 1
  return full_name
end


text:
using Documenter
using ClimateSatellite

makedocs(
    modules  = [ClimateSatellite],
    doctest  = false,
    format   = Documenter.HTML(
        collapselevel = 1,
        prettyurls    = false
    ),
    authors  = "Nathanael Wong",
    sitename = "ClimateSatellite.jl",
    pages    = [
        "Home"      => "index.md",
        "Tutorials" => [
            "Data Downloads"       => "tutorials/downloads.md",
            "Using GeoRegions"     => "tutorials/georegions.md",
            # "Region Extraction"    => "tutorials/extract.md",
            # "Preliminary Analysis" => "tutorials/analysis.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/JuliaClimate/ClimateSatellite.jl.git",
)


text:
# libpurex:
# oracle computations

# a Pure Exploration problem (pep) is parameterised by
# - a domain, \mathcal M, representing the prior knowledge about the structure
# - a query, as embodied by a correct-answer function istar
#
# To specify a PE problem here, we need to compute the following things:
# - nanswers: number of possible answers
# - istar: correct answer for feasible μ
# - glrt: value and best response (λ and ξ) to (N, ̂μ) or (w, μ)
# - oracle: characteristic time and oracle weights at μ

struct BestArm
    expfam; # common exponential family
end

nanswers(pep::BestArm, μ) = length(μ);
istar(pep::BestArm, μ) = argmax(μ);
getexpfam(pep::BestArm, k) = pep.expfam;

function glrt(pep::BestArm, w, μ)
    @assert length(size(μ)) == 1

    ⋆ = argmax(μ); # index of best arm among μ

    val, k, θ = minimum(
        begin
        # transport ⋆ and k to weighted midpoint
        θ = (w[⋆]*μ[⋆]+w[k]*μ[k])/(w[⋆]+w[k]);
        w[⋆]*d(pep.expfam, μ[⋆], θ) + w[k]*d(pep.expfam, μ[k], θ), k, θ
        end
        for k in eachindex(μ)
        if k != ⋆
    );

    λ = copy(μ);
    λ[⋆] = θ;
    λ[k] = θ;

    # note: ξs = μ here

    val, (k, λ), (⋆, μ);
end

# solve for x such that d(μ1, μx) + x*d(μi, μx) == v
# where μx = (μ1+x*μi)/(1+x)
function X(expfam, μ1, μi, v)
    kl1i = d(expfam, μ1, μi); # range of V(x) is [0, kl1i]
    @assert 0 ≤ v ≤ kl1i "0 ≤ $v ≤ $kl1i";
    α = binary_search(
        z -> let μz = (1-z)*μ1+z*μi
        (1-z)*d(expfam, μ1, μz) + z*d(expfam, μi, μz) - (1-z)*v
        end,
        0, 1);
    α/(1-α);
end

# oracle problem
function oracle(pep::BestArm, μs)
    μstar = maximum(μs);

    if all(μs .== μstar) # yes, this happens
        return Inf, ones(length(μs))/length(μs);
    end

    # determine upper range for subsequent binary search
    hi = minimum(
        d(pep.expfam, μstar, μ)
        for μ in μs
        if μ != μstar
    );

    val = binary_search(
        z -> sum(
            let x = X(pep.expfam, μstar, μ, z),
               μx = (μstar + x*μ)/(1+x);
            d(pep.expfam, μstar, μx) / d(pep.expfam, μ, μx)
            end
            for μ in μs
            if μ != μstar
            ) - 1,
        0, hi);

    ws = [(μ == μstar) ? 1. : X(pep.expfam, μstar, μ, val) for μ in μs];
    Σ = sum(ws);
    Σ/val, ws ./ Σ;
end

# oracle problem for best μ in confidence interval
function optimistic_oracle(pep::BestArm, hμ, N)

    t = sum(N);

    μdn = [ddn(pep.expfam, hμ[k], log(t)/N[k]) for k in eachindex(hμ)];
    μup = [dup(pep.expfam, hμ[k], log(t)/N[k]) for k in eachindex(hμ)];

    # try to make each arm the best-looking arm so far,
    # then move everybody else down as far as possible
    minimum(oracle(pep,
                   ((k == j) ? μup[k] : μdn[k]
                    for k in eachindex(hμ)))
            for j in eachindex(hμ)
            if μup[j] > maximum(μdn)
            );
end

struct LargestProfit
    expfam;
end

nanswers(pep::LargestProfit, μ) = length(μ)/2;

function istar(pep::LargestProfit, _μs)
    μs = reshape(_μs, 2, :);
    argmax(μs[1,:]-μs[2,:]); # current best
end
getexpfam(pep::LargestProfit, k) = pep.expfam;

# NOTE this does NOT work for divergences that are non-convex in the
# second argument. Moreover, the problem here seems fundamental.
function glrt(pep::LargestProfit, _ws, _μs)

    # convert to 2×K indexing
    ws = reshape(_ws, 2, :);
    μs = reshape(_μs, 2, :);

    K = size(ws, 2);

    ⋆ = argmax(μs[1,:]-μs[2,:]); # current best

    val, k, λ1k, λ2k, λ1s, λ2s = minimum(
        begin
        # arm k needs to be the best

        xmax = 1000; # voodoo
        x = binary_search(x ->
                          + invh(pep.expfam, μs[1,k], +x/ws[1,k])
                          - invh(pep.expfam, μs[2,k], -x/ws[2,k])
                          - invh(pep.expfam, μs[1,⋆], -x/ws[1,⋆])
                          + invh(pep.expfam, μs[2,⋆], +x/ws[2,⋆]), 0, xmax);

        λ1k = invh(pep.expfam, μs[1,k], +x/ws[1,k]);
        λ2k = invh(pep.expfam, μs[2,k], -x/ws[2,k]);
        λ1s = invh(pep.expfam, μs[1,⋆], -x/ws[1,⋆]);
        λ2s = invh(pep.expfam, μs[2,⋆], +x/ws[2,⋆]);

        ws[1,⋆]*d(pep.expfam, μs[1,⋆], λ1s) +
        ws[2,⋆]*d(pep.expfam, μs[2,⋆], λ2s) +
        ws[1,k]*d(pep.expfam, μs[1,k], λ1k) +
        ws[2,k]*d(pep.expfam, μs[2,k], λ2k), k, λ1k, λ2k, λ1s, λ2s
        end
        for k in 1:K
        if k != ⋆
    );
    _λs = copy(_μs);
    λs = reshape(_λs, 2, :);
    λs[1,⋆] = λ1s;
    λs[2,⋆] = λ2s;
    λs[1,k] = λ1k;
    λs[2,k] = λ2k;
    # note: ξs = μs here
    val, (k, _λs), (⋆, _μs);
end

function oracle(pep::LargestProfit, μs)
    @warn "oracle for LargestProfit not yet implemented";
    nothing, nothing;
end

struct MinimumThreshold
    expfam;
    γ;
end

nanswers(pep::MinimumThreshold, μ) = 2;
istar(pep::MinimumThreshold, μs) = 1+(minimum(μs) > pep.γ);
getexpfam(pep::MinimumThreshold, k) = pep.expfam;

# answers are encoded as 1 = "<", 2 = ">"
function glrt(pep::MinimumThreshold, ws, μs)
    if minimum(μs) < pep.γ
        #everybody moves up
        λs = max.(μs, pep.γ);
        val = sum(k -> ws[k]*d(pep.expfam, μs[k], λs[k]), eachindex(ws));
        val, (2, λs), (1, μs);
    else
        #someone moves down
        val, k = minimum(
            (ws[k]*d(pep.expfam, μs[k], pep.γ), k)
            for k in eachindex(ws));
        λs = copy(μs);
        λs[k] = pep.γ;
        val, (1, λs), (2, μs);
    end
end

function oracle(pep::MinimumThreshold, μs)
    μstar = minimum(μs);
    if μstar < pep.γ
        # everything on lowest arm
        ws = Float64.(μs .== μstar);
        Tstar = 1/d(pep.expfam, μstar, pep.γ);
    else
        ws = [1/d(pep.expfam, μ, pep.γ) for μ in μs];
        Tstar = sum(ws);
    end
    Tstar, ws ./ sum(ws);
end

# oracle problem for best μ in confidence interval
function optimistic_oracle(pep::MinimumThreshold, hμ, N)

    t = sum(N);

    μdn = (ddn(pep.expfam, hμ[k], log(t)/N[k]) for k in eachindex(hμ));
    μup = (dup(pep.expfam, hμ[k], log(t)/N[k]) for k in eachindex(hμ));

    if minimum(μup) < pep.γ
        # entire conf. interval below
        oracle(pep, μdn);
    elseif minimum(μdn) > pep.γ
        # entire conf. interval above
        oracle(pep, μup);
    else
        # both options in conf. interval. Take best
        min(oracle(pep, μup),
            oracle(pep, μdn));
    end
end

# Best arm with a simplex constraint on μ
struct BestArmSimplex
end

# all arms on same side of threshold
# interesting because ̂μ can be outside
struct SameSign
    expfam;
    γ; # threshold
end

nanswers(pep::SameSign, μ) = 2;
# answer 1 is everyone below,
# answer 2 is everyone above

function istar(pep::SameSign, μ)
    if all(μ .≤ pep.γ)
        1
    elseif all(μ .≥ pep.γ)
        2
    else
        @assert false "infeasible bandit model $μ";
    end
end
getexpfam(pep::SameSign, k) = pep.expfam;

function glrt(pep::SameSign, w, μ)
    μlo = min.(pep.γ, μ); # lower everyone
    μhi = max.(pep.γ, μ); # raise everyone
    dlo = sum(k -> w[k]*d(pep.expfam, μ[k], μlo[k]), eachindex(μ));
    dhi = sum(k -> w[k]*d(pep.expfam, μ[k], μhi[k]), eachindex(μ));

    if dlo < dhi
        dhi-dlo, (2, μhi), (1, μlo);
    else
        dlo-dhi, (1, μlo), (2, μhi);
    end
end

function oracle(pep::SameSign, μ)
    @assert all(μ .≤ pep.γ) || all(μ .≥ pep.γ) "infeasible bandit model $μ";

    v, i = findmax([d(pep.expfam, μ[k], pep.γ) for k in eachindex(μ)]);
    1/v, Float64.(eachindex(μ).==i)
end

# two arms with equal means but different exponential family (e.g. variance).
# question if they are above/below the threshold γ
# do algorithms pull the low variance one too often?
struct HeteroSkedastic
    expfams; # two(!) exponential families
    γ;
end

nanswers(pep::HeteroSkedastic, μ) = 2;
function istar(pep::HeteroSkedastic, μ)
    @assert length(μ) == 2;
    @assert abs(μ[1] - μ[2]) < 1e-4 "unequal $μ";
    1 + (μ[1] > pep.γ);
end

getexpfam(pep::HeteroSkedastic, k) = pep.expfams[k];

function glrt(pep::HeteroSkedastic, w, μ)

    # compute ξ
    xmax = 400; # voodoo
    if minimum(w) > 1e-5
        x = binary_search(x ->
                          + invh(pep.expfams[1], μ[1], +x/w[1])
                          - invh(pep.expfams[2], μ[2], -x/w[2]), -xmax, xmax);

        ξ = [invh(pep.expfams[1], μ[1], +x/w[1]),
             invh(pep.expfams[2], μ[2], -x/w[2])];
    else
        ξ = [w'μ, w'μ]; # two copies of the average
    end

    ⋆ = istar(pep, ξ);
    λ = [pep.γ, pep.γ];
    val = sum(k -> w[k]*(d(pep.expfams[k], μ[k], λ[k]) -
                         d(pep.expfams[k], μ[k], ξ[k])), 1:2);
    val, (3-⋆, λ), (⋆, ξ);
end

function oracle(pep::HeteroSkedastic, μ)
    v, i = findmax([d(pep.expfams[k], μ[k], pep.γ) for k in eachindex(μ)]);
    1/v, Float64.(eachindex(μ).==i)
end


text:
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot
using Printf
include("utilities.jl")

series_1 = Polylogarithms.polylog_series_1
series_2 = Polylogarithms.polylog_series_2
series_3 = Polylogarithms.polylog_series_3
polylog_reciprocal = Polylogarithms.polylog_reciprocal

# test behaviour along the branch
s = -1.5
z = 3.0 .+ collect( -0.01: 0.001: 0.01 )*im

for i=1:length(z)
    println(" ", imag(z[i]), ", ", imag( series_2(s, z[i])[1]), ", ", imag( polylog_reciprocal(s, z[i])[1]) )
end


text:
@testset "Preprocessing" begin

    sample_text1 = "This is 1 MESSED υπ string!"
    sample_text1_wo_punctuation = "This is 1 MESSED υπ string"
    sample_text1_wo_punctuation_numbers = "This is  MESSED υπ string"
    sample_text1_wo_punctuation_numbers_case = "this is  messed υπ string"
    sample_text1_wo_punctuation_numbers_case_az = "this is  messed  string"

    sample_texts = [
        sample_text1,
        sample_text1_wo_punctuation,
        sample_text1_wo_punctuation_numbers,
        sample_text1_wo_punctuation_numbers_case,
        sample_text1_wo_punctuation_numbers_case_az
    ]

    # This idiom is _really_ ugly since "OR" means "AND" here.
    for str in sample_texts
        sd = StringDocument(str)
        prepare!(
            sd,
            strip_punctuation | strip_numbers | strip_case | strip_whitespace | strip_non_letters
        )
        @test isequal(strip(sd.text), "this is messed string")
    end

    # Need to only remove words at word boundaries
    doc = Document("this is sample text")
    remove_words!(doc, ["sample"])
    @test isequal(doc.text, "this is  text")

    doc = Document("this is sample text")
    prepare!(doc, strip_articles)
    @test isequal(doc.text, "this is sample text")

    doc = Document("this is sample text")
    prepare!(doc, strip_definite_articles)
    @test isequal(doc.text, "this is sample text")

    doc = Document("this is sample text")
    prepare!(doc, strip_indefinite_articles)
    @test isequal(doc.text, "this is sample text")

    doc = Document("this is sample text")
    prepare!(doc, strip_prepositions)
    @test isequal(doc.text, "this is sample text")

    doc = Document("this is sample text")
    prepare!(doc, strip_pronouns)
    @test isequal(doc.text, "this is sample text")

    doc = Document("this is sample text")
    prepare!(doc, strip_stopwords)
    @test isequal(strip(doc.text), "sample text")

    doc = Document("this is sample text")
    prepare!(doc, strip_whitespace)
    @test isequal(doc.text, "this is sample text")

    # stem!(sd)

    # Do preprocessing on TokenDocument, NGramDocument, Corpus
    d = NGramDocument("this is sample text")
    @test haskey(d.ngrams, "sample")
    remove_words!(d, ["sample"])
    @test !haskey(d.ngrams, "sample")

    d = StringDocument(
        """

Hello
 world

        """
    )
    remove_html_tags!(d)
    @test "Hello world" == strip(d.text)

    style_html_doc = StringDocument(
      """

Hello
 world

      """
     )
    remove_html_tags!(style_html_doc)
    @test "Hello world" == strip(style_html_doc.text)

    #Test #62
    remove_corrupt_utf8("abc") == "abc"
    remove_corrupt_utf8(String([0x43, 0xf0])) == "C "

    #Test frequent_terms sparse_terms
    crps = Corpus(StringDocument.(sample_texts))
    @test isempty(setdiff(frequent_terms(crps),["string","is"]))
    @test isempty(setdiff(sparse_terms(crps,0.3),["!"]))

    #Tests strip_punctuation regex conditions
    str = Document("These punctuations should be removed [-.,:;,!?'\"[](){}|\`#\$%@^&*_+<>")
    answer = Document("These punctuations should be removed ")
    prepare!(str, strip_punctuation)
    @test isequal(str.text, answer.text)

    str = Document("Intel(tm) Core i5-3300k, is a geat CPU! ")
    answer = Document("Inteltm Core i53300k is a geat CPU ")   #tests old implementation
    prepare!(str, strip_punctuation)
    @test isequal(str.text, answer.text)

    #Tests no whitespace at end or begining
    doc = Document("   this is sample text   ")
    prepare!(doc, strip_whitespace)
    @test isequal(doc.text, "this is sample text")

    doc = Document("   this is sample text   ")
    prepare!(doc, strip_whitespace)
    @test isequal(doc.text, "this is sample text")

    crps = Corpus(
            [StringDocument("         Document      1"),
            StringDocument("       Document      2              ")]
           )
    prepare!(crps, strip_whitespace)
    @test isequal(crps[1].text, "Document 1")
    @test isequal(crps[2].text, "Document 2")

    crps = Corpus([StringDocument("     Hi     there    !     ")])
    prepare!(crps, strip_html_tags | strip_whitespace | strip_non_letters)
    @test isequal(crps[1].text, "Hi there")
end


text:
isapproxzero(x::T; kws...) where {T<:Real} = x == zero(T)
isapproxzero(x::T; ztol=Base.rtoldefault(T)) where {T<:AbstractFloat} = abs(x) < ztol
isapproxzero(x::AbstractVector{T}; kws...) where {T<:Real} = isapproxzero(maximum(abs.(x)); kws...)
#isapproxzero(x::Union{SymPoint, Ray, Line}; kws...) = isapproxzero(coord(x); kws...)
isapproxzero(x::Union{Ray, Line}; kws...) = isapproxzero(coord(x); kws...)
isapproxzero(h::HRepElement; kws...) = isapproxzero(h.a; kws...) && isapproxzero(h.β; kws...)

_isapprox(x::Union{T, AbstractArray{T}}, y::Union{T, AbstractArray{T}}) where {T<:Union{Integer, Rational}} = x == y
# I check with zero because isapprox(0, 1e-100) is false...
# but isapprox(1e-100, 2e-100) should be false
_isapprox(x, y) = (isapproxzero(x) ? isapproxzero(y) : (isapproxzero(y) ? isapproxzero(x) : isapprox(x, y)))

#Base.isapprox(r::SymPoint, s::SymPoint) = _isapprox(coord(r), coord(s)) || _isapprox(coord(r), -coord(s))
function _scaleray(r::Union{Line, Ray}, s::Union{Line, Ray})
    cr = coord(r)
    cs = coord(s)
    cr * sum(abs.(cs)), cs * sum(abs.(cr))
end
function Base.isapprox(r::Line, s::Line)
    rs, ss = _scaleray(r, s)
    _isapprox(rs, ss) || _isapprox(rs, -ss)
end
function Base.isapprox(r::Ray, s::Ray)
    _isapprox(_scaleray(r, s)...)
end
# TODO check that Vec in GeometryTypes also does that
function _scalehp(h1, h2)
    s1 = sum(abs.(h1.a)) + abs(h1.β)
    s2 = sum(abs.(h2.a)) + abs(h2.β)
    (h1.a*s2, h1.β*s2), (h2.a*s1, h2.β*s1)
end
function Base.:(==)(h1::HyperPlane, h2::HyperPlane)
    (a1, β1), (a2, β2) = _scalehp(h1, h2)
    (a1 == a2 && β1 == β2) || (a1 == -a2 && β1 == -β2)
end
function Base.isapprox(h1::HyperPlane, h2::HyperPlane)
    (a1, β1), (a2, β2) = _scalehp(h1, h2)
    (_isapprox(a1, a2) && _isapprox(β1, β2)) || (_isapprox(a1, -a2) && _isapprox(β1, -β2))
end
function Base.:(==)(h1::HalfSpace, h2::HalfSpace)
    (a1, β1), (a2, β2) = _scalehp(h1, h2)
    a1 == a2 && β1 == β2
end
function Base.isapprox(h1::HalfSpace, h2::HalfSpace)
    (a1, β1), (a2, β2) = _scalehp(h1, h2)
    _isapprox(a1, a2) && _isapprox(β1, β2)
end

_lt(x::T, y::T) where {T<:Real} = x < y
_lt(x::T, y::T) where {T<:AbstractFloat} = x < y && !_isapprox(x, y)
_lt(x::S, y::T) where {S<:Real,T<:Real} = _lt(promote(x, y)...)
_gt(x::S, y::T) where {S<:Real, T<:Real} = _lt(y, x)
_leq(x::T, y::T) where {T<:Real} = x <= y
_leq(x::T, y::T) where {T<:AbstractFloat} = x <= y || _isapprox(x, y)
_leq(x::S, y::T) where {S<:Real,T<:Real} = _leq(promote(x, y)...)
_geq(x::T, y::T) where {T<:Real} = _leq(y, x)
_pos(x::T) where {T<:Real} = _gt(x, zero(T))
_neg(x::T) where {T<:Real} = _lt(x, zero(T))
_nonneg(x::T) where {T<:Real} = _geq(x, zero(T))
_nonpos(x::T) where {T<:Real} = _leq(x, zero(T))

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
using Documenter, LibPQ2, Memento

setlevel!(getlogger(LibPQ2), "critical")

DocMeta.setdocmeta!(LibPQ2, :DocTestSetup, quote
    using LibPQ2
    DATABASE_USER = get(ENV, "LIBPQJL_DATABASE_USER", "postgres")
    conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER")
end; recursive=true)

makedocs(;
    modules=[LibPQ2],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Type Conversions" => "pages/type-conversions.md",
        "API" => "pages/api.md",
        "FAQ" => "pages/faq.md",
    ],
    repo="https://github.com/arkoniak/LibPQ2.jl/blob/{commit}{path}#L{line}",
    sitename="LibPQ2.jl",
    checkdocs=:exports,
    linkcheck=true,
    linkcheck_timeout=60,
    strict=true,
    authors="Eric Davies",
)

deploydocs(;
    repo="github.com/arkoniak/LibPQ2.jl",
)


text:
include("BUDE.jl")
using CUDA, StaticArrays

function devices()::Vector{DeviceWithRepr}
  return !CUDA.functional(false) ? [] :
         map(d -> (d, "$(CUDA.name(d)) ($(repr(d)))", "CUDA.jl"), CUDA.devices())
end

function run(params::Params, deck::Deck, device::DeviceWithRepr)
  # show_reason is set to true here so it dumps CUDA info
  # for us regardless of whether it's functional
  if !CUDA.functional(true)
    error("Non-functional CUDA configuration")
  end

  # so CUDA's device is 0 indexed, so -1 from Julia
  CUDA.device!(device[1])
  CUDA.math_mode!(CUDA.FAST_MATH)

  println("Using CUDA device: $(CUDA.name(device[1])) ($(repr(device[1])))")

  protein = CuArray{Atom}(deck.protein)
  ligand = CuArray{Atom}(deck.ligand)
  forcefield = CuArray{FFParams}(deck.forcefield)
  poses = CuArray{Float32,2}(deck.poses)

  nprotein::Int = length(deck.protein)
  nligand::Int = length(deck.ligand)
  nforcefield::Int = length(deck.forcefield)
  nposes::Int = size(deck.poses)[2]

  etotals = CuArray{Float32}(undef, nposes)
  nblocks = ceil(UInt, ceil(nposes / params.ppwi) / params.wgsize)
  nthreads = params.wgsize

  println("Using kernel parameters: <<<$(nblocks),$(nthreads)>>> 1:$nposes")

  shared = false
  shmem_bytes = shared ? sizeof(FFParams) * nforcefield : 0

  # warmup
  @cuda blocks = nblocks threads = nthreads shmem = shmem_bytes fasten_main(
    nprotein,
    nligand,
    nforcefield,
    nposes,
    protein,
    ligand,
    forcefield,
    poses,
    etotals,
    Val(shared),
    Val(convert(Int, params.ppwi)),
  )
  CUDA.synchronize()
  elapsed = @elapsed begin
    for _ = 1:params.iterations
      @cuda blocks = nblocks threads = nthreads shmem = shmem_bytes fasten_main(
        nprotein,
        nligand,
        nforcefield,
        nposes,
        protein,
        ligand,
        forcefield,
        poses,
        etotals,
        Val(shared),
        Val(convert(Int, params.ppwi)),
      )
    end
    CUDA.synchronize()
  end

  (Array{Float32}(etotals), elapsed, params.ppwi)
end

# const NUM_TD_PER_THREAD = 8

@fastmath function fasten_main(
  nprotein::Int,
  nligand::Int,
  nforcefield::Int,
  nposes::Int,
  protein_::CuDeviceVector{Atom},
  ligand_::CuDeviceVector{Atom},
  forcefield_::CuDeviceVector{FFParams},
  poses_::CuDeviceMatrix{Float32},
  etotals::CuDeviceVector{Float32},
  ::Val{Shared},
  ::Val{PPWI},
) where {Shared,PPWI}

  protein = CUDA.Const(protein_)
  ligand = CUDA.Const(ligand_)
  global_forcefield = CUDA.Const(forcefield_)
  poses = CUDA.Const(poses_)

  #  Get index of first TD
  ix = (blockIdx().x - 1) * blockDim().x * PPWI + threadIdx().x

  #  Have extra threads do the last member intead of return.
  #  A return would disable use of barriers, so not using return is better
  ix = ix <= nposes ? ix : nposes - PPWI

  etot = MArray{Tuple{PPWI},Float32}(undef)
  transform = MArray{Tuple{PPWI,3},Vec4f32}(undef)

  if Shared
    forcefield = @cuDynamicSharedMem(FFParams, (nforcefield))
    if ix < nforcefield
      @inbounds forcefield[ix] = global_forcefield[ix]
    end
  else
    forcefield = global_forcefield
  end

  lsz = blockDim().x
  @simd for i = 1:PPWI
    index = (ix) + (i - 1) * lsz
    @inbounds sx::Float32 = sin(poses[1, index])
    @inbounds cx::Float32 = cos(poses[1, index])
    @inbounds sy::Float32 = sin(poses[2, index])
    @inbounds cy::Float32 = cos(poses[2, index])
    @inbounds sz::Float32 = sin(poses[3, index])
    @inbounds cz::Float32 = cos(poses[3, index])
    @inbounds transform[i, 1] = #
      Vec4f32(cy * cz, sx * sy * cz - cx * sz, cx * sy * cz + sx * sz, poses[4, index])
    @inbounds transform[i, 2] = #
      Vec4f32(cy * sz, sx * sy * sz + cx * cz, cx * sy * sz - sx * cz, poses[5, index])
    @inbounds transform[i, 3] = #
      Vec4f32(-sy, sx * cy, cx * cy, poses[6, index])
    @inbounds etot[i] = 0
  end

  if Shared
    CUDA.sync_threads()
  end

  for il = 1:nligand

    @inbounds l_atom::Atom = ligand[il]
    @inbounds l_params::FFParams = forcefield[l_atom.type+1]
    lhphb_ltz::Bool = l_params.hphb < Zero
    lhphb_gtz::Bool = l_params.hphb > Zero

    lpos = MArray{Tuple{PPWI},Vec3f32}(undef)

    @simd for i = 1:PPWI
      # Transform ligand atom
      @inbounds lpos[i] = Vec3f32(
        transform[i, 1].w +
        l_atom.x * transform[i, 1].x +
        l_atom.y * transform[i, 1].y +
        l_atom.z * transform[i, 1].z,
        transform[i, 2].w +
        l_atom.x * transform[i, 2].x +
        l_atom.y * transform[i, 2].y +
        l_atom.z * transform[i, 2].z,
        transform[i, 3].w +
        l_atom.x * transform[i, 3].x +
        l_atom.y * transform[i, 3].y +
        l_atom.z * transform[i, 3].z,
      )
    end

    for ip = 1:nprotein
      @inbounds p_atom = protein[ip]
      @inbounds p_params::FFParams = forcefield[p_atom.type+1]

      radij::Float32 = p_params.radius + l_params.radius
      r_radij::Float32 = One / radij

      elcdst::Float32 = (p_params.hbtype == HbtypeF && l_params.hbtype == HbtypeF) ? Four : Two
      elcdst1::Float32 = (p_params.hbtype == HbtypeF && l_params.hbtype == HbtypeF) ? Quarter : Half
      type_E::Bool = ((p_params.hbtype == HbtypeE || l_params.hbtype == HbtypeE))

      phphb_ltz::Bool = p_params.hphb < Zero
      phphb_gtz::Bool = p_params.hphb > Zero
      phphb_nz::Bool = p_params.hphb != Zero
      p_hphb::Float32 = p_params.hphb * (phphb_ltz && lhphb_gtz ? -One : One)
      l_hphb::Float32 = l_params.hphb * (phphb_gtz && lhphb_ltz ? -One : One)
      distdslv::Float32 =
        (phphb_ltz ? (lhphb_ltz ? Npnpdist : Nppdist) : (lhphb_ltz ? Nppdist : -Float32Max))
      r_distdslv::Float32 = One / distdslv

      chrg_init::Float32 = l_params.elsc * p_params.elsc
      dslv_init::Float32 = p_hphb + l_hphb

      @simd for i = 1:PPWI
        @inbounds x::Float32 = lpos[i].x - p_atom.x
        @inbounds y::Float32 = lpos[i].y - p_atom.y
        @inbounds z::Float32 = lpos[i].z - p_atom.z

        distij::Float32 = sqrt(x ^ 2 + y ^ 2 + z ^ 2)

        # Calculate the sum of the sphere radii
        distbb::Float32 = distij - radij
        zone1::Bool = (distbb < Zero)

        # Calculate steric energy
        @inbounds etot[i] += (One - (distij * r_radij)) * (zone1 ? Two * Hardness : Zero)

        # Calculate formal and dipole charge interactions
        chrg_e::Float32 =
          chrg_init * ((zone1 ? One : (One - distbb * elcdst1)) * (distbb < elcdst ? One : Zero))
        neg_chrg_e::Float32 = -abs(chrg_e)
        chrg_e = type_E ? neg_chrg_e : chrg_e
        @inbounds etot[i] += chrg_e * Cnstnt

        # Calculate the two cases for Nonpolar-Polar repulsive interactions
        coeff::Float32 = (One - (distbb * r_distdslv))
        dslv_e::Float32 = dslv_init * ((distbb < distdslv && phphb_nz) ? One : Zero)
        dslv_e *= (zone1 ? One : coeff)
        @inbounds etot[i] += dslv_e
      end
    end
  end

  td_base = (blockIdx().x - 1) * blockDim().x * PPWI + threadIdx().x
  if td_base <= nposes
    @simd for i = 1:PPWI
      @inbounds etotals[td_base+(i-1)*blockDim().x] = etot[i] * Half
    end
  end

  nothing
end

main()

max_stars_repo_name: wsmoses/Enzyme-BUDE

max_issues_repo_name: wsmoses/Enzyme-BUDE

max_forks_repo_name: UoB-HPC/miniBUDE


text:
"
    Iterating over an AbstractGroup is the same as iterating over the set.
"
abstract type AbstractGroup end

# TODO: convert `Set` to `AbstractSet` where possible to support OrderedSets et al
"""
    Stucture consisting of a set and a binary operation. No constraints are put on either expression.
"""
struct Groupoid{T} <: AbstractGroup
    set::Set{T}
    operation::Function
end

"""
    Group.
    The following axioms must hold. Assuming a group G = (S, ∘), then
      - Closure: the group must be closed under the binary operation.
        ∀x,y ∈ S: x ∘ y ∈ S
      - Associativity:
        ∀x,y,z ∈ S: x ∘ (y ∘ z) =  (x ∘ y) ∘ z
      - Identity:
        ∃e ∈ S: ∀x ∈ S, x ∘ e = x = e ∘ x
      - Inverses:
        ∀x ∈ S, ∃ x⁻¹ ∈ S : x ∘ x⁻¹ = e = x⁻¹ ∘ x
"""
struct Group{T} <: AbstractGroup where {T}
    set::Set{T}
    operation::Function
    Group{T}(set, operation) where T = if _validate_group(Groupoid(set, operation)) new(set, operation) end
end
Group(set::Set{T}, operation::Function) where {T} = Group{T}(set, operation)

"""
    Subgroup of a group.
    Create using SubGroup(group, subset).
    Must be closed under the group operation, include the identity element and include inverses for each element.
"""
struct SubGroup{T} <: AbstractGroup where {T}
    group::Group{T}
    set::Set{T}
    operation::Function
end

function SubGroup(group::Group{T}, subset::Set{<:T})::SubGroup{T} where {T}
    @assert subset ⊆ group.set
    g = Groupoid(subset, group.operation)
    assert_closure(g)
    _assert_inverses(get_identity_element(g), g.set, g.operation)
    return SubGroup(group, subset, group.operation)
end

abstract type AbstractGroupAction end

struct _GroupActionLike{T} <: AbstractGroupAction where {T}
    group::Group{<:T}
    set::Set{<:T}
    action::Function
end

struct GroupAction{T} <: AbstractGroupAction where {T}
    group::Group{<:T}
    set::Set{<:T}
    action::Function
    #GroupAction{T}
end

# TODO: fix convert below, then simplify code as necessary
# import Base: convert
# convert(T::Set{S}, x::Array{S, 1}) where {S} = T(_ for _ in x)
# convert(Set, Array{Int, 1}([1]))
function coset(left::T, set::Set{T}, operation::Function)::Set{T} where {T}
    #Set(x for x in operation.(left, set))
    Set(operation(left, right) for right in set)
end

function coset(set::Set{T}, right::T, operation::Function)::Set{T} where {T}
#    Set(x for x in operation.(set, right))
    Set(operation(left, right) for left in set)

end

function coset(left, group::AbstractGroup)
    coset(left, group.set, group.operation)
end

function coset(group::AbstractGroup, right)
    coset(group.set, right, group.operation)
end

function set_composition(left_set::Set{<:T}, right_set::Set{<:T}, operation::Function)::Set  where {T}
    # `map` is not defined on sets
    # reduce(union, coset(left, right_set, operation) for left in left_set)
    Set(operation(left, right) for left in left_set for right in right_set)
end

function assert_closure(set::Set, closed_in_set::Set, operation::Function)::Nothing
    # We only check left closure.
    @assert set_composition(set, closed_in_set, operation) ⊆ closed_in_set
end

function assert_closure(group::AbstractGroup)::Nothing
    assert_closure(group.set, group.set, group.operation)
end

function assert_closure(ga::AbstractGroupAction)::Nothing
    assert_closure(ga.group.set, ga.set, ga.action)
end

function assert_associativity(set::Set, operation::Function)::Nothing
    for x in set
        for y in set
            for z in set
                a = operation(x, operation(y, z))
                b = operation(operation(x, y), z)
                @assert a == b "$x ∘ ($y ∘ $z) = $a ≠ $b = ($x ∘ $y) ∘ $z"
            end
        end
    end
end

function assert_associativity(group::AbstractGroup)::Nothing
    assert_associativity(group.set, group.operation)
end
function _assert_identity_element(element::T, set::Set{T}, operation::Function)::Nothing where {T}
    # we assume broadcasting works on the operation
    # @assert operation.(element, set) == set  # should be using ordered sets or vectors here
    for x in set
         @assert x == operation(element, x) == operation(x, element)
    end
end

function get_identity_element(group::AbstractGroup) # TODO: return type
    for x in group.set
        try
            _assert_identity_element(x, group.set, group.operation)
            return x
        catch AssertionError
            continue
        end
    end
    throw(AssertionError("No identity element found"))
end

function _assert_inverses(identity_element::T, set::Set{T}, operation::Function)::Nothing where {T}
    for item in set
        @assert identity_element in coset(item, set, operation)
        @assert identity_element in coset(set, item, operation)
    end
end

function _validate_group(group::AbstractGroup)::Bool
    # Validation
    assert_closure(group)
    assert_associativity(group)
    e = get_identity_element(group)  # This raises on failure
    _assert_inverses(e, group.set, group.operation)
    return true
end

function naive_isequal(x::T, y::T) where T
    for f in fieldnames(T)
        if getfield(x, f) != getfield(y, f)
            return false
        end
    end
    return true
end

import Base: ==
==(x::T, y::T) where T<:AbstractGroup = naive_isequal(x, y)
# ==(x::Group, y::Group) = naive_isequal(x, y)
# ==(x::SubGroup, y::SubGroup) = naive_isequal(x, y)
function ==(x::AbstractGroup, y::SubGroup)::Bool
    if y.set == x.set && x.operation == y.operation
        return true
    end
    return false
end
==(x::SubGroup, y::AbstractGroup) = ==(y, x)
# Prevent `MethodError: ==(::SubGroup{Int64}, ::SubGroup{Int64}) is ambiguous`
==(x::SubGroup, y::SubGroup) = naive_isequal(x, y)

import Base: iterate
iterate(group::AbstractGroup) = iterate(group.set)
iterate(group::AbstractGroup, t::T) where {T} = iterate(group.set, t)
import Base: length
length(x::AbstractGroup) = length(x.set)

"""
    N is a _normal_ subgroup of G if ∀g ∈ G, gN = Ng
    N is a normal subgroup of G iff N is a subgroup of G and N is a union of conjugacy classes of G
"""
function isnormal(subgroup::SubGroup)
    for g in subgroup.group
        left_coset = coset(g, subgroup)
        right_coset = coset(subgroup, g)
        if coset(g, subgroup) != coset(subgroup, g)
            return false
        end
    end
    return true
end

function iscyclic(group::AbstractGroup)
    if length(find_generators(group)) > 0
        return true
    end
    return false
end

function quotient_group(subgroup::SubGroup)
    @assert isnormal(subgroup) "Subgroup must be normal"
    # The set here is equal to the "partition" of the group into cosets created by group elements with the subgroup
    f(x::Set, y::Set) = set_composition(x, y, subgroup.operation)
    return Group(Set(coset(g, subgroup) for g in subgroup.group), f)
end

function generate_subgroup(group::AbstractGroup, generator)::SubGroup
    e = get_identity_element(group)
    generated_set = Set([e generator])
    p = group.operation(generator, generator)
    while p != e
        push!(generated_set, p)
        p = group.operation(p, generator)
    end
    SubGroup(group, generated_set, group.operation)
end

function find_generators(group::AbstractGroup)::Set
    filter(g->generate_subgroup(group, g) == group, group.set)
end

function _cayley_table(group::AbstractGroup) # TODO: move to utils
    elements = [x for x in group.set]
    inside = map(x->
            map(y->group.operation(x, y),
            elements
            ), elements)

    return elements, inside
end
# _cayley_table((generate_subgroup(G, 6) |> quotient_group))[2]

"""

Conjugate x by g, i.e. perform gxg⁻¹

"""
function conjugate(group::AbstractGroup, x::T, g::T) where {T}
    group.operation(group.operation(g, x), inv(group, g))
end

"""
Return { gxg⁻¹: ∀g ∈ G}
"""
function conjugacy_class(group::AbstractGroup, x)::Set
    Set(conjugate(group, x, g) for g in group)
end

"""
Return the distinct conjugacy classes in `group`.
The set of distinct conjugacy classes forms a partition of the group.
"""
function conjugacy_classes(group::AbstractGroup)::Set
    Set(conjugacy_class(group, x) for x in group)
end

"""
Obtain the inverse of `x` in `group`.
"""
function inv(group::AbstractGroup, x)
    e = get_identity_element(group)
    for h in group
        if group.operation(x, h) == group.operation(h, x) == e
            return h
        end
    end
    return missing  # TODO: missing or nothing?
end

abstract type AbstractHomomorphism end
"""
An isomorphism ϕ: (G, ∘) → (H, ⋆) is a mapping which satisfies the following properties
- ϕ is one-to-one and onto

See also: Homomorphism
"""
struct Isomorphism <: AbstractHomomorphism
    from_group::AbstractGroup
    to_group::AbstractGroup
    mapping::Dict
end

"""
An isomorphism ϕ: (G, ∘) → (H, ⋆) is a mapping which satisfies the following property:
    ∀x,y ∈ G, ϕ(x ∘ y) = ϕ(x) ⋆ ϕ(y)   (i.e. it preserves composites)
"""
struct Homomorphism <: AbstractHomomorphism
    from_group::AbstractGroup
    to_group::AbstractGroup
    mapping::Dict
end

function order(group::AbstractGroup)
    return length(group.set)
end

function order(group::AbstractGroup, x)
    return length(generate_subgroup(group, x))
end

function get_isomorphism(a::AbstractGroup, b::AbstractGroup)
    if a == b  # Automorphism
         return Dict(k=>k for k in a)
    end
    if iscyclic(a) && iscyclic(b) && order(a) == order(b)  # Cyclic subgroups of the same order
        gen_a, gen_b = find_generators(a) |> first, find_generators(b) |> first
        mapping = Dict(gen_a => gen_b)
        p_a = gen_a
        p_b = gen_b
        for i in 1:order(a)
            p_a = a.operation(p_a, gen_a)
            p_b = b.operation(p_b, gen_b)
            mapping[p_a] = p_b
        end
        return mapping
    end
    throw("Not implemented")
end

"""
Assert the following property, where mapping=ϕ, from=(G, ∘), to=(H, ⋆)
    ∀ x,y ∈ G; ϕ(x, y) = ϕ(x) ⋆ ϕ(y)
"""
function _assert_homomorphism_property(mapping::Dict, from::AbstractGroup, to::AbstractGroup)
   for x in from
        for y in from
            @assert mapping[from.operation(x, y)] == to.operation(mapping[x], mapping[y])
        end
    end
end

"""
Let (G, ∘) be a group, X be a set and ^ be a group action.
    ∀ g,h ∈ G; ∀x in X; g ^ (h ^ x) = (g ∘ h) ^ x
"""
function _assert_homomorphism_property(ga::AbstractGroupAction)
    for g in ga.group
        for h in ga.group
            for x in ga.set
                @assert ga.action(g, ga.action(h, x)) == ga.action(ga.group.operation(g, h), x)
            end
        end
    end
end

function _transform(ϕ::Function, arr::Union{Vector, AbstractSet})::Union{Vector, AbstractSet}
    return ϕ.(arr)
end

function _transform(ϕ::Dict, arr::Union{Vector, AbstractSet})::Union{Vector, AbstractSet}
    return _transform(x-> ϕ[x], arr)
end

"""
Let G be a group and ϕ: G → H be a homomorphism. Then,
    Ker ϕ = {g ∈ G: ϕ(g) = e_H}

Let V, W be vector subspaces and t: V → W be a linear transfomation. Then,
    Ker t = {v⃗ ∈ V: ϕ(v⃗) = 0⃗}
"""
function kernel(transformation::Union{Function, Dict}, from::Union{AbstractSet, AbstractGroup}, identity)
    v = [g for g in from]
    t = _transform(transformation, v)
    mask = t .== identity
    return Set(v[mask])
end

function image(transformation, from::AbstractSet)::Set
   _transform(transformation, from) |> Set
end
function image(transformation, from::AbstractGroup)::Set
   _transform(transformation, from.set) |> Set
end

function GroupAction(group::Group{<:T}, set::Set{<:T}, action::Function)::GroupAction{T} where {T}
    ga = _GroupActionLike(group, set, action)
    assert_closure(ga)

    e = get_identity_element(ga.group)
    for x in ga.set
        @assert ga.action(e, x) == x
    end

    _assert_homomorphism_property(ga)

    return GroupAction{T}(group, set, action)
end

"""
    Orb x = {∀g ∈ G, g ^ x}
"""
function orbit(ga::AbstractGroupAction, x)::Set  # TODO: define type
    @assert x in ga.set
    return ga.action.(ga.group, x) |> Set
end

"""
    Get the set of all orbits for the group action
"""
function orbits(ga::AbstractGroupAction)::Set
    return Set(orbit(ga, x) for x in ga.set)
end

"""
    Stab x = {g ∈ G: g^x = x}
"""
function stabilizer(ga::AbstractGroupAction, x)
    @assert x in ga.set
    v = [g for g in ga.group.set]
    t = ga.action.(v, x)
    stable_mask = t .== x
    return Set(v[stable_mask])
end

# function stabilizers(ga::AbstractGroupAction)
#     return Set(stabilizer(ga, x) for x in ga.set)
# end

"""
    Fix g = {x ∈ X: g^x = x}
"""
function fixed_set(ga::AbstractGroupAction, g)
    @assert g in ga.group.set
    v = [x for x in ga.set]
    t = ga.action.(g, v)
    return Set(v[t .== v])
end


text:
module IRViewer
import InteractiveUtils
using Cascadia
using Gumbo

include("color.jl")
include("ir2graph.jl")
include("macro.jl")
export @llvmir2svg,@llvmir2html

end # module

max_stars_repo_name: ChenNingCong/IRViewer

max_issues_repo_name: ChenNingCong/IRViewer

max_forks_repo_name: ChenNingCong/IRViewer


text:
# Pseudea.jl, daughter of Eris

export Sn,Zn,Cu,S
export randCZTS,randCZTSlattice,moduloCZTSlattice

const Sn=-3
const Zn=-2
const Cu=-1
const S=+2

randCZTS() = [Cu,Sn,Zn,Zn,S,S,S,S][rand(1:8)]
randCZTSlattice(n::Int) =  [randCZTS() for i in 1:n, j in 1:n, k in 1:n]

function moduloCZTSlattice(n::Int)
    lattice=Array{Int}(n,n,n)

    for i in 1:n
        for j in 1:n
            for k in 1:n
                if ((i+j+k)%2==0) lattice[i,j,k]=+2 end # paint the sulphur in; i+j+k = EVEN
                if ((i+j+k)%4==1) lattice[i,j,k]=-2 end
                if ((i+j+k)%8==3) lattice[i,j,k]=-3 end
                if ((i+j+k)%8==7) lattice[i,j,k]=-1 end
            end
        end
    end
    return(lattice)
end


text:
function setup_etc(backend::GPUBackend, layer::RandomNormalLayer)
  cuda_rand_states = CuPtr
  kernel = backend.mocha.stdnormal_init
  rnd_state_size_blob = make_blob(backend, Float64, 1, 1, 1, 1)
  CUDA.launch(backend.mocha.stdnormal_alloc_size, 1, 1, (rnd_state_size_blob.ptr.p, ))
  rnd_state_size = Float64[0]
  copy!(rnd_state_size, rnd_state_size_blob)
  destroy(rnd_state_size_blob)
  rnd_state_size = int(rnd_state_size[1])

  etc = Any[]
  outlen = prod(layer.output_dims)
  for i = 1:length(layer.tops)
      len = outlen*layer.batch_sizes[i]
      cuda_rand_states = CUDA.cualloc(UInt8, rnd_state_size*len)
      x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
      seed = rand(UInt)
      println("launching stdnormal init on bloc $i with len $len")
      CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X, (cuda_rand_states, seed))
      push!(etc, cuda_rand_states)
  end
  return etc
end

function destroy_etc(backend::GPUBackend, state::RandomNormalLayerState)
    for i = 1:length(state.etc)
        CUDA.free(state.etc[i])
    end
end

function forward(backend::GPUBackend, state::RandomNormalLayerState, inputs::Vector{Blob})
    for i = 1:length(state.blobs)
        len = length(state.blobs[i])
        x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
        data_type = state.layer.eltype
        if data_type == Float32
            kernel = backend.mocha.stdnormal_forward_float
        elseif data_type == Float64
            kernel = backend.mocha.stdnormal_forward_double
        end

        CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X,
                    (state.etc[i], state.blobs[i].ptr.p, len,))
    end
end

function backward(backend::GPUBackend, state::RandomNormalLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})

end


text:
using Flux
using Plots
using DifferentialEquations

ω0 = 0.5
γ = 0.1
x₀ = 10.0
v₀ = 0.0
tspan = (0.0, 20.0)
acceleration(dx, x, p, t) = -ω0^2*x-2*γ*dx
prob = SecondOrderODEProblem(acceleration, v₀, x₀, tspan)
sol = solve(prob)

plot(sol,label=["Velocity" "Position"])

# generate the training dataset
plot_t = 0:0.01:10
data_plot = sol(plot_t)
positions_plot = [state[2] for state in data_plot]
force_plot = [force(state[1],state[2],k,t) for state in data_plot]

t = 0:3.3:10
dataset = sol(t)
position_data = [state[2] for state in sol(t)]
force_data = [force(state[1],state[2],k,t) for state in sol(t)]

plot(plot_t,force_plot,xlabel="t",label="True Force")
scatter!(t,force_data,label="Force Measurements")

# create NN
# create NN
NNForce = Chain(x -> [x],
                Dense(1,32,tanh),
                Dense(32,1),
                first)

loss() = sum(abs2,NNForce(position_data[i]) - force_data[i] for i in 1:length(position_data))
loss()

opt = Flux.Descent(0.01)
data = Iterators.repeated((), 5000)
iter = 0
cb = function () #callback function to observe training
    global iter += 1
    if iter % 500 == 0
        display(loss())
    end
end
display(loss())

# train the NN
Flux.train!(loss, Flux.params(NNForce), data, opt; cb=cb)

# visualize the result
learned_force_plot = NNForce.(positions_plot)
plot(plot_t,force_plot,xlabel="t",label="True Force")
plot!(plot_t,learned_force_plot,label="Predicted Force")
scatter!(t,force_data,label="Force Measurements")

# solve SHO problem to justify PINN
acceleration2(dx,x,k,t) = -ω0^2*x
prob_simplified = SecondOrderODEProblem(acceleration2, v₀, x₀, tspan)
sol_simplified = solve(prob_simplified)
p2 = plot(sol_simplified, label=["Velocity" "Position"])
plot!(p2, sol_simplified,label=["Velocity Simplified" "Position Simplified"])

# create Physics Loss Term to force solution to satisfy SHO ode as well
random_positions = [2rand()-1 for i in 1:100] # random values in [-1,1]
loss_ode() = sum(abs2,NNForce(x) - (-(ω0^2)*x) for x in random_positions)
loss_ode()

λ = 0.025  # factor controlling physics regularization
composed_loss() = loss() + λ*loss_ode()

opt = Flux.Descent(0.01)
data = Iterators.repeated((), 5000)
iter = 0
cb = function () #callback function to observe training
    global iter += 1
    if iter % 500 == 0
        display(composed_loss())
    end
end
display(composed_loss())

# train the new model
Flux.train!(composed_loss, Flux.params(NNForce), data, opt; cb=cb)

loss()
λ * loss_ode()

# visualize the final result!
learned_force_plot = NNForce.(positions_plot)
plot(plot_t,force_plot,xlabel="t",label="True Force")
plot!(plot_t,learned_force_plot,label="Predicted Force")
scatter!(t,force_data,label="Force Measurements")


text:
module NormalVectorsTests

using Test
using Gridap

tags = [23,24,25]
model = CartesianDiscreteModel(partition=(3,4,2))

btrian = BoundaryTriangulation(model,tags)

n = NormalVector(btrian)

# @santiagobadia : Francesc can you take a look at this error?
# sr = "NormalVector object:\n length: 20"
# sprint(show,"text/plain",n)
# @test sprint(show,"text/plain",n) == sr

# @test string(n) == "NormalVector object"

bquad = CellQuadrature(btrian,degree=2)
q = coordinates(bquad)

n_q = evaluate(n,q)
a = collect(n_q)

@test a[1] ≈ VectorValue{3,Float64}[(3.1225e-17, -1.0, -0.0), (0.0, -1.0, -1.11022e-16), (4.16334e-17, -1.0, -0.0), (-0.0, -1.0, -0.0)]
@test a[14] ≈ VectorValue{3,Float64}[(-1.04083e-17, -1.0, -1.38778e-17), (-0.0, -1.0, -0.0), (-0.0, -1.0, -0.0), (-0.0, -1.0, -0.0)]

bphi = CellGeomap(btrian)

x = evaluate(bphi,q)

#writevtk(btrian,"btrian")
#writevtk(x,"x",pointdata=["n"=>n_q])

strian = SkeletonTriangulation(model,"interior")
squad = CellQuadrature(strian,degree=2)
qs = coordinates(squad)

sphi = CellGeomap(strian)

n = NormalVector(strian)
x = evaluate(sphi,qs)
n_q = evaluate(n,qs)

#writevtk(strian,"strian")
#writevtk(x,"x",pointdata=["n"=>n_q])

end # module


text:
scripts_dir = @__DIR__
reports_dir = (length(ARGS) > 0 ? ARGS[1] : "examples/benchmarking")

nbfile = joinpath(reports_dir, Libc.strftime("benchmark_runs_%y%m%d.csv", time()))

if isfile(nbfile)
    @info("Deleting existing benchmark file: ", nbfile)
    rm(nbfile)
end

nreps = (length(ARGS) > 0 ? parse(Int, ARGS[2]) : 10)

cmd = `$(Sys.BINDIR)/julia --color=yes -L $scripts_dir/activate_pwd_env.jl -- $scripts_dir/compare_optimizers.jl update --benchmarkfile $nbfile --numreps $nreps --optimizers all`

run(cmd)


text:
function packages()
    if RDatasets.__packages === nothing
        path = joinpath(dirname(@__FILE__), "..", "doc", "packages.csv")
        global __packages = CSV.read(path)
    end
    return RDatasets.__packages::DataFrame
end

max_stars_repo_name: SCXFPROTOCOL/GameRewards-System

max_issues_repo_name: SCXFPROTOCOL/GameRewards-System

max_forks_repo_name: SCXFPROTOCOL/GameRewards-System


text:
function add_node!(bg::BondGraph, nodes)
    for node in nodes
        add_node!(bg, node)
    end
end

function add_node!(bg::BondGraph, node::AbstractNode)
    g.add_vertex!(bg, node) || @warn "$node already in model"
end

function remove_node!(bg::BondGraph, nodes)
    for node in nodes
        remove_node!(bg, node)
    end
end

function remove_node!(bg::BondGraph, node::AbstractNode)
    g.rem_vertex!(bg, node) || @warn "$node not in model"
    for bond in filter(bond -> node in bond, bg.bonds)
        g.rem_edge!(bg, srcnode(bond), dstnode(bond))
    end
end

function connect!(bg::BondGraph, srcnode::AbstractNode, dstnode::AbstractNode;
        srcportindex=nextfreeport(srcnode), dstportindex=nextfreeport(dstnode))
    srcnode in bg.nodes || error("$srcnode not found in bond graph")
    dstnode in bg.nodes || error("$dstnode not found in bond graph")
    srcport = Port(srcnode, srcportindex)
    dstport = Port(dstnode, dstportindex)
    return g.add_edge!(bg, srcport, dstport)
end

function disconnect!(bg::BondGraph, node1::AbstractNode, node2::AbstractNode)
    # rem_edge! removes the bond regardless of the direction of the bond
    return g.rem_edge!(bg, node1, node2)
end

# TODO
# Flip bond function

function swap!(bg::BondGraph, oldnode::AbstractNode, newnode::AbstractNode)
    _check_port_number(oldnode,newnode)

    # may be a redundant check
    if !g.has_vertex(bg, newnode)
        add_node!(bg, newnode)
    end

    srcnodes = g.inneighbors(bg, oldnode)
    dstnodes = g.outneighbors(bg, oldnode)
    remove_node!(bg, oldnode)

    for src in srcnodes
        connect!(bg, src, newnode)
    end
    for dst in dstnodes
        connect!(bg, newnode, dst)
    end
end

_check_port_number(oldnode::AbstractNode, newnode::AbstractNode) =
    numports(newnode) >= numports(oldnode) || error("New node must have a greater or equal number of ports to the old node")
_check_port_number(oldnode::AbstractNode, newnode::Junction) = true

# TODO implement according to https://bondgraphtools.readthedocs.io/en/latest/api.html#BondGraphTools.expose
# function expose!()

# end

# Inserts an AbstractNode between two connected (bonded) nodes
# The direction of the original bond is preserved by this action
function insert_node!(bg::BondGraph, bond::Bond, newnode::AbstractNode)
    src = srcnode(bond)
    dst = dstnode(bond)

    disconnect!(bg, src, dst)

    try
        add_node!(bg, newnode)
        connect!(bg, src, newnode)
        connect!(bg, newnode, dst)
    catch e
        # if connection fails, reconnect original bond
        disconnect!(bg, src, newnode)
        disconnect!(bg, newnode, dst)
        connect!(bg, src, dst)
        error(e)
    end
end
function insert_node!(bg::BondGraph, tuple::Tuple, newnode::AbstractNode)
    bonds = getbonds(bg, tuple)
    isempty(bonds) && error("$(tuple[1]) and $(tuple[2]) are not connected")
    insert_node!(bg, bonds[1], newnode)
end

function merge_nodes!(bg::BondGraph, node1::AbstractNode, node2::AbstractNode; junction=EqualEffort())
    node1.type == node2.type || error("$(node1.name) must be the same type as $(node2.name)")

    # node1 taken as the node to keep
    for nb in g.all_neighbors(bg, node1)
        junc = deepcopy(junction)
        bond = getbonds(bg, node1, nb)[1]
        insert_node!(bg, bond, junc)
        swap!(bg, node2, junc)
    end
end
function merge_nodes!(bg::BondGraph, node1::Junction, node2::Junction)
    # node1 taken as the node to keep
    # remove conflicting connections between junctions if they exist
    disconnect!(bg, node1, node2)
    shared_neighbors = intersect(g.all_neighbors(bg, node1), g.all_neighbors(bg, node2))
    for shared_neighbor in shared_neighbors
        disconnect!(bg, node2, shared_neighbor)
    end
    swap!(bg, node2, node1)
end

function simplify_junctions!(bg::BondGraph; remove_redundant=true, squash_identical=true)
    junctions = filter(n -> n isa Junction, bg.nodes)

    # Removes junctions with 2 or less connected ports
    if remove_redundant
        for j in junctions
            n_nbrs = length(g.all_neighbors(bg, j))
            if n_nbrs == 2
                #srcnode = g.inneighbors(bg, j)[1]
                #dstnode = g.outneighbors(bg, j)[1]
                node1, node2 = g.all_neighbors(bg, j)
                remove_node!(bg, j)
                # bond direction may not be preserved here
                connect!(bg, node1, node2)
            elseif n_nbrs < 2
                remove_node!(bg, j)
            end
        end
    end

    # Squashes identical copies of the same junction type into one junction
    if squash_identical
        for j in junctions, nbr in g.all_neighbors(bg, j)
            g.has_vertex(bg, j) || continue # in case j was removed
            if type(j) == type(nbr)
                merge_nodes!(bg, j, nbr)
            end
        end
    end
    bg
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
@doc """
    struct LaserGeometry

The geometrical aspects of the laser can be described through its orientation
with respect to a given reference system and by its propagation.
"""
LaserGeometry

struct LaserGeometry{O,R}
    orientation::O
    rotation_matrix::R
end

function LaserGeometry(orientation; isversor=false)
    v₁ = !isversor ? normalize(orientation[1]) : orientation[1]
    v₃ = !isversor ? normalize(orientation[2]) : orientation[2]
    v₂ = v₃ × v₁

    # original basis
    e₁ = SVector{3}(1,0,0)
    e₂ = SVector{3}(0,1,0)
    e₃ = SVector{3}(0,0,1)

    # original basis matrix
    E = [e₁ e₂ e₃]

    # desired direction matrix
    V = [v₁ v₂ v₃]

    # rotation matrix
    R = E * transpose(V)

    LaserGeometry(orientation, R)
end

function sym2vec(dir::Symbol)
    if dir == :x
        SVector{3}(1,0,0)
    elseif dir == :y
        SVector{3}(0,1,0)
    else
        SVector{3}(0,0,1)
    end
end

function LaserGeometry(orientation::Tuple{Symbol,Symbol})
    if orientation[1] == :x && orientation[2] == :z
        return LaserGeometry((:x, :z), I)
    end
    v₁ = sym2vec(orientation[1])
    v₃ = sym2vec(orientation[2])

    LaserGeometry((v₁, v₃), isversor=true)
end


text:
"""
    lower, upper, xmin, xmax = fixedrange!(X[, lower, upper, xmin, xmax; obsdim, operate_on])

or

    lower, upper, xmin, xmax = fixedrange!(D[, lower, upper, xmin, xmax; operate_on])

where `X` is of type Matrix or Vector and `D` of type DataFrame.
Normalize `X` or `D` along `obsdim` to the interval [lower:upper].
If `lower` and `upper`  are omitted the default range is [0:1].

`lower`     :  (Scalar) Lower limit of new range.
               Defaults to 0.

`upper`     :  (Scalar) Upper limit of new range.
               Defaults to 1.

`xmin`      :  (Vector) Minimum values of data before normalization. `xmin` will
               correspond to `lower` after transformation.
               Defaults to `minimum(X, obsdim)`.

`xmin`      :  (Vector) Maximum value of data before normalization. `xmax` will
               correspond to `upper` after transformation.
               Defaults to `maximum(X, obsdim)`.

`obsdim`    :  Specify which axis corresponds to observations.
               Defaults to obsdim=2 (observations are columns of matrix)
               For DataFrames `obsdim` is obsolete and rescaling occurs
               column wise.

`operate_on`:  Specify the indices of columns or rows to be centered.
               Defaults to all columns/rows.
               For DataFrames this must be a vector of symbols, not indices.
               E.g. `operate_on`=[1,3] will perform centering on columns
               with index 1 and 3 only (if obsdim=1, else rows 1 and 3)

Note on DataFrames:
Columns containing `missing` values are skipped.
Columns containing non numeric elements are skipped.

Examples:

    X = rand(4, 100)
    x = rand(10)
    D = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    lower, upper, xmin, xmax = fixedrange!(X)
    lower, upper, xmin, xmax = fixedrange!(X, -1, 1)
    lower, upper, xmin, xmax = fixedrange!(X, -1, 1, obsdim=1)
    lower, upper, xmin, xmax = fixedrange!(X, -1, 1, obsdim=1, operate_on=[1,4])

    lower, upper, xmin, xmax = fixedrange!(D)
    lower, upper, xmin, xmax = fixedrange!(D, -1, 1)
    lower, upper, xmin, xmax = fixedrange!(D, -1, 1, operate_on=[:A,:B])
"""
function fixedrange!(X; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    fixedrange!(X, convert(ObsDimension, obsdim), operate_on)
end

function fixedrange!(X::AbstractArray{T,N}, ::ObsDim.Last, operate_on) where {T,N}
    fixedrange!(X, ObsDim.Constant{N}(), operate_on)
end

function fixedrange!(X, obsdim::ObsDim.Constant{M}, operate_on) where {M}
    lower = 0
    upper = 1
    xmin = minimum(X; dims=M)[operate_on]
    xmax = maximum(X; dims=M)[operate_on]
    fixedrange!(X, lower, upper, xmin, xmax, obsdim, operate_on)
end

function fixedrange!(X, lower, upper; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    fixedrange!(X, lower, upper, convert(ObsDimension, obsdim), operate_on)
end

function fixedrange!(X, lower, upper, obsdim::ObsDim.Constant{M}, operate_on) where {M}
    xmin = minimum(X; dims=M)[operate_on]
    xmax = maximum(X; dims=M)[operate_on]
    fixedrange!(X, lower, upper, xmin, xmax, obsdim, operate_on)
end

function fixedrange!(X::AbstractArray{T,M}, lower, upper, obsdim::ObsDim.Last, operate_on) where {T,M}
    fixedrange!(X, lower, upper, ObsDim.Constant{M}(), operate_on)
end

function fixedrange!(X, lower, upper, xmin, xmax; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    fixedrange!(X, lower, upper, xmin, xmax, convert(ObsDimension, obsdim), operate_on)
end

function fixedrange!(X::AbstractMatrix, lower::Real, upper::Real, xmin::AbstractVector, xmax::AbstractVector, ::ObsDim.Constant{1}, operate_on::AbstractVector)
    @assert length(xmin) == length(xmax) == length(operate_on)
    xrange = xmax .- xmin
    scale = upper - lower
    nObs, nVars = size(X)

    for (i, iVar) in enumerate(operate_on)
        @inbounds for iObs in 1:nObs
            X[iObs, iVar] = lower + (X[iObs, iVar] - xmin[i]) / xrange[i] * scale
        end
    end
    lower, upper, xmin, xmax
end

function fixedrange!(X::AbstractMatrix, lower::Real, upper::Real, xmin::AbstractVector, xmax::AbstractVector, ::ObsDim.Constant{2}, operate_on::AbstractVector)
    @assert length(xmin) == length(xmax) == length(operate_on)
    xrange = xmax .- xmin
    scale = upper - lower
    nVars, nObs = size(X)

    for iObs in 1:nObs
        @inbounds for (i, iVar) in enumerate(operate_on)
            X[iVar, iObs] = lower + (X[iVar, iObs] - xmin[i]) / xrange[i] * scale
        end
    end
    lower, upper, xmin, xmax
end

function fixedrange!(X::AbstractMatrix, lower::Real, upper::Real, xmin::AbstractVector, xmax::AbstractVector, ::ObsDim.Last, operate_on::AbstractVector)
    fixedrange!(X, lower, upper, xmin, xmax, ObsDim.Constant{2}(), operate_on)
end

function fixedrange!(x::AbstractVector; lower=0.0, upper=1.0)
    xmin = minimum(x)
    xmax = maximum(x)
    fixedrange!(x, lower, upper, xmin, xmax)
end

function fixedrange!(x::AbstractVector, lower::Real, upper::Real)
    xmin = minimum(x)
    xmax = maximum(x)
    fixedrange!(x, lower, upper, xmin, xmax)
end

function fixedrange!(x::AbstractVector, lower::Real, upper::Real, xmin::AbstractVector, xmax::AbstractVector, ::ObsDim.Constant{M}, operate_on::AbstractVector) where {M}
    @assert length(xmin) == length(xmax) == length(operate_on)
    xrange = xmax .- xmin
    scale = upper - lower
    nVars = length(x)
    @inbounds for (i, iVar) in enumerate(operate_on)
        x[iVar] = lower + (x[iVar] - xmin[i]) / xrange[i] * scale
    end
    lower, upper, xmin, xmax
end

function fixedrange!(x::AbstractVector, lower::Real, upper::Real, xmin::AbstractVector, xmax::AbstractVector, ::ObsDim.Last, operate_on::AbstractVector)
    fixedrange!(x, lower, upper, xmin, xmax, ObsDim.Constant{1}(), operate_on)
end

function fixedrange!(x::AbstractVector, lower::Real, upper::Real, xmin::Real, xmax::Real)
    xrange = xmax - xmin
    scale = upper - lower
    n = length(x)
    @inbounds for i in 1:n
        x[i] = lower + (x[i] - xmin) / xrange * scale
    end
    lower, upper, xmin, xmax
end

# --------------------------------------------------------------------

function fixedrange!(D::AbstractDataFrame; operate_on=default_scaleselection(D))
    fixedrange!(D, 0, 1, operate_on)
end

function fixedrange!(D::AbstractDataFrame, lower, upper; operate_on=default_scaleselection(D))
    fixedrange!(D, lower, upper, operate_on)
end

function fixedrange!(D::AbstractDataFrame, lower::Real, upper::Real, operate_on::AbstractArray)
    xmin = Float64[]
    xmax = Float64[]

    for colname in operate_on
        if eltype(D[colname]) <: Union{Real,Missing}
            minval = minimum(D[colname])
            maxval = maximum(D[colname])
            if ismissing(minval)
                @warn("Skipping \"$colname\" because it contains missing values")
                continue
            end
            fixedrange!(D, lower, upper, minval, maxval, colname)
            push!(xmin, minval)
            push!(xmax, maxval)
        else
            @warn("Skipping \"$colname\" because data is not of type T <: Real.")
        end
    end
    lower, upper, xmin, xmax
end

function fixedrange!(D::AbstractDataFrame, lower, upper, xmin, xmax; operate_on=default_scaleselection(D))
    fixedrange!(D, lower, upper, xmin, xmax, operate_on)
end

function fixedrange!(D::AbstractDataFrame, lower::Real, upper::Real, xmin::AbstractArray, xmax::AbstractArray, operate_on::AbstractVector)
    @assert length(xmin) == length(xmax) == length(operate_on)
    for (iVar, colname) in enumerate(operate_on)
        fixedrange!(D, lower, upper, xmin[iVar], xmax[iVar], colname)
    end
    lower, upper, xmin, xmax, operate_on
end

function fixedrange!(D::AbstractDataFrame, lower::Real, upper::Real, xmin::Real, xmax::Real, colname::Symbol)
    if any(ismissing, D[colname]) | !(eltype(D[colname]) <: Union{Real,Missing})
        @warn("Skipping \"$colname\" because it contains missing values or is not of type <: Real")
    else
        newcol::Vector{Float64} = convert(Vector{Float64}, D[colname])
        fixedrange!(newcol, lower, upper, xmin, xmax)
        D[colname] = newcol
    end
    lower, upper, xmin, xmax, colname
end

"""
`FixedRangeScaler` is used with the functions `fit()`, `transform()` and `fit_transform()`
to scale data in a Matrix `X` or DataFrame to a fixed range [lower:upper].
After fitting a `FixedRangeScaler` to one data set, it can be used to perform the same
transformation to a new set of data. E.g. fit the `FixedRangeScaler` to your training
data and then apply the scaling to the test data at a later stage. (See examples below).

    fit(FixedRangeScaler, X[, lower, upper; obsdim, operate_on])

    fit_transform(FixedRangeScaler, X[, lower, upper; obsdim, operate_on])

`X`         :  Data of type Matrix or `DataFrame`.

`lower`     :  (Scalar) Lower limit of new range.
               Defaults to 0.

`upper`     :  (Scalar) Upper limit of new range.
               Defaults to 1.

`obsdim`    :  Specify which axis corresponds to observations.
               Defaults to obsdim=2 (observations are columns of matrix)
               For DataFrames `obsdim` is obsolete and rescaling occurs
               column wise.

`operate_on`:  Specify the indices of columns or rows to be centered.
               Defaults to all columns/rows.
               For DataFrames this must be a vector of symbols, not indices.
               E.g. `operate_on`=[1,3] will perform centering on columns
               with index 1 and 3 only (if obsdim=1, else rows 1 and 3)

Note on DataFrames:
Columns containing `missing` values are skipped.
Columns containing non numeric elements are skipped.

Examples:

    Xtrain = rand(100, 4)
    Xtest  = rand(10, 4)
    x = rand(10)
    D = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    scaler = fit(FixedRangeScaler, Xtrain)
    scaler = fit(FixedRangeScaler, Xtrain, -1, 1)
    scaler = fit(FixedRangeScaler, Xtrain, -1, 1, obsdim=1)
    scaler = fit(FixedRangeScaler, Xtrain, -1, 1, obsdim=1, operate_on=[1,3])
    scaler = fit(FixedRangeScaler, D, -1, 1, operate_on=[:A,:B])

    Xscaled = transform(Xtest, scaler)
    transform!(Xtest, scaler)

    Xscaled, scaler = fit_transform(FixedRangeScaler, X, -1, 1, obsdim=1, operate_on=[1,2,4])
    scaler = fit_transform!(FixedRangeScaler, X, -1, 1, obsdim=1, operate_on=[1,2,4])

Note that for `transform!` the data matrix `X` has to be of type <: AbstractFloat
as the scaling occurs inplace. (E.g. cannot be of type Matrix{Int64}). This is not
the case for `transform` however.
For `DataFrames` `transform!` can be used on columns of type <: Integer.
"""
struct FixedRangeScaler{T<:Real,U<:Real,V<:Real,W<:Real,M,I}
    lower::T
    upper::U
    xmin::Vector{V}
    xmax::Vector{W}
    obsdim::ObsDim.Constant{M}
    operate_on::Vector{I}
end

function FixedRangeScaler(X::AbstractArray{T,N}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    FixedRangeScaler(X, convert(ObsDimension, obsdim), operate_on)
end

function FixedRangeScaler(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}, operate_on) where {T<:Real,N,M}
    xmin = vec(minimum(X; dims=M))[operate_on]
    xmax = vec(maximum(X; dims=M))[operate_on]
    FixedRangeScaler(0, 1, xmin, xmax, obsdim, operate_on)
end

function FixedRangeScaler(X::AbstractArray{T,N}, ::ObsDim.Last, operate_on) where {T<:Real,N}
    FixedRangeScaler(X, ObsDim.Constant{N}(), operate_on)
end

function FixedRangeScaler(X::AbstractArray{T,N}, lower, upper; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    FixedRangeScaler(X, lower, upper, convert(ObsDimension, obsdim), operate_on)
end

function FixedRangeScaler(X::AbstractArray{T,N}, lower, upper, obsdim::ObsDim.Constant{M}, operate_on) where {T<:Real,N,M}
    xmin = vec(minimum(X; dims=M))[operate_on]
    xmax = vec(maximum(X; dims=M))[operate_on]
    FixedRangeScaler(lower, upper, xmin, xmax, obsdim, operate_on)
end

function FixedRangeScaler(X::AbstractArray{T,N}, lower, upper, ::ObsDim.Last, operate_on) where {T<:Real,N}
    FixedRangeScaler(X, lower, upper, ObsDim.Constant{N}(), operate_on)
end

function FixedRangeScaler(D::AbstractDataFrame; operate_on=default_scaleselection(D))
    FixedRangeScaler(D, 0, 1, operate_on)
end

function FixedRangeScaler(D::AbstractDataFrame, lower::Real, upper::Real; operate_on=default_scaleselection(D))
    FixedRangeScaler(D, lower, upper, operate_on)
end

function FixedRangeScaler(D::AbstractDataFrame, lower::Real, upper::Real, operate_on::AbstractVector{Symbol})
    xmin = Float64[]
    xmax = Float64[]
    colnames = valid_columns(D, operate_on)
    for colname in colnames
        push!(xmin, minimum(D[colname]))
        push!(xmax, maximum(D[colname]))
    end
    FixedRangeScaler(lower, upper, xmin, xmax, ObsDim.Constant{1}(), colnames)
end

function StatsBase.fit(::Type{FixedRangeScaler}, X::AbstractArray{T,N}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    FixedRangeScaler(X, convert(ObsDimension, obsdim), operate_on)
end

function fit_transform(::Type{FixedRangeScaler}, X::AbstractArray{T,N}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    scaler = FixedRangeScaler(X, convert(ObsDimension, obsdim), operate_on)
    Xnew = transform(X, scaler)
    return Xnew, scaler
end

function fit_transform!(::Type{FixedRangeScaler}, X::AbstractArray{T,N}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    scaler = FixedRangeScaler(X, convert(ObsDimension, obsdim), operate_on)
    transform!(X, scaler)
    return scaler
end

function StatsBase.fit(::Type{FixedRangeScaler}, X::AbstractArray{T,N}, lower, upper; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    FixedRangeScaler(X, lower, upper, convert(ObsDimension, obsdim), operate_on)
end

function fit_transform(::Type{FixedRangeScaler}, X::AbstractArray{T,N}, lower, upper; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    scaler = FixedRangeScaler(X, lower, upper, convert(ObsDimension, obsdim), operate_on)
    Xnew = transform(X, scaler)
    return Xnew, scaler
end

function fit_transform!(::Type{FixedRangeScaler}, X::AbstractArray{T,N}, lower, upper; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,N}
    scaler = FixedRangeScaler(X, lower, upper, convert(ObsDimension, obsdim), operate_on)
    transform!(X, scaler)
    return scaler
end

function StatsBase.fit(::Type{FixedRangeScaler}, D::AbstractDataFrame; operate_on=default_scaleselection(D))
    FixedRangeScaler(D, 0, 1, operate_on)
end

function fit_transform(::Type{FixedRangeScaler}, D::AbstractDataFrame; operate_on=default_scaleselection(D))
    scaler = FixedRangeScaler(D, 0, 1, operate_on)
    Dnew = transform(D, scaler)
    return Dnew, scaler
end

function fit_transform!(::Type{FixedRangeScaler}, D::AbstractDataFrame; operate_on=default_scaleselection(D))
    scaler = FixedRangeScaler(D, 0, 1, operate_on)
    transform!(D, scaler)
    return scaler
end

function StatsBase.fit(::Type{FixedRangeScaler}, D::AbstractDataFrame, lower, upper; operate_on=default_scaleselection(D))
    FixedRangeScaler(D, lower, upper, operate_on)
end

function fit_transform(::Type{FixedRangeScaler}, D::AbstractDataFrame, lower, upper; operate_on=default_scaleselection(D))
    scaler = FixedRangeScaler(D, lower, upper, operate_on)
    Dnew = transform(D, scaler)
    return Dnew, scaler
end

function fit_transform!(::Type{FixedRangeScaler}, D::AbstractDataFrame, lower, upper; operate_on=default_scaleselection(D))
    scaler = FixedRangeScaler(D, lower, upper, operate_on)
    transform!(D, scaler)
    return scaler
end

function transform!(X::AbstractArray{T,N}, cs::FixedRangeScaler) where {T<:AbstractFloat,N}
    fixedrange!(X, cs.lower, cs.upper, cs.xmin, cs.xmax, cs.obsdim, cs.operate_on)
end

function transform!(x::AbstractVector{T}, cs::FixedRangeScaler) where {T<:AbstractFloat}
    fixedrange!(x, cs.lower, cs.upper, cs.xmin, cs.xmax, cs.obsdim, cs.operate_on)
end

function transform!(D::AbstractDataFrame, cs::FixedRangeScaler)
    fixedrange!(D, cs.lower, cs.upper, cs.xmin, cs.xmax, cs.operate_on)
end

function transform(X::AbstractArray{T,N}, cs::FixedRangeScaler) where {T<:AbstractFloat,N}
    Xnew = copy(X)
    fixedrange!(Xnew, cs.lower, cs.upper, cs.xmin, cs.xmax, cs.obsdim, cs.operate_on)
    Xnew
end

function transform(X::AbstractArray{T,N}, cs::FixedRangeScaler) where {T<:Real,N}
    Xnew = convert(AbstractArray{Float64, N}, X)
    fixedrange!(Xnew, cs.lower, cs.upper, cs.xmin, cs.xmax, cs.obsdim, cs.operate_on)
    Xnew
end

function transform(D::AbstractDataFrame, cs::FixedRangeScaler)
    Dnew = deepcopy(D)
    fixedrange!(Dnew, cs.lower, cs.upper, cs.xmin, cs.xmax, cs.operate_on)
    Dnew
end


text:
@pywrap struct IB
end

IB() = IB(ib_insync.IB())
Base.isopen(ib::IB) = ib.isConnected()
function Base.open(ib::IB, args...; kwargs...)
    ib.connect(args...; kwargs...)
    return ib
end
function Base.close(ib::IB)
    isopen(ib) && ib.disconnect()
    return nothing
end

####
#### IBAccount
####

# TODO currency field?
# TODO portfolio margin doesn't have SMA? ELV?
struct IBAccount{C<:CashQuantity} <: Account
    "Account name or number"
    name::String

    "ELV – Initial Margin. Liquidation occurs if < 0 at any time."
    excess_liquidity::C
    "Maximum of (Prior Day SMA +/- Change in Day's Cash +/- Today's Trades Initial Margin Requirements) and
    (ELV - Reg T Margin). Liquidation occurs if < 0 at end of trading day (1550 ET)."
    sma::C
    "Funds available for trading. Equal to ELV – Initial Margin."
    available_funds::C

    "GPV: Absolute value of Long Stock Value + Short Stock Value + Long Option Value + Short Option Value + Long SSF Notional Value + Short SSF Notional Value + Fund Value."
    gross_position_value::C
    "ELV: Equals cash + stock value + bond value + mutual fund value + European and Asian options value (excludes market value U.S. securities & futures options and cash maintained in futures segment)."
    equity_with_loan_value::C
    "NLV: For a securities account equals total cash value + stock value + securities options value + bond value + fund value."
    net_liquidation_value::C
    "Total cash balance. Negative value represented borrowed money."
    total_cash_value::C
end

function IBAccount(ib::IB, account::String)
    sum = ib.accountSummary(account)
    function get(fname, tag)
        idx = findfirst(x -> x[2] == tag, sum)
        idx == nothing && error("Unknown tag $tag")
        _, tag, value, currency, model_code = sum[idx]
        T = fieldtype(IBAccount, fname)
        if T <: CashQuantity
            @assert currency != ""
            value = parse(Float64, value) * parse_currency(currency)
        elseif T <: Number
            value = parse(T, value)
        end
        return value
    end
    IBAccount(
        account,
        get(:excess_liquidity, "ExcessLiquidity"),
        get(:sma, "SMA"),
        get(:available_funds, "AvailableFunds"),
        get(:gross_position_value, "GrossPositionValue"),
        get(:equity_with_loan_value, "EquityWithLoanValue"),
        get(:net_liquidation_value, "NetLiquidation"),
        get(:total_cash_value, "TotalCashValue")
    )
end

accounts(ib::IB) = map(n -> IBAccount(ib, n), managed_accounts(ib))

managed_accounts(ib::IB) = ib.managedAccounts()

function positions(ib::IB, account::IBAccount)
    ps = Dict{String,Position}()
    for (c, position, mkt_price, mkt_value, avg_cost, _, _, _account) in ib.portfolio()
        account.name == _account || continue

        Tname = pyname(c)
        currency = parse_currency(c.currency)
        c = Tname == "Stock" ? Stock(c.symbol, c.exchange, currency) :
            error("Unknown position type: $Tname")
        sym = fullsymbol(c)

        @assert !haskey(ps, sym)
        ps[sym] = Position(c, position, mkt_price * currency, mkt_value * currency, avg_cost * currency)
    end
    return ps
end

####
#### Conversion
####

function PyCall.PyObject(x::Stock)
    ib_insync.Stock(x.symbol, x.exchange, string(x.currency))
end

# function PyCall.PyObject(x::Stock)
#     ib_insync.Stock(x.symbol, x.exchange, string(x.currency))
# end
#
# Base.convert(::Type{<:Stock}, o::PyObject) = Stock(o.symbol, o.exchange, parse_currency(o.currency))


text:
# test_local_channels_pids.jl

@info("NOT FULLY TESTED YET!!!")

using Pkg
Pkg.activate(".")
Pkg.precompile()

# start docker containers
cmd = `docker start test_sshd`
run(cmd)

cmd = `docker start test_sshd2`
run(cmd)

cmd = `docker ps`
run(cmd)

sleep(5)

# enable distrbuted computing
using Distributed

#addprocs(4; exeflags=`--project=$(Base.active_project())`)
addprocs([("rob@172.17.0.2", 1), ("rob@172.17.0.3", 1)]; exeflags=`--project=$(Base.active_project())`, tunnel=true, dir="/home/rob")

#addprocs([("pi@192.168.2.3", 1)]; exename=`/home/pi/julia/julia-1.3.1/bin/julia`, dir="/home/pi")

# remove processes > 3
while length(procs()) ≥ 4
    rmprocs(procs()[length(procs())])
end

gl_pid = procs()[2] # general ledger
ar_pid = procs()[3] # accounts receivable (orders/bankstatements)

# activate the packages
@everywhere begin
    using AppliSales
    using AppliGeneralLedger
    using AppliAR
    using Query
end;

# get tasks
include("./api/api3.jl")

# start dispatcher
rx = dispatcher(gl_pid, ar_pid)

# start application
put!(rx, "START")

# processing the uppaid invoices
stms = AppliAR.read_bank_statements(PATH_CSV)
sleep(15)
put!(rx, stms)
sleep(15)

# unkown type
test = "Test unkown type"
put!(rx, test)

# print aging report
using DataFrames
r1 = @fetchfrom ar_pid report(path_unpaid=PATH_DB, path_paid=PATH_DB_PAID)
result = DataFrame(r1)
println("\nUnpaid invoices\n===============")
show(result)

# print general ledger
r2 = @fetchfrom gl_pid AppliGeneralLedger.read_from_file(PATH_LEDGER)
df = DataFrame(r2)
println("\nGeneral Ledger mutations\n========================")
show(df)

df2 = r2 |> @filter(_.accountid == 1300) |> DataFrame
println("\nAccounts receivable\n===================")
show(df2)

println("\nBalance Accounts Receivable: $(sum(df2.debit) - sum(df2.credit))")

# open shell in container
cmd = `ssh rob@172.17.0.2`
@info("after run(cmd) is activated: goto console, press Enter, and rm test3_* files. Leave the container with Ctrl-D")
run(cmd)

# open shell in container
cmd = `ssh rob@172.17.0.3`
@info("after run(cmd) is activated: goto console, press Enter, and rm test3_invoicing.sqlite*. Leave the container with Ctrl-D")
run(cmd)
@info("Ctrl-L to clean the consule. Close julia with Ctrl-D.")


text:
#=
 Algorithms for triple-word arithmetic
 Nicolas Fabiano, Jean-Michel Muller, Joris Picot
=#

function vecsum(x0::T,x1::T,x2::T) where {T}
   x0 === 0.0 && return (two_sum(x1, x2)..., 0.0)
   s1, e2 = two_sum(x1, x2)
   e0, e1 = two_sum(x0, s1)
   return (e0, e1, e2)
end

function vseb(a0::T, a1::T, a2::T, a3::T, a4::T, a5::T) where {T}
    a2 === 0.0 && return (a0, a1, a2)
	r0, e1 = two_sum(a0, a1)
	if e1 != 0.0
	  y0 = r0
	else
	  e1 = r0
	end
	y1, y2  = two_sum(e1, a2)
	return (y0, y1, y2)
end

function vseb(a0::T, a1::T, a2::T) where {T}
    a2 === 0.0 && return (a0, a1, a2)
	r0, e1 = two_sum(a0, a1)
	if e1 != 0.0
	  y0 = r0
	else
	  e1 = r0
	end
	y1, y2  = two_sum(e1, a2)
	return (y0, y1, y2)
end

function vecseb(x0::T, x1::T, x2::T) where {T}
   a0, a1, a2 = vecsum(x0, x1, x2)
   return vseb(a0, a1, a2)
end

function tripleword(a::T, b::T, c::T) where {T}
    d0, d1 = two_sum(a,b)
    # e0, e1, e2 = vecsum(d0, d1, c)
    s1, e2 = two_sum(d1, c)
    e0, e1 = two_sum(d0, s1)
    # s0, s1, s2 = vseb(e0, e1, e2)
    r0, e1 = two_sum(e0, e1)
	if e1 != 0.0
	  s0 = r0
	else
	  e1 = r0
	end
	s1, s2  = two_sum(e1, e2)
	return (s0, s1, s2)
end

using SetRounding

function Float64(x0::Float64, x1::Float64, x2::Float64)
    s, e = two_hilo_sum(x0, x1+x1)
    if e !== 0.0
        y = x0 + x1
    elseif x2 > 0.0
        setrounding(Float64, RoundUp)
        y = x0 + x1
        setrounding(Float64, RoundNearest)
    elseif x2 < 0.0
        setrounding(Float64, RoundDown)
        y = x0 + x1
        setrounding(Float64, RoundNearest)
    else
        y = x0 + x1
    end
    return y
end


text:
abstract TreeNode

type BranchNode <: TreeNode
    children :: Array{BranchNode, 1}
    data
    extrainfo
end

type NullNode <: TreeNode
end

const nullnode = NullNode()

function leaves_of(root::TreeNode)
    code = Int64[]
    function traverse(node::TreeNode)
        if node == nullnode
            return
        end
        if length(node.children) == 0
            produce((node.data, copy(code)))    # notice that we should copy the current state of code
        end
        for (index, child) in enumerate(node.children)
            push!(code, index)
            traverse(child)
            pop!(code)
        end
    end
    Task(() -> traverse(root))
end

function internal_nodes_of(root::TreeNode)
    function traverse(node::TreeNode)
        if node == nullnode
            return
        end
        if length(node.children) != 0
            produce(node)
        end
        for child in node.children
            traverse(child)
        end
    end
    Task(() -> traverse(root))
end

function average_height(tree::TreeNode)
    (h, c) = (0, 0)
    for (_, path) in leaves_of(tree)
        h += length(path)
        c += 1
    end
    h / c
end


text:
using CartesianGeneticProgramming
using PyCall
using Cambrian
using ArgParse
import Distances
import Random
import Formatting
import Base.GC

```
Continuous gym problems using pybullet
```

"Gym BERL environment"
mutable struct BulletEnv <: BERLenv
    name::String
    memory::Array # Array of storable data
    gym
    pybullet
    env_name
    gen
    n_steps::Int64
end

environments["pybullet"]=BulletEnv

function BulletEnv(cfg::Dict)
    pybullet_envs = pyimport("pybullet_envs")
    gym = pyimport("gym")
    cfg["env"] = get!(cfg, "gym_env", "AntBulletEnv-v0") #'HopperBulletEnv-v0'
    env = gym.make(cfg["env"])
    cfg["n_out"] = length(env.action_space.sample())
    cfg["n_in"] = length(env.observation_space.sample())
    gen = get!(cfg, "seed", 0)
    Random.seed!(gen)
    cfg["nsteps"] = 0

    cfg["env"]=

    BulletEnv(cfg["env"], [], gym, pybullet_envs, cfg["env"], gen, 0)
end

"Returns the fitness of one Cambrian Individual in this environment"
function fitness(indiv::Cambrian.Individual, env::BulletEnv)
    game = env.gym.make(env.env_name)
    game.seed(env.gen)
    obs = game.reset()
    total_reward = 0.0
    done = false

    max_obs = 2*pi

    while ~done
        action = process(indiv, obs ./ max_obs)
        obs, reward, done, _ = game.step(action)
        newmax = maximum(abs.(obs))
        if newmax > max_obs
            println("Increased max_obs from ", max_obs, " to ", newmax)
            max_obs = newmax
        end
        total_reward += reward
        env.n_steps += 1
    end

    game.close()
    game = nothing
    Base.GC.gc()
    [total_reward]
end

"Sets up the env for the next generation"
function new_gen!(env::BulletEnv, cfg::Dict=nothing)
    #env.gen += 1
end