text:
using RungeKutta: name, order, eachstage, nstages, coefficients, weights, nodes, to_array, to_file, from_file

@testset "$(rpad("Tableau",80))" begin

    for s in 1:5
        for T ∈ (Float64, BigFloat)

            a = rand(T,s,s)
            b = rand(T,s)
            c = rand(T,s)

            tab1 = Tableau(:Test, 2s, a, b, c)
            tab2 = Tableau(:Test, 2s, s, a, b, c)
            tab3 = Tableau{T}(:Test, 2s, a, b, c)
            tab4 = Tableau{T}(:Test, 2s, s, a, b, c)
            tabϵ = Tableau(:Test, 2s, a .+ 4eps(T) .* rand(T,s,s),
                                      b .+ 4eps(T) .* rand(T,s),
                                      c .+ 4eps(T) .* rand(T,s))

            @test tab1 ≈ tab2 ≈ tab3 ≈ tab4
            @test tab1 == tab2 == tab3 == tab4
            @test hash(tab1) == hash(tab2) == hash(tab3) == hash(tab4)

            @test tab1 ≈ tabϵ
            @test tab1 != tabϵ
            @test hash(tab1) != hash(tabϵ)

            @test isequal(tab1, tab2)
            @test isequal(tab1, tab3)
            @test isequal(tab1, tab4)
            @test !isequal(tab1, tabϵ)

            @test eltype(tab1) == eltype(tab2) == eltype(tab3) == eltype(tab4) == T

            @test name(tab1) == :Test
            @test order(tab1) == 2s
            @test eachstage(tab1) == 1:s
            @test nstages(tab1) == s
            @test coefficients(tab1) == a
            @test weights(tab1) == b
            @test nodes(tab1) == c

            @test tab1.a == tab2.a == tab3.a == tab4.a == a
            @test tab1.b == tab2.b == tab3.b == tab4.b == b
            @test tab1.c == tab2.c == tab3.c == tab4.c == c

            @test tab1.â == tab2.â == tab3.â == tab4.â == zero(a)
            @test tab1.b̂ == tab2.b̂ == tab3.b̂ == tab4.b̂ == zero(b)
            @test tab1.ĉ == tab2.ĉ == tab3.ĉ == tab4.ĉ == zero(c)

            @test tab1 == Tableau(tab1.name, tab1.o, to_array(tab1))
            @test tab1 == Tableau(tab1.name, tab1.o, convert(Matrix, tab1))
            @test tab1 == Tableau(tab1.name, tab1.o, convert(Array, tab1))
            @test tab1 == convert(Tableau, convert(Matrix, tab1); name=tab1.name, o=tab1.o)
            @test tab1 == convert(Tableau, convert(Array,  tab1); name=tab1.name, o=tab1.o)
            @test convert(Matrix, tab1) == convert(Matrix{T}, tab1)
            @test convert(Array,  tab1) == convert(Array{T},  tab1)

            @test startswith(repr(tab1), "\nRunge-Kutta Tableau")
            @test startswith(repr(MIME("text/markdown"), tab1), "Runge-Kutta Tableau")

        end

        T1 = BigFloat
        T2 = Float64

        a = rand(T1,s,s)
        b = rand(T1,s)
        c = rand(T1,s)

        tab1 = Tableau{T1}(:Test, 2s, s, a, b, c)
        tab2 = Tableau{T2}(:Test, 2s, s, a, b, c)

        @test tab2.â == T2.(tab1.a .- tab2.a)
        @test tab2.b̂ == T2.(tab1.b .- tab2.b)
        @test tab2.ĉ == T2.(tab1.c .- tab2.c)

        tmp = mktempdir()
        to_file(tmp, tab2)
        tabf = from_file(tmp, "Test")
        rm(tmp, recursive=true)

        @test tabf.o == tab2.o
        @test tabf.s == tab2.s
        @test tabf.a == tab2.a
        @test tabf.b == tab2.b
        @test tabf.c == tab2.c

    end

end


text:
module TestOnlineStats

include("preamble.jl")
using OnlineStats: Mean

@testset "reducingfunction" begin
    @test reducingfunction(Mean()) isa Function
    @test reducingfunction(Map(identity), Mean()) isa Function
end

end

ext: jl


text:
function genlib(cpucompiler::String, gpucompiler::String, coredir::String)
# Examples: genlib("g++","");
#           genlib("g++","nvcc");

mydir = pwd();
cd(coredir);

if length(cpucompiler)>0
    str = `$cpucompiler -fPIC -O3 -c commonCore.cpp`;
    run(str);
    str = `ar rvs commonCore.a commonCore.o`;
    run(str);

    str = `$cpucompiler -fPIC -O3 -c opuCore.cpp`;
    run(str);
    str = `ar rvs opuCore.a opuCore.o`;
    run(str);

    # str = `$cpucompiler -fPIC -O3 -c cpuCore.cpp -fopenmp`;
    # run(str);
    # str = `ar rvs cpuCore.a cpuCore.o`;
    # run(str);
end

if length(gpucompiler)>0
    str = `$gpucompiler -D_FORCE_INLINES -O3 -c --compiler-options '-fPIC' gpuCore.cu`;
    run(str);
    str = `ar rvs gpuCore.a gpuCore.o `;
    run(str);
end

if Sys.isapple()
#    mv("*.a",pwd() * "/Mac");
    if !isdir("Mac")
        mkdir("Mac")
    end
    if isfile("commonCore.a")
        run(`mv commonCore.a Mac`);
    end
    if isfile("opuCore.a")
        run(`mv opuCore.a Mac`);
    end
    if isfile("cpuCore.a")
        run(`mv cpuCore.a Mac`);
    end
    if isfile("gpuCore.a")
        run(`mv gpuCore.a Mac`);
    end
elseif Sys.isunix()
    if !isdir("Linux")
        mkdir("Linux")
    end
    if isfile("commonCore.a")
        run(`mv commonCore.a Linux`);
    end
    if isfile("opuCore.a")
        run(`mv opuCore.a Linux`);
    end
    if isfile("cpuCore.a")
        run(`mv cpuCore.a Linux`);
    end
    if isfile("gpuCore.a")
        run(`mv gpuCore.a Linux`);
    end
elseif Sys.iswindows()
    if !isdir("Windows")
        mkdir("Windows")
    end
    if isfile("commonCore.a")
        run(`mv commonCore.a Windows`);
    end
    if isfile("opuCore.a")
        run(`mv opuCore.a Windows`);
    end
    if isfile("cpuCore.a")
        run(`mv cpuCore.a Windows`);
    end
    if isfile("gpuCore.a")
        run(`mv gpuCore.a Windows`);
    end
end

rm("commonCore.a", force=true);
rm("opuCore.a", force=true);
rm("cpuCore.a", force=true);
rm("gpuCore.a", force=true);
rm("commonCore.o", force=true);
rm("opuCore.o", force=true);
rm("cpuCore.o", force=true);
rm("gpuCore.o", force=true);

cd(mydir);

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# This file investigates MSE with respect to free variables

using FractalTools
using Makie

# Construct interpolation data
f(x, y) = x^2 + y^2 + 1
precisions = 2 .^ (7 : 11)

npts = 100
mse = map(precisions) do prec
    # Set precision
    setprecision(prec) do
        freevar = BigFloat(0.001)

        # Construct interpolation data
        ngon = [
            BigFloat.([0.0, 0.0]),
            BigFloat.([1.0, 0.0]),
            BigFloat.([0.5, 1.0])
        ]
        pts = getdata(f, ngon, npts)

        # Construct interpolant
        interp = interpolate(pts, Interp2D(freevar))

        # Construct test data
        tpts = getdata(ngon, npts)
        ntpts = length(tpts)
        @info tpts[1][1].prec
        @info interp.ifs.ws[1].A[1, 1].prec

        # Compute MSE
        fvals = map(pt -> f(pt...), tpts)
        ivals = map(pt -> interp(pt...), tpts)
        sum((fvals - ivals).^2) / ntpts
    end
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Free Variable", ylabel="MSE", title="2D Interpolation MSE")
stem!(ax, precisions, mse, color=:black)
save(joinpath(@__DIR__, "interp2d_mse.png"), fig)
display(fig)


text:
# Global constant allowing to chose between ForwardDiff and Zygote
# for hyperparameter optimization
const ADBACKEND = Ref(:Zygote)

function setadbackend(ad_backend::Symbol)
    (ad_backend == :ForwardDiff || ad_backend == :Zygote) ||
        error("Wrong AD Backend symbol, options are :ForwardDiff or :Zygote")
    return ADBACKEND[] = ad_backend
end

# Type piracy but will be introduced in a future version
# See https://github.com/FluxML/Optimisers.jl/pull/24
Optimisers.init(::Any, ::Any) = nothing

## Generic fallback when gradient is nothing
update_kernel!(::Any, ::Any, ::Nothing, ::Any) = nothing
update_Z!(::Any, ::Any, ::Nothing, ::Any) = nothing

## Generic fallback when optimizer is nothing
update_kernel!(::Nothing, ::Kernel, ::NamedTuple, ::NamedTuple) = nothing
update_kernel!(::Nothing, ::AbstractVector, ::NamedTuple, ::NamedTuple) = nothing
update_Z!(::Nothing, ::AbstractVector, ::AbstractVector) = nothing
update_Z!(::Nothing, ::AbstractVector, ::NamedTuple) = nothing

## Updating prior mean parameters ##
# function update!(μ::PriorMean, g::AbstractVector, X::AbstractVector, state)
# return update!(μ, g, X, state)
# end

## Updating kernel parameters ##

## ForwardDiff.jl approach (with destructure())
function update!(opt, k::Kernel, Δ::AbstractVector)
    ps = params(k)
    i = 1
    for p in ps
        d = length(p)
        Δp = Δ[i:(i + d - 1)]
        Δlogp = Optimise.apply!(opt, p, p .* Δp)
        @. p = exp(log(p) + Δlogp)
        i += d
    end
end

## Zygote.jl approach with named tuple
function update_kernel!(opt, k::Union{Kernel,Transform}, g::NamedTuple, state::NamedTuple)
    return NamedTuple(
        map(keys(g)) do fieldname
            Pair(
                fieldname,
                update_kernel!(
                    opt,
                    getfield(k, fieldname),
                    getfield(g, fieldname),
                    getfield(state, fieldname),
                ),
            )
        end,
    )
end

function update_kernel!(opt, x::AbstractArray, g::AbstractArray, state)
    state, Δ = Optimisers.apply(opt, state, x, x .* g)
    @. x = exp(log(x) + Δ) # Always assume that parameters need to be positive
    return state
end

## Updating inducing points
function update_Z!(opt, Z::AbstractVector, Z_grads::AbstractVector, state)
    return map(Z, Z_grads, state) do z, zgrad, st
        st, ΔZ = Optimisers.apply(opt, st, z, zgrad)
        z .+= ΔZ
        return st
    end
end

function update_Z!(opt, Z::Union{ColVecs,RowVecs}, Z_grads::NamedTuple, state)
    st, Δ = Optimisers.apply(opt, state, Z.X, Z_grads.X)
    Z.X .+= Δ
    return st
end


text:
using Test
using TopOpt
import Ferrite
import GeometryBasics
import Makie

# Define problem, can also be imported from .inp files
# nels = (30,10)
# sizes = (1.0,1.0)
nels = (20, 10, 10)
sizes = (1.0,1.0,1.0)
E = 1.0;
ν = 0.3;
force = -1.0;
# problem = PointLoadCantilever(nels, sizes, E, ν, force)
problem = PointLoadCantilever(Val{:Linear}, nels, sizes, E, ν, force);
# @show TopOpt.TopOptProblems.getdim(problem)

# Build element stiffness matrices and force vectors
# einfo = ElementFEAInfo(problem);

# Assemble global stiffness matrix and force vector
# ginfo = assemble(problem, einfo);

# Solve for node displacements
# u = ginfo.K \ ginfo.f
# @show u

solver = FEASolver(Displacement, Direct, problem)
solver()
u = solver.u

# ncells = Ferrite.getncells(problem)
# result_mesh = GeometryBasics.Mesh(problem, ones(ncells));
# GLMakie.mesh(result_mesh);

if get(ENV, "CI", nothing) != "true"
    @eval begin
        using TopOpt.TopOptProblems.Visualization: visualize
        # fig = visualize(problem)
        fig = visualize(problem, u; default_exagg_scale=0.07, scale_range=10.0, vector_linewidth=3, vector_arrowsize=0.5)
        # Makie.display(fig)
    end
end

# using GeometryBasics, Makie
# using Ferrite
# full_topology = ones(Float64, Ferrite.getncells(problem))
# result_mesh = GeometryBasics.Mesh(problem, full_topology);
# mesh(result_mesh);

# # Get the stiffness matrices for buckling analysis
# K, Kσ = TopOpt.TopOptProblems.buckling(problem, ginfo, einfo)

# using IterativeSolvers

# # Find the maximum eigenvalue of system Kσ x = 1/λ K x
# r = lobpcg(Kσ.data, K.data, true, 1)

# # Minimum eigenvalue of the system K x = λ Kσ x
# λ = 1/r.λ[1]

# https://mohamed82008.github.io/ScienceLounge/text/2018/07/11/The-buckle-dance/


text:
module TestBasicGraphs
using Test

import Graphs as SimpleGraphs
import MetaGraphs

using Catlab.Graphs.BasicGraphs

# Graphs
########

g = Graph()
add_vertex!(g)
add_vertices!(g, 2)
@test nv(g) == 3
@test ne(g) == 0

add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
@test ne(g) == 2
@test ne(g, 1, 2) == 1
@test has_edge(g, 1, 2)
@test !has_edge(g, 1, 3)
@test outneighbors(g, 2) == [3]
@test inneighbors(g, 2) == [1]
@test degree(g, 2) == 2
@test collect(all_neighbors(g, 2)) == [1,3]

add_edge!(g, 1, 2)
@test ne(g) == 3
@test ne(g, 1, 2) == 2
@test collect(edges(g, 1, 2)) == [1,3]
@test outneighbors(g, 1) == [2,2]
@test inneighbors(g, 1) == []
@test degree(g, 1) == 2
@test SimpleGraphs.DiGraph(g) == SimpleGraphs.path_digraph(3)

g = Graph(4)
add_edges!(g, [1,2,3], [2,3,4])
@test SimpleGraphs.DiGraph(g) == SimpleGraphs.path_digraph(4)
@test Graph(SimpleGraphs.path_digraph(4)) == g

rem_edge!(g, 3, 4)
@test ne(g) == 2
@test src(g) == [1,2]
@test tgt(g) == [2,3]
rem_vertex!(g, 2)
@test nv(g) == 3
@test ne(g) == 0

# Error handling.
g = Graph(2)
add_edge!(g, 1, 2)
@test_throws Exception add_edge!(g, 2, 3) # nonexistent target vertex
@test ne(g) == 1
@test isempty(outneighbors(g, 2))

# Induced subgraphs.
g = Graph(3)
add_edges!(g, [1,1,2], [3,3,3])
sub = induced_subgraph(g, [1,3])
@test nv(sub) == 2
@test src(sub) == [1,1]
@test tgt(sub) == [2,2]

# Symmetric graphs
##################

g = SymmetricGraph()

add_vertices!(g, 3)
@test nv(g) == 3
@test ne(g) == 0

add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
@test ne(g) == 4
@test collect(edges(g, 1, 2)) == [1]
@test neighbors(g, 1) == [2]
@test neighbors(g, 2) == [1,3]
@test neighbors(g, 3) == [2]
@test SimpleGraphs.Graph(g) == SimpleGraphs.path_graph(3)
@test SymmetricGraph(SimpleGraphs.path_graph(3)) == g

g = SymmetricGraph(4)
add_edges!(g, [1,2,3], [2,3,4])
lg = SimpleGraphs.DiGraph(map(SimpleGraphs.Edge, [1,2,3,2,3,4], [2,3,4,1,2,3]))
@test SimpleGraphs.DiGraph(g) == lg

rem_edge!(g, 3, 4)
@test ne(g) == 4
@test neighbors(g, 3) == [2]
@test neighbors(g, 4) == []
rem_vertex!(g, 2)
@test nv(g) == 3
@test ne(g) == 0

# Reflexive graphs
##################

g = ReflexiveGraph()
add_vertex!(g)
add_vertices!(g, 2)
@test nv(g) == 3
@test ne(g) == 3
@test refl(g,1) == 1
@test refl(g) == [1,2,3]
@test src(g) == [1,2,3]
@test tgt(g) == [1,2,3]

add_edges!(g, [1,2], [2,3])
add_edge!(g, 1, 3)
@test ne(g) == 6
@test src(g, 4:6) == [1,2,1]
@test tgt(g, 4:6) == [2,3,3]

g = ReflexiveGraph(4)
add_edges!(g, [1,2,3], [2,3,4])
rem_edge!(g, 3, 4)
@test ne(g) == 6
@test src(g, 5:6) == [1,2]
@test tgt(g, 5:6) == [2,3]
rem_vertex!(g, 2)
@test nv(g) == 3
@test ne(g) == 3
@test refl(g) == [1,2,3]
@test src(g) == [1,2,3]
@test tgt(g) == [1,2,3]

# Symmetric reflexive graphs
############################

g = SymmetricReflexiveGraph()
add_vertex!(g)
add_vertices!(g, 2)
@test nv(g) == 3
@test refl(g) == [1,2,3]
@test inv(g) == [1,2,3]
@test src(g) == [1,2,3]
@test tgt(g) == [1,2,3]

add_edge!(g, 1, 3)
@test ne(g) == 5
@test src(g, 4:5) == [1,3]
@test tgt(g, 4:5) == [3,1]
@test inv(g, 4:5) == [5,4]

g = SymmetricReflexiveGraph(4)
add_edges!(g, [1,2,3], [2,3,4])
rem_edge!(g, 3, 4)
@test ne(g) == 8
rem_vertex!(g, 2)
@test nv(g) == 3
@test ne(g) == 3

# Half-edge graphs
##################

g = HalfEdgeGraph()

add_vertices!(g, 2)
@test nv(g) == 2
@test isempty(half_edges(g))

add_edge!(g, 1, 2)
add_edge!(g, 2, 1)
@test collect(half_edges(g)) == [1,2,3,4]
@test vertex(g) == [1,2,2,1]
@test inv(g) == [2,1,4,3]
@test half_edges(g, 1) == [1,4]
@test half_edges(g, 2) == [2,3]

add_dangling_edge!(g, 1)
add_dangling_edges!(g, [2,2])
@test length(half_edges(g)) == 7
@test vertex(g, 5:7) == [1,2,2]
@test inv(g, 5:7) == [5,6,7]

g = HalfEdgeGraph(4)
add_edges!(g, [1,2,3], [2,3,4])
lg = SimpleGraphs.Graph(map(SimpleGraphs.Edge, [1,2,3], [2,3,4]))
@test SimpleGraphs.Graph(g) == lg

rem_edge!(g, 3, 4)
@test Set(zip(vertex(g), vertex(g,inv(g)))) == Set([(1,2),(2,1),(2,3),(3,2)])
add_edge!(g, 3, 4)
rem_edge!(g, last(half_edges(g)))
@test Set(zip(vertex(g), vertex(g,inv(g)))) == Set([(1,2),(2,1),(2,3),(3,2)])
rem_vertex!(g, 2)
@test nv(g) == 3
@test isempty(half_edges(g))

# Weighted graphs
#################

g = WeightedGraph{Float64}(4)
add_edges!(g, 1:3, 2:4, weight=[0.25, 0.5, 0.75])
@test weight(g, 1) == 0.25
@test weight(g) == [0.25, 0.5, 0.75]

mg = MetaGraphs.MetaDiGraph{Int,Float64}(4)
MetaGraphs.add_edge!(mg, 1, 2, :weight, 0.25)
MetaGraphs.add_edge!(mg, 2, 3, :weight, 0.5)
MetaGraphs.add_edge!(mg, 3, 4, :weight, 0.75)
@test MetaGraphs.MetaDiGraph(g) == mg

end

max_forks_repo_name: epatters/CompCat


text:
# Different ways to simulate molecules

export
    accelerations,
    VelocityVerlet,
    simulate!,
    VelocityFreeVerlet

"""
    accelerations(simulation, neighbours; parallel=true)

Calculate the accelerations of all atoms using the general and specific
interactions and Newton's second law.
"""
function accelerations(s::Simulation, neighbours; parallel::Bool=true)
    n_atoms = length(s.coords)

    if parallel && nthreads() > 1 && n_atoms >= 100
        forces_threads = [zero(s.coords) for i in 1:nthreads()]

        # Loop over interactions and calculate the acceleration due to each
        for inter in values(s.general_inters)
            if inter.nl_only
                @threads for ni in 1:length(neighbours)
                    i, j = neighbours[ni]
                    force!(forces_threads[threadid()], inter, s, i, j)
                end
            else
                @threads for i in 1:n_atoms
                    for j in 1:i
                        force!(forces_threads[threadid()], inter, s, i, j)
                    end
                end
            end
        end

        forces = sum(forces_threads)
    else
        forces = zero(s.coords)

        for inter in values(s.general_inters)
            if inter.nl_only
                for ni in 1:length(neighbours)
                    i, j = neighbours[ni]
                    force!(forces, inter, s, i, j)
                end
            else
                for i in 1:n_atoms
                    for j in 1:i
                        force!(forces, inter, s, i, j)
                    end
                end
            end
        end
    end

    for inter_list in values(s.specific_inter_lists)
        for inter in inter_list
            force!(forces, inter, s)
        end
    end

    for i in 1:n_atoms
        forces[i] /= s.atoms[i].mass
    end

    return forces
end

"""
    VelocityVerlet()

The velocity Verlet integrator.
"""
struct VelocityVerlet <: Simulator end

"""
    simulate!(simulation; parallel=true)
    simulate!(simulation, n_steps; parallel=true)
    simulate!(simulation, simulator, n_steps; parallel=true)

Run a simulation according to the rules of the given simulator.
Custom simulators should implement this function.
"""
function simulate!(s::Simulation,
                    ::VelocityVerlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    # See https://www.saylor.org/site/wp-content/uploads/2011/06/MA221-6.1.pdf for
    #   integration algorithm - used shorter second version
    n_atoms = length(s.coords)
    neighbours = find_neighbours(s, nothing, s.neighbour_finder, 0,
                                    parallel=parallel)
    accels_t = accelerations(s, neighbours, parallel=parallel)
    accels_t_dt = zero(s.coords)

    @showprogress for step_n in 1:n_steps
        for logger in values(s.loggers)
            log_property!(logger, s, step_n)
        end

        # Update coordinates
        for i in 1:length(s.coords)
            s.coords[i] += s.velocities[i] * s.timestep + accels_t[i] * (s.timestep ^ 2) / 2
            s.coords[i] = adjust_bounds.(s.coords[i], s.box_size)
        end

        accels_t_dt = accelerations(s, neighbours, parallel=parallel)

        # Update velocities
        for i in 1:length(s.velocities)
            s.velocities[i] += (accels_t[i] + accels_t_dt[i]) * s.timestep / 2
        end

        apply_thermostat!(s, s.thermostat)
        neighbours = find_neighbours(s, neighbours, s.neighbour_finder, step_n,
                                        parallel=parallel)

        accels_t = accels_t_dt
        s.n_steps_made[1] += 1
    end
    return s
end

"""
    VelocityFreeVerlet()

The velocity-free Verlet integrator, also known as the Störmer method.
In this case the `velocities` given to the `Simulator` act as the previous step
coordinates for the first step.
"""
struct VelocityFreeVerlet <: Simulator end

function simulate!(s::Simulation,
                    ::VelocityFreeVerlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    n_atoms = length(s.coords)
    neighbours = find_neighbours(s, nothing, s.neighbour_finder, 0,
                                    parallel=parallel)
    coords_last = s.velocities

    @showprogress for step_n in 1:n_steps
        for logger in values(s.loggers)
            log_property!(logger, s, step_n)
        end

        accels_t = accelerations(s, neighbours, parallel=parallel)

        # Update coordinates
        coords_copy = s.coords
        for i in 1:length(s.coords)
            s.coords[i] = 2 * s.coords[i] - coords_last[i] + accels_t[i] * s.timestep ^ 2
            s.coords[i] = adjust_bounds.(s.coords[i], s.box_size)
        end
        coords_last = coords_copy

        apply_thermostat!(s, s.thermostat)
        neighbours = find_neighbours(s, neighbours, s.neighbour_finder, step_n,
                                        parallel=parallel)

        s.n_steps_made[1] += 1
    end
    return s
end

function simulate!(s::Simulation, n_steps::Integer; parallel::Bool=true)
    simulate!(s, s.simulator, n_steps, parallel=parallel)
end

function simulate!(s::Simulation; parallel::Bool=true)
    simulate!(s, s.n_steps - first(s.n_steps_made), parallel=parallel)
end


text:
include("transit_poly.jl")

# Randomizer seed
using Random
Random.seed!(42)

nc = 2+ceil(Int64,rand()*20)
r=rand(); b=rand()*(1+r); c = rand(nc)
#r=1.2*rand(); b=0.0; c = rand(nc)
#r = rand(); b = rand(); c = [1.0,1.0,10.0]
#r=0.7821623254303613; b= 0.5; c = [0.0,1.0]

flux = transit_poly_c(r,b,c)
@time flux = transit_poly_c(r,b,c)

# Now, integrate by hand:
function transit_poly_int_c(r,b,c)
s_1 = maximum([0.0,b-r]); s_2=minimum([1.0,b+r])
ns = 500; ds=(s_2-s_1)/ns
s = linspace(s_1+.5*ds,s_2-.5*ds,ns)
fobs = zero(r)
for j=1:ns
  if s[j] < r-b
    dphi = 2pi
  else
    dphi = 2*acos((s[j]^2+b^2-r^2)/(2*b*s[j]))
  end
  z = sqrt(1-s[j]^2)
  imu = c[1] + c[2]*z
  for n=2:length(c)-1
    imu += c[n+1]*(z^n*(n+2)-n*z^(n-2))
  end
  fobs += s[j]*ds*dphi*imu
end
fobs /= pi*(c[1]+2/3*c[2])
fobs = 1-fobs
#println("r: ",r," b: ",b,"f_an: ",flux," flux: ",flux0," f_num: ",fobs," f_exp: ",pi*(c[1]+2/3*c[2]))
return fobs
end
f_num = transit_poly_int_c(r,b,c)
@time f_num = transit_poly_int_c(r,b,c)
println("r: ",r," b: ",b," f_an: ",flux," f_num: ",f_num)

#nphi = 500; ns = 500; dphi = 2pi/nphi; ds=1/ns
#phi = linspace(.5*dphi,2pi-.5*dphi,nphi); s = linspace(.5*ds,1-.5*ds,ns)
#flux0 = zero(r); fobs = zero(r)
#for i=1:nphi, j=1:ns
#  x=s[j]*cos(phi[i]); y = s[j]*sin(phi[i])
#  z = sqrt(1-s[j]^2)
#  imu = c[1]+c[2]*z
#  for k=2:length(c)-1
#    imu += c[k+1]*(z^k*(k+2)-k*z^(k-2))
#  end
#  flux0 += s[j]*ds*dphi*imu
#  if (y-b)^2+x^2 > r^2
#    fobs += s[j]*ds*dphi*imu
#  end
#end
##println("r: ",r," b: ",b,"f_an: ",flux," flux: ",flux0," f_num: ",fobs/flux0," f_exp: ",pi*(c[1]+2/3*c[2]))
#println("r: ",r," b: ",b," f_an: ",flux," f_num: ",fobs/flux0)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using Base: max_ccall_threads
@testset "695.max-area-of-island.jl" begin
    @test max_area_of_island(hcat([[0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
                                   [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
                                   [0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                                   [0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0],
                                   [0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0],
                                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
                                   [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
                                   [0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0]]...)) == 6
    @test max_area_of_island(reshape([0,0,0,0,0,0,0,0], :, 1)) == 0
end


text:
cd(@__DIR__); include("setups/path128.jl")
gr(dpi = 200)

l = 64
𝛍 = nat_spec_filter(l, distDCT; σ = 0.05 * maximum(distDCT))
plot(0:N-1, 𝛍, legend = true, c = :blue,
     grid = false, frame = :box, marker = :circle, ms = 3, lw = 1, mswidth = 0,
     lab = latexstring("\\sigma=0.05 \\cdot d_{\\max}"))

𝛍 = nat_spec_filter(l, distDCT; σ = 0.1 * maximum(distDCT))
plot!(0:N-1, 𝛍, legend = true, c = :red,
     grid = false, frame = :box, marker = :circle, ms = 3, lw = 1, mswidth = 0,
     lab = latexstring("\\sigma=0.1 \\cdot d_{\\max}"))

𝛍 = nat_spec_filter(l, distDCT; σ = 0.2 * maximum(distDCT))
plot!(0:N-1, 𝛍, ylim = [0, 0.1], legend = true, c = :green,
     grid = false, frame = :box, marker = :circle, ms = 3, lw = 1, mswidth = 0,
     lab = latexstring("\\sigma=0.2 \\cdot d_{\\max}"), tickfontsize=10, legendfontsize = 12)
xticks!([0; 15:16:127], vcat(string("DC"), [string(k) for k in 15:16:127]))
plot!(xlab = latexstring("i"), ylab = latexstring("\\mu_{63}"), guidefontsize = 16)
plt = current()
savefig(plt, "../figs/Path128_NSGfilters_mu$(l-1).png")
display(plt)


text:
using QMTK
using BenchmarkTools

@generated function _kronprod(list...)
    ex = Expr(:call, :kron, :(list[1]), :(list[2]))
    for i in 3:length(list)
        ex = Expr(:call, kron, ex, :(list[$i]))
    end
    return ex
end

const sigmax = sprand(10, 10, 0.5)
⊗ = kron

println("--------------------")
println("baseline")
display(@benchmark kron(kron(kron(sigmax, sigmax), sigmax), sigmax))
println("--------------------")
println("QMTK")
display(@benchmark kronprod(sigmax, sigmax, sigmax, sigmax))
println("--------------------")
println("generated")
display(@benchmark _kronprod(sigmax, sigmax, sigmax, sigmax))
println("--------------------")
println("product")
display(@benchmark sigmax ⊗ sigmax ⊗ sigmax ⊗ sigmax)


text:
using Revise, Test
# This test should only be run if you have a very small inotify limit

@testset "inotify" begin
    logs, _ = Test.collect_test_logs() do
        Revise.track("common.jl")
    end
    sleep(0.1)
    @test !isempty(logs)
    @test any(rec->occursin("inotify", rec.message), logs)
end


text:
@inline UNITLESS_ABS2(x::Number) = abs2(x)
@inline UNITLESS_ABS2(x::AbstractArray) = sum(UNITLESS_ABS2, x)
@inline UNITLESS_ABS2(x::RecursiveArrayTools.ArrayPartition) = sum(UNITLESS_ABS2, x.x)
Base.mapreduce_empty(::typeof(UNITLESS_ABS2), op, T) = abs2(Base.reduce_empty(op, T))

@inline recursive_length(u::AbstractArray{<:Number}) = length(u)
@inline recursive_length(u::Number) = length(u)
@inline recursive_length(u::AbstractArray{<:AbstractArray}) = sum(recursive_length, u)
@inline recursive_length(u::RecursiveArrayTools.ArrayPartition) = sum(recursive_length, u.x)
@inline recursive_length(u::RecursiveArrayTools.VectorOfArray) = sum(recursive_length, u.u)

@inline ODE_DEFAULT_NORM(u::Union{AbstractFloat,Complex},t) = @fastmath abs(u)

@inline function ODE_DEFAULT_NORM(u::Array{T},t) where T<:Union{AbstractFloat,Complex}
    x = abs2(u[1])
    @inbounds for i in 2:length(u)
        x += abs2(u[i])
    end
    sqrt(real(x) / length(u))
end

@inline ODE_DEFAULT_NORM(u::StaticArrays.StaticArray{T},t) where T<:Union{AbstractFloat,Complex} =
                                            sqrt(real(sum(abs2,u)) / length(u))

@inline ODE_DEFAULT_NORM(u::AbstractArray,t) = sqrt(UNITLESS_ABS2(u)/recursive_length(u))
@inline ODE_DEFAULT_NORM(u,t) = norm(u)

@inline ODE_DEFAULT_ISOUTOFDOMAIN(u,p,t) = false
@inline function ODE_DEFAULT_PROG_MESSAGE(dt,u::Array,p,t)
    tmp = u[1]
    for i in eachindex(u)
        tmp = ifelse(u[i] > tmp,u[i],tmp)
    end
    "dt="*string(dt)*"\nt="*string(t)*"\nmax u="*string(tmp)
end
@inline ODE_DEFAULT_PROG_MESSAGE(dt,u,p,t) =
           "dt="*string(dt)*"\nt="*string(t)*"\nmax u="*string(maximum(abs.(u)))

@inline NAN_CHECK(x::Number) = isnan(x)
@inline NAN_CHECK(x::Float64) = isnan(x) || (x>1e50)
@inline NAN_CHECK(x::Enum) = false
@inline NAN_CHECK(x::AbstractArray) = any(NAN_CHECK, x)
@inline NAN_CHECK(x::RecursiveArrayTools.ArrayPartition) = any(NAN_CHECK, x.x)

@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u,p,t) = false
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::Union{Number,AbstractArray},p,t) = NAN_CHECK(u)


text:
using SequentialMonteCarlo
using RNGPool
import SMCExamples.LinearGaussian.defaultLGModel

include("test.jl")

setRNGs(0)

model, theta, ys, ko = defaultLGModel(10)
println(ko.logZhats)

numParticles = 1024*1024
numTrials = 2

## just run the algorithm a few times

testSMC(model, numParticles, numTrials, false)
testSMC(model, numParticles, numTrials, false, 0.5)
testSMCParallel(model, numParticles, numTrials, false)
testSMCParallel(model, numParticles, numTrials, false, 0.5)


text:
"""
    run_opt(ts_data::ClustData,opt_data::OptDataCEP,config::Dict{String,Any},optimizer::DataType)
Organizing the actual setup and run of the CEP-Problem. This function shouldn't be called by a user, but from within the other `run_opt`-functions
Required elements are:
- `ts_data`: The time-series data.
- `opt_data`: In this case the OptDataCEP that contains information on costs, nodes, techs and for transmission also on lines.
- `config`: This includes all the settings for the design optimization problem formulation.
- `optimizer`: The used optimizer, which could e.g. be Clp: `using Clp` `optimizer=Clp.Optimizer` or Gurobi: `using Gurobi` `optimizer=Gurobi.Optimizer`.
"""
function run_opt(ts_data::ClustData,
                opt_data::OptDataCEP,
                config::Dict{String,Any},
                optimizer::DataType)
  #Check the consistency of the data provided
  check_opt_data_cep(opt_data)
  #Setup the basic elements
  cep=setup_opt_basic(ts_data, opt_data, config, optimizer, config["optimizer_config"])
  #Setup the basic variables
  setup_opt_basic_variables!(cep, ts_data, opt_data)
  # Setup the demand
  if config["demand"]
    setup_opt_demand!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup the non dispachable electricity generation
  if config["non_dispatchable_generation"]
    setup_opt_non_dispatchable_generation!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup the dispatchable electricity generation
  if config["dispatchable_generation"]
    setup_opt_dispatchable_generation!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup conversion
  if config["conversion"]
    setup_opt_conversion!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup storage
  if config["storage"]
    setup_opt_storage!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup seasonal storage
  if config["seasonalstorage"]
    setup_opt_seasonalstorage!(cep, ts_data, opt_data, config["scale"])
  #Or intra-day storage (same level in first and last time step for each period)
  elseif config["storage"] && !(config["seasonalstorage"])
    setup_opt_simplestorage!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup transmission
  if config["transmission"]
      setup_opt_transmission!(cep, ts_data, opt_data, config["scale"])
  end
  #If lost load cost exist
  if !isempty(config["lost_load_cost"])
    setup_opt_lost_load!(cep, ts_data, opt_data, config["scale"])
  end
  #If lost emission exist
  if !isempty(config["lost_emission_cost"])
    setup_opt_lost_emission!(cep, ts_data, opt_data)
  end
  #If limit of emissions
  if !isempty(config["limit_emission"])
    setup_opt_limit_emission!(cep, ts_data, opt_data, config["scale"]; limit_emission=config["limit_emission"],  lost_emission_cost=config["lost_emission_cost"])
  end
  #If fixed_design_variables are provided, fix the installed capacities to them
  if haskey(config,"fixed_design_variables")
    setup_opt_fix_design_variables!(cep, ts_data, opt_data, config["scale"], config["fixed_design_variables"])
  end
  #Setup constraints that bind the capacities of different capacities with each other
  #If design variables are provided, this is already done by the provided variables
  if !haskey(config,"fixed_design_variables")
    setup_opt_intertech_cap!(cep, ts_data, opt_data, config["scale"])
  end
  # Add existing infrastructure to
  setup_opt_existing_infrastructure!(cep, ts_data, opt_data, config["scale"])
  # Limit the infrastructure expansion
  setup_opt_limit_infrastructure!(cep, ts_data, opt_data, config["scale"])
  if config["transmission"]
    # Setup the energy balences taking transmission into account
    setup_opt_energy_balance_transmission!(cep, ts_data, opt_data, config["scale"])
  else
    # Setup the energy balences with a copperplate assumption without any transmission restrictions between nodes
    setup_opt_energy_balance_copperplate!(cep, ts_data, opt_data, config["scale"])
  end
  #Setup the objective
  setup_opt_objective!(cep, ts_data, opt_data, config["scale"]; lost_load_cost=config["lost_load_cost"], lost_emission_cost=config["lost_emission_cost"])
  # solve and return the CEP
  return solve_opt_cep(cep, ts_data, opt_data, config)
end

"""
     run_opt(ts_data::ClustData,opt_data::OptDataCEP,config::Dict{String,Any},fixed_design_variables::Dict{String,Any},optimizer::DataTyple;lost_el_load_cost::Number=Inf,lost_CO2_emission_cost::Number)
This problem runs the operational optimization problem only, with fixed design variables.
provide the fixed design variables and the `config` of the previous step (design run or another opterational run)
Required elements are:
- `ts_data`: The time-series data, which should be be the original time-series data for this operational run. The `keys(ts_data.data)` need to match the `[time_series_name]-[node]`
- `opt_data`: In this case the OptDataCEP that contains information on costs, nodes, techs and for transmission also on lines. - Should be the same as in the design run.
- `config`: This includes all the previous settings for the design optimization problem formulation and ensures that the configuration is the same.
- `fixed_design_variables`: All the design variables that are determined by the previous design run.
- `optimizer`: The used optimizer, which could e.g. be Clp: `using Clp` `optimizer=Clp.Optimizer` or Gurobi: `using Gurobi` `optimizer=Gurobi.Optimizer`.
What you can change in the `config`:
- `lost_load_cost`: Dictionary with numbers indicating the lost load price per carrier (e.g. `electricity` in price/MWh should be greater than 1e6), give Inf for no SLACK and LL (Lost Load - a variable for unmet demand by the installed capacities)
- `lost_emission_cost`: Dictionary with numbers indicating the emission price/kg-emission (Suggestion: around 700), give Inf for no LE (Lost Emissions - a variable for emissions that will exceed the limit in order to provide the demand with the installed capacities)
"""
function run_opt(ts_data::ClustData,
                opt_data::OptDataCEP,
                config::Dict{String,Any},
                fixed_design_variables::Dict{String,Any},
                optimizer::DataType;
                lost_load_cost::Dict{String,Number}=Dict{String,Number}(),
                lost_emission_cost::Dict{String,Number}=Dict{String,Number}())
  # Add the fixed_design_variables and new setting for slack costs to the existing config
  set_config_cep!(config;fixed_design_variables=fixed_design_variables, lost_load_cost=lost_load_cost, lost_emission_cost=lost_emission_cost)

  return run_opt(ts_data,opt_data,config,optimizer)
end

"""
    run_opt(ts_data::ClustData,
           opt_data::OptDataCEP,
           optimizer::DataType;
           descriptor::String="",
           storage_type::String="none",
           demand::Bool=true,
           dispatchable_generation::Bool=true,
           non_dispatchable_generation::Bool=true,
           conversion::Bool=false,
           transmission::Bool=false,
           lost_load_cost::Dict{String,Number}=Dict{String,Number}(),
           lost_emission_cost::Dict{String,Number}=Dict{String,Number}(),
           limit_emission::Dict{String,Number}=Dict{String,Number}(),
           infrastructure::Dict{String,Array}=Dict{String,Array}("existing"=>["demand"],"limit"=>Array{String,1}()),
           scale::Dict{Symbol,Int}=Dict{Symbol,Int}(:COST => 1e9, :CAP => 1e3, :GEN => 1e3, :SLACK => 1e3, :INTRASTOR => 1e3, :INTERSTOR => 1e6, :FLOW => 1e3, :TRANS =>1e3, :LL => 1e6, :LE => 1e9),
           print_flag::Bool=true,
           optimizer_config::Dict{Symbol,Any}=Dict{Symbol,Any}(),
           round_sigdigits::Int=9)
Wrapper function for type of optimization problem for the CEP-Problem (NOTE: identifier is the type of `opt_data` - in this case OptDataCEP - so identification as CEP problem).
Required elements are:
- `ts_data`: The time-series data, which could either be the original input data or some aggregated time-series data. The `keys(ts_data.data)` need to match the `[time_series_name]-[node]`
- `opt_data`: The OptDataCEP that contains information on costs, nodes, techs and for transmission also on lines.
- `optimizer`: The used optimizer, which could e.g. be Clp: `using Clp` `optimizer=Clp.Optimizer` or Gurobi: `using Gurobi` `optimizer=Gurobi.Optimizer`.
Options to tweak the model are:
- `descriptor`: A name for the model
- `storage_type`: String `"none"` for no storage, `"simple"` to include simple (only intra-day storage), or `"seasonal"` to include seasonal storage (inter-day)
- `demand`: Bool `true` or `false` for technology-group
- `dispatchable_generation`: Bool `true` or `false` for technology-group
- `non_dispatchable_generation`: Bool `true` or `false` for technology-group
- `conversion`: Bool `true` or `false` for technology-group
- `transmission`:Bool `true` or `false` for technology-group. If no transmission should be modeled, a 'copperplate' is assumed with no transmission restrictions between the nodes
- `limit`: Dictionary with numbers limiting the kg.-emission-eq./MWh (e.g. `CO2` normally in a range from 5-1250 kg-CO2-eq/MWh), give Inf or no kw if unlimited
- `lost_load_cost`: Dictionary with numbers indicating the lost load price per carrier (e.g. `electricity` in price/MWh should be greater than 1e6), give Inf for no SLACK and LL (Lost Load - a variable for unmet demand by the installed capacities). Example: lost_load_cost=Dict{String,Number}("electricity"=>1e6)
- `lost_emission_cost`: Dictionary with numbers indicating the emission price/kg-emission (Suggestion: around 700), give Inf for no LE (Lost Emissions - a variable for emissions that will exceed the limit in order to provide the demand with the installed capacities). Example: lost_emission_cost=Dict{String,Number}("CO2"=>700)
- `infrastructure` : Dictionary with Arrays indicating which technology groups should have `existing` infrastructure (`"existing" => ["demand","dispatchable_generation"]`) and which technology groups should have infrastructure `limit`ed (`"limit" => ["non_dispatchable_generation"]`)
- `scale`: Dict{Symbol,Int} with a number for each variable (like `:COST`) to scale the variables and equations to similar quantities. Try to acchieve that the numerical model only has to solve numerical variables in a scale of 0.01 and 100. The following equation is used as a relationship between the real value, which is provided in the solution (real-VAR), and the numerical variable, which is used within the model formulation (VAR): real-VAR [`EUR`, `MW` or `MWh`] = scale[:VAR] ⋅ VAR.
- `descriptor`: String with the name of this paricular model like "kmeans-10-co2-500"
- `print_flag`: Bool to decide if a summary of the Optimization result should be printed.
- `optimizer_config`: Each Symbol and the corresponding value in the Dictionary is passed on to the `with_optimizer` function in addition to the `optimizer`. For Gurobi an example Dictionary could look like `Dict{Symbol,Any}(:Method => 2, :OutputFlag => 0, :Threads => 2)` more information can be found in the optimizer specific documentation.
- `round_sigdigits`: Can be used to round the values of the result to a certain number of `sigdigits`.
"""
function run_opt(ts_data::ClustData,
                 opt_data::OptDataCEP,
                 optimizer::DataType;
                 descriptor::String="",
                 storage_type::String="none",
                 demand::Bool=true,
                 dispatchable_generation::Bool=true,
                 non_dispatchable_generation::Bool=true,
                 conversion::Bool=false,
                 transmission::Bool=false,
                 lost_load_cost::Dict{String,Number}=Dict{String,Number}(),
                 lost_emission_cost::Dict{String,Number}=Dict{String,Number}(),
                 limit_emission::Dict{String,Number}=Dict{String,Number}(),
                 infrastructure::Dict{String,Array}=Dict{String,Array}("existing"=>["demand"],"limit"=>Array{String,1}()),
                 scale::Dict{Symbol,Int}=Dict{Symbol,Int}(:COST => 1e9, :CAP => 1e3, :GEN => 1e3, :SLACK => 1e3, :INTRASTOR => 1e3, :INTERSTOR => 1e6, :FLOW => 1e3, :TRANS =>1e3, :LL => 1e6, :LE => 1e9),
                 print_flag::Bool=true,
                 optimizer_config::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                 round_sigdigits::Int=9,
                 time_series_config::Dict{String,Any}=Dict{String,Any}())
   # Activated seasonal or simple storage corresponds with storage
   if storage_type=="seasonal"
       storage=true
       seasonalstorage=true
   elseif storage_type=="simple"
       storage=true
       seasonalstorage=false
   elseif storage_type =="none"
       storage=false
       seasonalstorage=false
  else
      storage=false
      seasonalstorage=false
      warn("String indicating `storage_type` not identified as 'none', 'seasonal' or 'simple' → no storage")
   end
  #The limit_dir is organized as two dictionaries in each other: limit_dir[impact][carrier]='impact/carrier' The first dictionary has the keys of the impacts, the second level dictionary has the keys of the carriers and value of the limit per carrier
  limit_emission=get_limit_dir(limit_emission)
  #Setup the config file based on the data input and
  config=set_config_cep(opt_data; descriptor=descriptor, limit_emission=limit_emission, lost_load_cost=lost_load_cost, lost_emission_cost=lost_emission_cost, infrastructure=infrastructure, demand=demand, non_dispatchable_generation=non_dispatchable_generation, dispatchable_generation=dispatchable_generation, storage=storage, conversion=conversion, seasonalstorage=seasonalstorage, transmission=transmission, scale=scale, print_flag=print_flag, optimizer_config=optimizer_config, round_sigdigits=round_sigdigits, region=opt_data.region, time_series=Dict{String,Any}("years" => ts_data.years, "K" => ts_data.K, "T"=> ts_data.T, "config" => time_series_config, "weights"=>ts_data.weights, "delta_t"=>ts_data.delta_t))
  #Run the optimization problem
  run_opt(ts_data, opt_data, config, optimizer)
end # run_opt


text:
"""
    Blade{T<:Number}
Type representing the exterior product of a set of orthogonal basis vectors, {eₖ}.
λ (e₁ ∧ e₂ ∧ ... ∧ eₙ)
Every multivector can be written as a linear combination of Blades, making this
type the building block for higher grade objects. The actual representation that
expresses the mathematical properties of the Blade is the rep::BitVector. This
representation contains 1's at the indices corresponding to the basis vectors.
eg. e1e2e3 = [1,1,1] e1e3 = [1,0,1], e3 = [0,0,1]
"""
struct Blade{T<:Number}
    val::T
    rep::BitVector
    Blade{T}() where {T} = new(undef,falses(0))
    Blade{T}(val::T,rep::BitVector) where {T<:Number}= new(val,rep)
end
Blade(v::T,rep::BitVector) where {T<:Number} = Blade{T}(v,rep)
Blade{T}(x::Number) where {T<:Number} = Blade{T}(x,falses(0))
Blade{T}(b::Blade) where {T<:Number} = Blade{T}(b.val,b.rep)

Blade(x::Number, label::Integer) = (rep=falses(label);rep[label]=1;Blade(x,rep))
#const I = Blade(true,trues()) #TODO:figure out implementation for unit pseudoscalar

Base.promote_rule(::Type{Blade{T}}, ::Type{S}) where {T<:Number,S<:Number} = Blade{promote_type(T,S)}
Base.promote_rule(::Type{Blade{T}}, ::Type{Blade{S}}) where {T<:Number,S<:Number} = Blade{promote_type(T,S)}
"""
    MultiVector <: AbstractSet{Blade}
Represents a MultiVector as a set of component Blades.
"""
struct MultiVector <: AbstractSet{Blade}
    dict::Dict{Blade,Nothing}
    MultiVector()= new(Dict{Blade,Nothing}())
end
MultiVector(itr) = union!(MultiVector(),itr)

Base.isempty(m::MultiVector) = isempty(m.dict)
Base.length(m::MultiVector) = length(m.dict)
Base.in(x::Blade, m::MultiVector) = haskey(m.dict,x)
Base.push!(m::MultiVector,x::Blade) = (m.dict[x]=nothing;m)

########################
###   Conversions   ####

Base.:(==)(q::Blade,r::Blade) = (q.val == r.val) & (q.rep == r.rep)
Base.:(==)(q::Blade,x::Real) = (q.val == x) & (isempty(q.rep))

################
#     Math     #
################
"""
    signature(x::Blade, y::Blade)
Determines the number of swaps that occur in a geometric product. Returns either 1 or -1.
An even number of swaps produces +1, and odd number results in -1. Captures the
antisymmetry under interchange of basis vectors. e1∧e2 = -e2∧e1
"""
function signature(x, y)
    xv = x >> -1
    yv = cumsum(y)
    return iseven(sum(xv .* yv)) ? 1.0 : -1.0
end
function padfalse!(x::BitVector,y::BitVector)
    if length(x) == length(y)
        return
    elseif length(x)>length(y)
        append!(y,falses(length(x)-length(y)))
    else
        append!(x,falses(length(y)-length(x)))
    end
end

# Blade product is binary XOR of the int representation.
function Base.:*(x::Blade, y::Blade)
    xp = x.rep
    yp = y.rep
    padfalse!(xp,yp)
    Blade(signature(xp,yp)*x.val*y.val, xor.(xp, yp))
end

Base.:*(x::Blade, y::Number) = Blade(x.val*y, x.rep)
Base.:*(x::Number, y::Blade) = Blade(x*y.val, y.rep)

@inline function Base.:+(b1::Blade,b2::Blade)
    if b1.rep == b2.rep
        return Blade(b1.val+b2.val, b1.rep)
    else
        return MultiVector([b1,b2])
    end
end
# Base.:+(x::Blade, y::Number) = Blade(x.val*y, x.rep)
# Base.:+(x::Number, y::Blade) = Blade(x*y.val, y.rep)

###################
# Representations #
###################
""" Parse an expression as a string of basis vectors, create blades of those
vectors"""
macro basis(x...)
    ex = Expr(:block)
    if length(x)==1 && isa(x[1],Expr)
        @assert x[1].head === :tuple "expected list of vectors."
        x = x[1].args
    end
    for (i,s) in enumerate(x)
        ex_aux = quote
            $s = Blade(1.0, $i)
        end
        append!(ex.args,ex_aux.args)
    end
    return esc(ex)
end

################
#    Display   #
################
function Base.repr(bx::Blade)
    indices = findall(bx.rep)
    str = "e" # Using this for basis vectors for now.
    if isempty(indices)
        return repr(bx.val)
    else
        return "( "*string(bx.val)*" ) ⋅ "*str*join(string.(indices),"e")
    end
end

function Base.show(io::IO, bx::Blade)
    print(io, repr(bx))
end


text:
using DelimitedFiles
using MPI

function runtrajs!(ntraj::Int64,nsteps::Int64,dt::Float64,t0::Float64,
                   x0::Vector{Float64},system::AbstractThermostattedSystem,
                   model::MixedModel,integrator::StochasticEuler,
                   mcvb::MCVBCallback)

  # run trajectories for evaluating gradients
  for traj in 1:ntraj
    system.x = copy(x0)
    system.t = t0
    initialize!(mcvb)
    runtraj!(nsteps,dt,integrator,mcvb)
    #println("Traj: $traj, ",system.x," ",mcvb.Afunc(system,model,1500,dt))
  end

end

"""
Train forces using stochastic gradient descent method
"""
# TODO add IO capabilities
function mcvbtrainsgd!(nepochs::Int64,ntraj::Int64,nsteps::Int64,dt::Float64,
                       t0::Float64,x0::Vector{Float64},lrf::Float64,lrv::Float64,
                       system::AbstractThermostattedSystem,model::MixedModel,
                       integrator::StochasticEuler,mcvb::MCVBCallback;
#                       varlr::Bool=false,lrscale::Float64=1.0,lrscaleevery::Int64=1000,
                       restart=true,printevery::Int64=1,Ffile::String="coeffF",
                       Vfile::String="coeffV")

  # start optimization
  for epoch in 1:nepochs

    # print update
    println("# Epoch: $epoch")
    flush(stdout)

    # zero out the averages
    initializeavgs!(mcvb)

    # run trajectories
    runtrajs!(ntraj,nsteps,dt,t0,x0,system,model,integrator,mcvb)

    # average them ove rthe number of trajectories
    average!(ntraj,mcvb)

    # TODO make nicer with IO
    # write average quantities to stdout
    println(mcvb.dkl," ",mcvb.aval)
    flush(stdout)

    # update coefficients for forces
    updateparams!(lrf,mcvb.gradf,model.potentials[end])

    # update coefficients for valuebaseline
    updateparams!(lrv,mcvb.gradv,mcvb.vbl)

    # TODO file system handling and outpout
    if restart
      writedlm("running_coeffF.txt",model.potentials[end].theta)
      writedlm("running_coeffV.txt",mcvb.vbl.theta)
    end

    if epoch % printevery == 0
      # write updated coefficients
      writedlm(Ffile*"_$epoch.txt",model.potentials[end].theta)
      writedlm(Vfile*"_$epoch.txt",model.potentials[end].theta)
    end

#    # update learning rate
#    if varlr
#      if epoch % lrscaleevery
#        lrf *= lrscale
#        lrv *= lrscale
#      end
#    end

  end

end

"""
Train forces using stochastic gradient descent method with parallelization
"""
# TODO add IO capabilities
function mcvbtrainsgdpar!(nepochs::Int64,totntraj::Int64,nsteps::Int64,dt::Float64,
                          t0::Float64,x0::Vector{Float64},lrf::Float64,lrv::Float64,
                          system::AbstractThermostattedSystem,model::MixedModel,
                          integrator::StochasticEuler,mcvb::MCVBCallback;
                          restart=true,printevery::Int64=1,Ffile::String="coeffF",
                          Vfile::String="coeffV")

  comm::MPI.Comm = MPI.COMM_WORLD
  rank::Int64 = MPI.Comm_rank(comm)
  nprocs::Int64 = MPI.Comm_size(comm)

  # calculate number of trajectories for this processor
  ntraj::Int64 = floor(Int64,totntraj/nprocs)
  extra::Int64 = totntraj % nprocs
  if (nprocs - rank - 1) < extra
    ntraj += 1
  end

  if rank == 0
    gradf_buf::Vector{Float64} = zeros(Float64, size(mcvb.gradf))
    gradv_buf::Vector{Float64} = zeros(Float64, size(mcvb.gradv))
  end

  # start optimization
  for epoch in 1:nepochs

    # print update
    if rank == 0
      println("# Epoch: $epoch")
      flush(stdout)
    end

    # zero out the averages
    initializeavgs!(mcvb)

    # run trajectories and average gradients
    runtrajs!(ntraj,nsteps,dt,t0,x0,system,model,integrator,mcvb)

    # average the gradients
    if rank == 0

      # receive information from each sub process and add to average
      for src in 1:(nprocs-1)
        MPI.Recv!(gradf_buf,src,100*src+0,comm)
        mcvb.gradf += gradf_buf
        MPI.Recv!(gradv_buf,src,100*src+1,comm)
        mcvb.gradv += gradv_buf
        mcvb.dkl += MPI.recv(src,100*src+2,comm)[1]
        mcvb.aval += MPI.recv(src,100*src+3,comm)[1]
      end

      # average them over the number of trajectories
      average!(totntraj,mcvb)

      # TODO make nicer with IO
      # write average quantities to stdout
      println(mcvb.dkl," ",mcvb.aval)
      flush(stdout)

      # update coefficients for forces
      updateparams!(lrf,mcvb.gradf,model.potentials[end])

      # update coefficients for valuebaseline
      updateparams!(lrv,mcvb.gradv,mcvb.vbl)

      # TODO file system handling and outpout
      if restart
        writedlm("running_coeffF.txt",model.potentials[end].theta)
        writedlm("running_coeffV.txt",mcvb.vbl.theta)
      end

      if epoch % printevery == 0
        # write updated coefficients
        writedlm(Ffile*"_$epoch.txt",model.potentials[end].theta)
        writedlm(Vfile*"_$epoch.txt",mcvb.vbl.theta)
      end

    else

      # send gradients to main processor
      MPI.Send(mcvb.gradf,0,100*rank+0,comm)
      MPI.Send(mcvb.gradv,0,100*rank+1,comm)
      MPI.send(mcvb.dkl,0,100*rank+2,comm)
      MPI.send(mcvb.aval,0,100*rank+3,comm)

    end

    # broadcast the updated parameters
    MPI.Barrier(comm)
    MPI.Bcast!(model.potentials[end].theta,0,comm)
    MPI.Bcast!(mcvb.vbl.theta,0,comm)
    MPI.Barrier(comm)

    # change the parameters for the model on the other processors
    if rank != 0
      setparams!(model.potentials[end])
      setparams!(mcvb.vbl)
    end

    # synchronize
    MPI.Barrier(comm)

#    # update learning rate
#    if varlr
#      if epoch % lrscaleevery
#        lrf *= lrscale
#        lrv *= lrscale
#      end
#    end

  end

end

"""
Train forces using ADAM optimizer
"""
# TODO add IO capabilities
function mcvbtrainadam!(nepochs::Int64,ntraj::Int64,nsteps::Int64,dt::Float64,
                        t0::Float64,x0::Vector{Float64},lrf::Float64,lrv::Float64,
                        beta1f::Float64,beta2f::Float64,beta1v::Float64,beta2v::Float64,eps::Float64,
                        system::AbstractThermostattedSystem,model::MixedModel,
                        integrator::StochasticEuler,mcvb::MCVBCallback;
                        restart=true,printevery::Int64=1,Ffile::String="coeffF",
                        Vfile::String="coeffV")

  # initialize extra vectors for adam algorithm
  mtf::Vector{Float64}   = zeros(Float64, size(mcvb.gradf))
  vtf::Vector{Float64}   = zeros(Float64, size(mcvb.gradf))
  mtv::Vector{Float64}   = zeros(Float64, size(mcvb.gradv))
  vtv::Vector{Float64}   = zeros(Float64, size(mcvb.gradv))
  t::Int64 = 1

  # start optimization
  for epoch in 1:nepochs

    # print update
    println("# Epoch: $epoch")
    flush(stdout)

    # zero out the averages
    initializeavgs!(mcvb)

    # run trajectories
    runtrajs!(ntraj,nsteps,dt,t0,x0,system,model,integrator,mcvb)

    # average them ove rthe number of trajectories
    average!(ntraj,mcvb)

    # TODO make nicer with IO
    # write average quantities to stdout
    println(mcvb.dkl," ",mcvb.aval)
    flush(stdout)

    # update bias corrections
    mtf .= (beta1f .* mtf) .+ ((1-beta1f) .* mcvb.gradf)
    vtf .= (beta2f .* vtf) .+ ((1-beta2f) .* mcvb.gradf.^2)
    mtv .= (beta1v .* mtv) .+ ((1-beta1v) .* mcvb.gradv)
    vtv .= (beta2v .* vtv) .+ ((1-beta2v) .* mcvb.gradv.^2)

    # add bias correction to gradients
    alphaf::Float64 = lrf * sqrt(1-beta2f^t) / (1-beta1f^t)
    alphav::Float64 = lrv * sqrt(1-beta2v^t) / (1-beta1v^t)
    mcvb.gradf .=  mtf ./ (sqrt.(vtf) .+ eps)
    mcvb.gradv .=  mtv ./ (sqrt.(vtv) .+ eps)
    t += 1

    # update coefficients for forces
    updateparams!(alphaf,mcvb.gradf,model.potentials[end])

    # update coefficients for valuebaseline
    updateparams!(alphav,mcvb.gradv,mcvb.vbl)

    # TODO file system handling and outpout
    if restart
      writedlm("running_coeffF.txt",model.potentials[end].theta)
      writedlm("running_coeffV.txt",mcvb.vbl.theta)
    end

    if epoch % printevery == 0
      # write updated coefficients
      writedlm(Ffile*"_$epoch.txt",model.potentials[end].theta)
      writedlm(Vfile*"_$epoch.txt",mcvb.vbl.theta)
    end

#    # update learning rate
#    if varlr
#      if epoch % lrscaleevery
#        lrf *= lrscale
#        lrv *= lrscale
#      end
#    end

  end

end

"""
Train forces using stochastic gradient descent method with parallelization
"""
# TODO add IO capabilities
function mcvbtrainadampar!(nepochs::Int64,totntraj::Int64,nsteps::Int64,dt::Float64,
                          t0::Float64,x0::Vector{Float64},lrf::Float64,lrv::Float64,
                          beta1f::Float64,beta2f::Float64,beta1v::Float64,beta2v::Float64,eps::Float64,
                          system::AbstractThermostattedSystem,model::MixedModel,
                          integrator::StochasticEuler,mcvb::MCVBCallback;
                          restart=true,printevery::Int64=1,Ffile::String="coeffF",
                          Vfile::String="coeffV")

  comm::MPI.Comm = MPI.COMM_WORLD
  rank::Int64 = MPI.Comm_rank(comm)
  nprocs::Int64 = MPI.Comm_size(comm)

  # calculate number of trajectories for this processor
  ntraj::Int64 = floor(Int64,totntraj/nprocs)
  extra::Int64 = totntraj % nprocs
  if (nprocs - rank - 1) < extra
    ntraj += 1
  end

  if rank == 0
    gradf_buf::Vector{Float64} = zeros(Float64, size(mcvb.gradf))
    gradv_buf::Vector{Float64} = zeros(Float64, size(mcvb.gradv))
    mtf::Vector{Float64} = zeros(Float64, size(mcvb.gradf))
    vtf::Vector{Float64} = zeros(Float64, size(mcvb.gradf))
    mtv::Vector{Float64} = zeros(Float64, size(mcvb.gradv))
    vtv::Vector{Float64} = zeros(Float64, size(mcvb.gradv))
    t::Int64 = 1
  end

  # start optimization
  for epoch in 1:nepochs

    # print update
    if rank == 0
      println("# Epoch: $epoch")
      flush(stdout)
    end

    # zero out the averages
    initializeavgs!(mcvb)

    # run trajectories and average gradients
    runtrajs!(ntraj,nsteps,dt,t0,x0,system,model,integrator,mcvb)

    # average the gradients
    if rank == 0

      # receive information from each sub process and add to average
      for src in 1:(nprocs-1)
        MPI.Recv!(gradf_buf,src,100*src+0,comm)
        mcvb.gradf += gradf_buf
        MPI.Recv!(gradv_buf,src,100*src+1,comm)
        mcvb.gradv += gradv_buf
        mcvb.dkl += MPI.recv(src,100*src+2,comm)[1]
        mcvb.aval += MPI.recv(src,100*src+3,comm)[1]
      end

      # average them over the number of trajectories
      average!(totntraj,mcvb)

      # TODO make nicer with IO
      # write average quantities to stdout
      println(mcvb.dkl," ",mcvb.aval," ")
      flush(stdout)

      # update bias corrections
      mtf .= (beta1f .* mtf) .+ ((1-beta1f) .* mcvb.gradf)
      vtf .= (beta2f .* vtf) .+ ((1-beta2f) .* mcvb.gradf.^2)
      mtv .= (beta1v .* mtv) .+ ((1-beta1v) .* mcvb.gradv)
      vtv .= (beta2v .* vtv) .+ ((1-beta2v) .* mcvb.gradv.^2)

      # add bias correction to gradients
      alphaf::Float64 = lrf * sqrt(1-beta2f^t) / (1-beta1f^t)
      alphav::Float64 = lrv * sqrt(1-beta2v^t) / (1-beta1v^t)
      #println(alphaf," ",alphav)
      println("grads before (f,v): ",sqrt(sum(mcvb.gradf.^2))," ",sqrt(sum(mcvb.gradv.^2)))
      gradfscale::Float64 = sqrt(sum(vtf)) + eps
      gradvscale::Float64 = sqrt(sum(vtv)) + eps
      mcvb.gradf .= gradfscale .* mtf ./ (sqrt.(vtf) .+ eps)
      mcvb.gradv .= gradvscale .* mtv ./ (sqrt.(vtv) .+ eps)
      println("mean things (f,v): ",sqrt(sum(mtf.^2))," ",sqrt(sum(mtv.^2)))
      println("variance things (f,v): ",sqrt(sum(vtf))," ",sqrt(sum(vtv)))
      println("grads after (f,v): ",sqrt(sum(mcvb.gradf.^2))," ",sqrt(sum(mcvb.gradv.^2)))
      t += 1
      #alphaf::Float64 = lrf
      #alphav::Float64 = lrv

      #println(sqrt(sum(mcvb.gradf.^2))," ",sqrt(sum(mcvb.gradv.^2)))
      #flush(stdout)

      # update coefficients for forces
      updateparams!(alphaf,mcvb.gradf,model.potentials[end])

      # update coefficients for valuebaseline
      updateparams!(alphav,mcvb.gradv,mcvb.vbl)

      # TODO file system handling and outpout
      if restart
        writedlm("running_coeffF.txt",model.potentials[end].theta)
        writedlm("running_coeffV.txt",mcvb.vbl.theta)
      end

      if epoch % printevery == 0
        # write updated coefficients
        writedlm(Ffile*"_$epoch.txt",model.potentials[end].theta)
        writedlm(Vfile*"_$epoch.txt",mcvb.vbl.theta)
      end

    else

      # send gradients to main processor
      MPI.Send(mcvb.gradf,0,100*rank+0,comm)
      MPI.Send(mcvb.gradv,0,100*rank+1,comm)
      MPI.send(mcvb.dkl,0,100*rank+2,comm)
      MPI.send(mcvb.aval,0,100*rank+3,comm)

    end

    # broadcast the updated parameters
    MPI.Barrier(comm)
    MPI.Bcast!(model.potentials[end].theta,0,comm)
    MPI.Bcast!(mcvb.vbl.theta,0,comm)
    MPI.Barrier(comm)

    # change the parameters for the model on the other processors
    if rank != 0
      setparams!(model.potentials[end])
      setparams!(mcvb.vbl)
    end

    # synchronize
    MPI.Barrier(comm)

  end

end


text:
# <editor-fold desc= create other elements of model"

# XXX create variables and capacity constraints for trade variables
function createTradeVarCns!(partTrd::OthPart,anyM::anyModel)
	for type in (:Buy, :Sell)
		trdPrc_sym = Symbol(:trd,type,:Prc)
		trd_sym = Symbol(:trd,type)
		if trdPrc_sym in keys(partTrd.par) && :C in namesSym(partTrd.par[trdPrc_sym].data)

			# <editor-fold desc="create trade variables"
			c_arr = unique(partTrd.par[trdPrc_sym].data[!,:C])

			# create dataframe with all potential entries for trade/sell variable
			var_df = createPotDisp(c_arr,anyM)

			# match all potential variables with defined prices
			var_df = matchSetParameter(var_df,partTrd.par[trdPrc_sym],anyM.sets)[!,Not(:val)]

			var_df = createVar(var_df,string(:trd,type),getUpBound(var_df,anyM.options.bound.disp / anyM.options.scaFac.dispTrd,anyM.supTs,anyM.sets[:Ts]),anyM.optModel,anyM.lock,anyM.sets, scaFac = anyM.options.scaFac.dispTrd)
			partTrd.var[trd_sym] = orderDf(var_df)
			produceMessage(anyM.options,anyM.report, 3," - Created variables for $(type == :Buy ? "buying" : "selling") carriers")
			# </editor-fold>

			# <editor-fold desc="create capacity constraint on variable"
			trdCap_sym = Symbol(trd_sym,:Cap)
			if trdCap_sym in keys(partTrd.par)
				cns_df = matchSetParameter(var_df,partTrd.par[trdCap_sym],anyM.sets,newCol = :cap)
				sca_arr = getResize(cns_df,anyM.sets[:Ts],anyM.supTs)
				cns_df[!,:cap] = cns_df[!,:cap] .* sca_arr

				# prepare, scale and create constraints
				cns_df[!,:cnsExpr] = map(x -> x.var - x.cap, eachrow(cns_df))
				scaleCnsExpr!(cns_df,anyM.options.coefRng,anyM.options.checkRng)
				partTrd.cns[trdCap_sym] = createCns(cnsCont(cns_df,:smaller),anyM.optModel)

				produceMessage(anyM.options,anyM.report, 3," - Created capacity restrictions for $(type == :Buy ? "buying" : "selling") carriers")
			end

			# </editor-fold>
		end
		produceMessage(anyM.options,anyM.report, 2," - Created variables and constraints for $(type == :Buy ? "buying" : "selling") carriers")
	end
	produceMessage(anyM.options,anyM.report, 1," - Created variables and constraints for trade")
end

# XXX create all energy balances (and curtailment variables if required)
function createEnergyBal!(techSym_arr::Array{Symbol,1},anyM::anyModel)

	partBal = anyM.parts.bal
	c_arr = filter(x -> x != 0,getfield.(values(anyM.sets[:C].nodes),:idx))
	allDim_df = createPotDisp(c_arr,anyM)
	bal_tup = (:C,:Ts_dis)
	agg_arr = [:Ts_dis, :R_dis, :C]

	# <editor-fold desc="create potential curtailment and loss loss load variables

	for varType in (:crt,:lss)
		# get defined entries
		var_df = DataFrame()
		for par in intersect(keys(partBal.par),vcat(varType == :crt ? :costCrt : :costLss, Symbol.(varType,[:Up,:Low,:Fix])...))
			append!(var_df,matchSetParameter(allDim_df,partBal.par[par],anyM.sets)[!,Not(:val)])
		end

		# obtain upper bound for variables and create them
		if !isempty(var_df)
			partBal.var[varType] = orderDf(createVar(var_df,string(varType),getUpBound(var_df,anyM.options.bound.disp / anyM.options.scaFac.dispTrd,anyM.supTs,anyM.sets[:Ts]),anyM.optModel,anyM.lock,anyM.sets, scaFac = anyM.options.scaFac.dispTrd))
		end
	end
	# </editor-fold>

	# <editor-fold desc="create actual balance"

	# finds all carriers that require an energy balance (not required for carriers that can only be shifted (temporal or spatial), e.g. that only have storage or exchange defined for them)
	relC_arr = Array{Int,1}()

	# if demand does not specify a carrier, it is assumed all carriers are relevant
	if !isempty(anyM.parts.bal.par[:dem].data)
		if :C in namesSym(anyM.parts.bal.par[:dem].data)
			append!(relC_arr,unique(anyM.parts.bal.par[:dem].data[!,:C]))
		else
			append!(relC_arr,c_arr)
		end
	end

	if :crt in keys(anyM.parts.bal.var) append!(relC_arr,unique(anyM.parts.bal.var[:crt][!,:C])) end
	if :lss in keys(anyM.parts.bal.var) append!(relC_arr,unique(anyM.parts.bal.var[:lss][!,:C])) end
	if :trdSell in keys(anyM.parts.trd.var) append!(relC_arr,unique(anyM.parts.trd.var[:trdSell][!,:C])) end
	if :trdBuy in keys(anyM.parts.trd.var) append!(relC_arr,unique(anyM.parts.trd.var[:trdBuy][!,:C])) end

	# add carriers beings generated or used
	append!(relC_arr,union(union(map(x -> anyM.parts.tech[x].carrier |> (y -> map(z -> getfield(y,z),intersect(keys(y),(:gen,:use)))),techSym_arr)...)...))
	relC_arr = unique(relC_arr)

	# create object to write constraint data too
	cns_arr = Array{Pair{Symbol,cnsCont}}(undef,length(relC_arr))
	itrC_arr = collect(enumerate(relC_arr))

	@threads for (idx,c) in itrC_arr

		subC_arr = unique([c,getDescendants(c,anyM.sets[:C],true)...])
		cRes_tup = anyM.cInfo[c] |> (x -> (Ts_dis = x.tsDis, R_dis = x.rDis, C = anyM.sets[:C].nodes[c].lvl))

		# XXX add demand and size it
		cns_df = matchSetParameter(filter(x -> x.C == c,allDim_df),partBal.par[:dem],anyM.sets)
		cns_df[!,:dem] = cns_df[!,:val] .* getResize(cns_df,anyM.sets[:Ts],anyM.supTs)
		select!(cns_df,Not(:val))

		# XXX get relevant variables
		src_df = cns_df[!,Not([:Ts_disSup,:dem])]

		# add tech variables
		cns_df[!,:techVar] = getTechEnerBal(c,subC_arr,src_df,anyM.parts.tech,anyM.cInfo,anyM.sets)

		# add curtailment variables
		for varType in (:crt,:lss)
			if varType in keys(partBal.var)
				cns_df[!,Symbol(varType,:Var)] = filterCarrier(partBal.var[varType],subC_arr) |> (x -> aggUniVar(x,src_df,agg_arr, cRes_tup,anyM.sets))
			else
				cns_df[!,Symbol(varType,:Var)] .= AffExpr()
			end
		end

		# add trade variables
		if !isempty(anyM.parts.trd.var)
			cns_df[!,:trdVar] = sum([filterCarrier(anyM.parts.trd.var[trd],subC_arr) |> (x -> aggUniVar(x,src_df,agg_arr,cRes_tup,anyM.sets) |> (y -> trd != :trdSell ? y : -1.0 * y)) for trd in keys(anyM.parts.trd.var)])
		else
			cns_df[!,:trdVar] .= AffExpr()
		end

		# add exchange variables
		if !isempty(anyM.parts.exc.var)
			excVarTo_df = filterCarrier(anyM.parts.exc.var[:exc],subC_arr)
			excVarFrom_df = convertExcCol(copy(excVarTo_df))

			# get loss values and apply them to variables
			excVarFrom_df = getExcLosses(excVarFrom_df,anyM.parts.exc.par,anyM.sets)
			excVarFrom_df[!,:var] = excVarFrom_df[!,:var] .* (1.0 .- excVarFrom_df[!,:loss])
			select!(excVarFrom_df,Not(:loss))

			balTo_tup, balFrom_tup = [tuple(replace(collect(bal_tup),:R_dis => x)...) for x in [:R_to, :R_from]]

			excFrom_arr = aggUniVar(convertExcCol(excVarFrom_df),rename(src_df,:R_dis => :R_to),[:Ts_dis,:R_to,:C],(Ts_dis = cRes_tup[1], R_to = cRes_tup[2], C = cRes_tup[3]),anyM.sets)
			excTo_arr  = aggUniVar(excVarTo_df,rename(src_df,:R_dis => :R_from),[:Ts_dis,:R_from,:C],(Ts_dis = cRes_tup[1], R_from = cRes_tup[2], C = cRes_tup[3]),anyM.sets)

			cns_df[!,:excVar] =  excFrom_arr .- excTo_arr
		else
			cns_df[!,:excVar] .= AffExpr()
		end

		# prepare, scale and save constraints to dictionary
		c_str = Symbol(anyM.sets[:C].nodes[c].val)
		cns_df[!,:cnsExpr] = map(x -> x.techVar + x.excVar + x.trdVar + x.lssVar - x.dem - x.crtVar, eachrow(cns_df))
		cns_df = orderDf(cns_df[!,[intCol(cns_df)...,:cnsExpr]])
		filter!(x -> x.cnsExpr != AffExpr(),cns_df)
		scaleCnsExpr!(cns_df,anyM.options.coefRng,anyM.options.checkRng)
		cns_arr[idx] = Symbol(c_str) => cnsCont(cns_df,anyM.cInfo[c].eq ? :equal : :greater)

		produceMessage(anyM.options,anyM.report, 2," - Prepared energy balance for $(c_str)")
	end

	# loops over stored constraints outside of threaded loop to create actual jump constraints
	for cns in cns_arr
		partBal.cns[cns[1]] = createCns(cns[2],anyM.optModel)
	end

	produceMessage(anyM.options,anyM.report, 1," - Created energy balances for all carriers")
	# </editor-fold>
end

# XXX aggregate all technology variables for energy balance
function getTechEnerBal(cBal_int::Int,subC_arr::Array{Int,1},src_df::DataFrame,tech_dic::Dict{Symbol,TechPart},
																				cInfo_dic::Dict{Int,NamedTuple{(:tsDis,:tsExp,:rDis,:rExp,:eq),Tuple{Int,Int,Int,Int,Bool}}},sets_dic::Dict{Symbol,Tree})
	techVar_arr = Array{Array{AffExpr,1}}(undef,length(subC_arr))

	# get temporal and spatial resolution for carrier being balanced
	cBalRes_tup = cInfo_dic[cBal_int] |> (x -> (x.tsDis, x.rDis))

	# loops over all carriers relevant for respective energy balance
	for (idx,c) in enumerate(subC_arr)

		# gets technologies relevant for respective filterCarrier
		relTech_arr = getRelTech(c,tech_dic,sets_dic[:C])

		# leaves loop for carrier, if no relevant technologies could be obtained
		if isempty(relTech_arr)
			techVar_arr[idx]  = fill(AffExpr(),size(src_df,1))
			continue
		end

		# prepare loop over tech for c by creating empty dataframe and get temporal and spatial resolution for carrier being balanced
		allVar_df = DataFrame(Ts_dis = Int[], R_dis = Int[], var = AffExpr[])
		cRes_tup = cInfo_dic[c] |> (x -> (x.tsDis, x.rDis))

		for x in relTech_arr
			# gets resolution and adjusts add_df in case of an agggregated technology
			add_df = select(filter(r -> r.C == c,tech_dic[x[1]].var[x[2]]),[:Ts_dis,:R_dis,:var])
			if isempty(add_df) continue end
			tRes_tup = tech_dic[x[1]].disAgg ? (cRes_tup[1], tech_dic[x[1]].balLvl.exp[2]) : cRes_tup
			checkTechReso!(tRes_tup,cBalRes_tup,add_df,sets_dic)

			# adds sign to variables and adds them to overall dataframe
			add_df[!,:var] = add_df[!,:var] .* (x[2] in (:use,:stExtIn) ? -1.0 : 1.0)
			append!(allVar_df, add_df)
		end

		# returns empty expression if no variales could be obtained
		if isempty(allVar_df)
			techVar_arr[idx]  = fill(AffExpr(),size(src_df,1))
		else
			grpVar_df = combine(groupby(allVar_df, [:Ts_dis, :R_dis]), :var => (x -> sum(x)) => :var)
			techVar_arr[idx] = joinMissing(src_df,grpVar_df, [:Ts_dis, :R_dis], :left, Dict(:var => AffExpr()))[!,:var]
		end
	end

	return map(x -> sum(x),eachrow(hcat(techVar_arr...)))
end

# XXX create constarints that enforce any type of limit (Up/Low/Fix) on any type of variable
function createLimitCns!(partLim::OthPart,anyM::anyModel)

	parLim_arr = String.(collectKeys(keys(partLim.par)))
	techLim_arr = filter(x ->  any(map(y -> occursin(y,x),["Up","Low","Fix"])),parLim_arr)
	limVar_arr = map(x -> map(k -> Symbol(k[1]) => Symbol(k[2][1]), filter(z -> length(z[2]) == 2,map(y -> y => split(x,y),["Up","Low","Fix"])))[1], techLim_arr)
	varToPar_dic = Dict(y => getindex.(filter(z -> z[2] == y,limVar_arr),1) for y in unique(getindex.(limVar_arr,2)))

	# loop over all variables that are subject to any type of limit (except emissions)
	allKeys_arr = collect(keys(varToPar_dic))
	cns_dic = Dict{Symbol,cnsCont}()
	signLim_dic= Dict(:Up => :smaller, :Low => :greater, :Fix => :equal)

	@threads for va in allKeys_arr

		varToPart_dic = Dict(:exc => :exc, :crt => :bal,:trdSell => :trd, :trdBuy => :trd)

		# obtain all variables relevant for limits
		allVar_df = getAllVariables(va,anyM)

		# check if acutally any variables were obtained
		if isempty(allVar_df)
			lock(anyM.lock)
			push!(anyM.report,(2,"limit",string(va),"limits for variable provided, but none of these variables are actually created"))
			unlock(anyM.lock)
			continue
		end

		allLimit_df = DataFrame(var = AffExpr[])
		# XXX loop over respective type of limits to obtain data
		for lim in varToPar_dic[va]
			par_obj = copy(partLim.par[Symbol(va,lim)])
			if va in (:capaExc,:oprCapaExc) && :R_a in namesSym(par_obj.data) && :R_b in namesSym(par_obj.data)
				par_obj.data = vcat(par_obj.data,rename(par_obj.data,:R_a => :R_b,:R_b => :R_a))
			end
			agg_tup = tuple(intCol(par_obj.data)...)

			# aggregate search variables according to dimensions in limit parameter
			if isempty(agg_tup)
				grpVar_df = allVar_df
			else
				grpVar_df = combine(groupby(convertExcCol(allVar_df),collect(agg_tup)), :var => (x -> sum(x)) => :var)
			end

			# try to aggregate variables to limits directly provided via inputs
			limit_df = copy(par_obj.data)
			if size(limit_df,2) != 1
				limit_df[!,:var] = aggDivVar(grpVar_df, limit_df[!,Not(:val)], agg_tup, anyM.sets, aggFilt = agg_tup)
			else
				limit_df[!,:var] .= sum(grpVar_df[!,:var])
			end

			# gets provided limit parameters, that no variables could assigned to so far and tests if via inheritance any could be assigned
			mtcPar_arr, noMtcPar_arr  = findall(map(x -> x != AffExpr(),limit_df[!,:var])) |>  (x -> [x, setdiff(1:size(par_obj.data,1),x)])
			# removes entries with no parameter assigned from limits
			limit_df = limit_df[mtcPar_arr,:]

			if !isempty(noMtcPar_arr)
				# tries to inherit values to existing variables only for parameters without variables aggregated so far
				aggPar_obj = copy(par_obj,par_obj.data[noMtcPar_arr,:])
				aggPar_obj.data = matchSetParameter(grpVar_df[!,Not(:var)],aggPar_obj,anyM.sets, useNew = false)
				# again performs aggregation for inherited parameter data and merges if original limits
				aggLimit_df = copy(aggPar_obj.data)
				if !isempty(aggLimit_df)
					aggLimit_df[!,:var]  = aggDivVar(grpVar_df, aggLimit_df, agg_tup, anyM.sets, aggFilt = agg_tup)
					limit_df = vcat(limit_df,aggLimit_df)
				end
			end

			# merge limit constraint to other limits for the same variables
			limit_df = convertExcCol(rename(limit_df,:val => lim))
			join_arr = [intersect(intCol(allLimit_df),intCol(limit_df))...,:var]
			miss_arr = [intCol(allLimit_df),intCol(limit_df)] |> (y -> union(setdiff(y[1],y[2]), setdiff(y[2],y[1])))
			allLimit_df = joinMissing(allLimit_df, limit_df, join_arr, :outer, merge(Dict(z => 0 for z in miss_arr),Dict(:Up => nothing, :Low => nothing, :Fix => nothing)))
		end

		# XXX check for contradicting values
		colSet_dic = Dict(x => Symbol(split(string(x),"_")[1]) for x in intCol(allLimit_df))
		limitCol_arr = intersect(namesSym(allLimit_df),(:Fix,:Up,:Low))
		entr_int = size(allLimit_df,1)
		if :Low in limitCol_arr || :Up in limitCol_arr
			# ! errors

			# upper and lower limit contradicting each other
			if :Low in limitCol_arr && :Up in limitCol_arr
				for x in findall(replace(allLimit_df[!,:Low],nothing => 0.0) .> replace(allLimit_df[!,:Up],nothing => Inf))
					dim_str = join(map(y -> allLimit_df[x,y] == 0 ?  "" : string(y,": ",join(getUniName(allLimit_df[x,y], anyM.sets[colSet_dic[y]])," < ")),intCol(allLimit_df)),"; ")
					lock(anyM.lock)
					push!(anyM.report,(3,"limit",string(va),"contradicting values for upper and lower limit detected for: " * dim_str))
					unlock(anyM.lock)
				end
			end

			# fix and upper limit contradicting each other
			if :Fix in limitCol_arr && :Up in limitCol_arr
				for x in findall(replace(allLimit_df[!,:Fix],nothing => 0.0) .> (replace(allLimit_df[!,:Up],nothing => Inf) .+ 0.0001))
					dim_str = join(map(y -> allLimit_df[x,y] == 0 ?  "" : string(y,": ",join(getUniName(allLimit_df[x,y], anyM.sets[colSet_dic[y]])," < ")),intCol(allLimit_df)),"; ")
					lock(anyM.lock)
					push!(anyM.report,(3,"limit",string(va),"fixed limit exceeds upper limit for: " * dim_str))
					unlock(anyM.lock)
				end
			end

			# fix and lower limit contradicting each other
			if :Fix in limitCol_arr && :Low in limitCol_arr
				for x in findall(replace(allLimit_df[!,:Fix],nothing => Inf) .< replace(allLimit_df[!,:Low],nothing => 0.0) .- 0.0001)
					dim_str = join(map(y -> allLimit_df[x,y] == 0 ?  "" : string(y,": ",join(getUniName(allLimit_df[x,y], anyM.sets[colSet_dic[y]])," < ")),intCol(allLimit_df)),"; ")
					lock(anyM.lock)
					push!(anyM.report,(3,"limit",string(va),"fixed limit is smaller than lower limit for: " * dim_str))
					unlock(anyM.lock)
				end
			end

			# residual values already violate limits
			resiVal_arr = getfield.(allLimit_df[!,:var],:constant)
			if :Up in limitCol_arr
				for x in findall(resiVal_arr .>  replace(allLimit_df[!,:Up],nothing => Inf))
					dim_str = join(map(y -> allLimit_df[x,y] == 0 ?  "" : string(y,": ",join(getUniName(allLimit_df[x,y], anyM.sets[colSet_dic[y]])," < ")),intCol(allLimit_df)),"; ")
					lock(anyM.lock)
					push!(anyM.report,(3,"limit",string(va),"residual values already exceed the upper limit for: " * dim_str))
					unlock(anyM.lock)
				end
			end

			# ! warnings
			# upper or lower limit of zero
			if !isempty(limitCol_arr |> (y -> filter(x -> collect(x[y]) |> (z -> any(isnothing.(z)) ? false : any(z .== 0)),allLimit_df))) && va != :emission
				lock(anyM.lock)
				push!(anyM.report,(2,"limit",string(va),"upper or lower limit of zero detected, please consider to use fix or omit limit instead"))
				unlock(anyM.lock)
				entr_int = size(allLimit_df,1)
			end

			# value is fixed, but still a upper a lower limit is provided
			if :Fix in limitCol_arr
				if !isempty(filter(x -> x != :Fix, limitCol_arr) |> (z -> filter(x -> all([!isnothing(x.Fix),any(.!isnothing.(collect(x[z])))]) ,eachrow(allLimit_df))))
					lock(anyM.lock)
					push!(anyM.report,(2,"limit",string(va),"upper and/or lower limit detected, although variable is already fixed"))
					unlock(anyM.lock)
				end
			end
		end

		# if installed capacities differ depending on the direction, because residual values were defined and at the same time fixed limits on the installed capacity were provided
		# an error will occur, because a value cannot be fixed but and the same time differ by direction, this is detected here
		if :Fix in limitCol_arr && va == :capaExc
			# filters fixed exchange capacities and extracts residual values
			fix_df = select(filter(x -> x.Fix != nothing, allLimit_df) ,intCol(allLimit_df,:var))
			fix_df[!,:resi] .= getfield.(fix_df[!,:var],:constant)
			# joins together capacities in both directions
			joinA_df = rename(select(fix_df,Not([:var])),:resi => :resiA)
			joinB_df = rename(joinA_df,:resiA => :resiB)
			comp_df = innerjoin(joinA_df,joinB_df, on = intCol(joinA_df) .=> replace(intCol(joinB_df),:R_from => :R_to,:R_to => :R_from))
			# finds cases that lead to contradiction and reports on them
			contraExc_df = filter(x -> x.resiA != x.resiB && x.R_from > x.R_to,comp_df)
			for x in eachrow(contraExc_df)
				dim_str = join(map(y -> x[y] == 0 ?  "" : string(y,": ",join(getUniName(x[y], anyM.sets[colSet_dic[y]])," < ")),intCol(contraExc_df)),"; ")
				lock(anyM.lock)
				push!(anyM.report,(3,"limit",string(va),"for the exchange capacity '" * dim_str * "' residual capacites differ by direction but at the same time the installed capacity in both directions is fixed to the same value by capaExcFix, this is a contradiction and  would lead to an infeasible model"))
				unlock(anyM.lock)
			end
		end

		# XXX check for suspicious entries for capacity where limits are provided for the sum of capacity over several years
		if occursin("capa",string(va))
			if !(:Ts_disSup in namesSym(allLimit_df))
				lock(anyM.lock)
				push!(anyM.report,(2,"limit","capacity","capacity limits were provided without specificing the superordinate dispatch timestep, this means the sum of capacity over all years was limited instead of enforcing the same limit for each year
																												(see https://leonardgoeke.github.io/AnyMOD.jl/stable/parameter_list/#Limits-on-quantities-dispatched)"))
				unlock(anyM.lock)
			elseif 0 in unique(allLimit_df[!,:Ts_disSup])
				relEntr_df = filter(x -> x.Ts_disSup == 0, allLimit_df)
				if :Te in namesSym(relEntr_df)
					allTe_arr = unique(relEntr_df[!,:Te])
					for tInt in allTe_arr
						push!(anyM.report,(2,"limit","capacity","capacity limits were provided for $(string(techSym(tInt,anyM.sets[:Te]))) without specificing the superordinate dispatch timestep, this means the sum of capacity over all years was limited instead of enforcing the same limit for each year
																														(see https://leonardgoeke.github.io/AnyMOD.jl/stable/parameter_list/#Limits-on-quantities-dispatched)"))
					end
				else
					lock(anyM.lock)
					push!(anyM.report,(2,"limit","capacity","capacity limits were provided without specificing the superordinate dispatch timestep, this means the sum of capacity over all years was limited instead of enforcing the same limit for each year
																													(see https://leonardgoeke.github.io/AnyMOD.jl/stable/parameter_list/#Limits-on-quantities-dispatched)"))
					unlock(anyM.lock)
				end
			end
		end

		# XXX write constraint containers
		for lim in limitCol_arr
			# filter respective limits (low, fix or up) out of the entire dataframe
			relLim_df = filter(x -> !isnothing(x[lim]),allLimit_df[!,Not(filter(x -> x != lim,limitCol_arr))])
			relLim_df = filter(x -> x.var != AffExpr() || x.Fix != 0.0, relLim_df)
            if isempty(relLim_df) continue end
			rename!(relLim_df,lim => :Lim)

			# prepare, scale and save constraints to dictionary
			relLim_df[!,:cnsExpr] = map(x -> x.var - x.Lim, eachrow(relLim_df))
			relLim_df = orderDf(relLim_df[!,[intCol(relLim_df)...,:cnsExpr]])
			scaleCnsExpr!(relLim_df,anyM.options.coefRng,anyM.options.checkRng)
			cns_dic[Symbol(va,lim)] = cnsCont(relLim_df,signLim_dic[lim])

			produceMessage(anyM.options,anyM.report, 3," - Created constraints for $(lim == :Up ? "upper" : (lim == :Low ? "lower" : "fixed")) limit of variable $va")
		end
		typeLim_sym = va in (:emission,) ? "term" : "variable"
		produceMessage(anyM.options,anyM.report, 2," - Prepared constraints to limit $typeLim_sym $va")
	end

	# loops over stored constraints outside of threaded loop to create actual jump constraints
	for cnsSym in keys(cns_dic)
		partLim.cns[cnsSym] = createCns(cns_dic[cnsSym],anyM.optModel)
	end

	produceMessage(anyM.options,anyM.report, 1," - Created all limiting constraints")
end

# </editor-fold>

# <editor-fold desc= utility functions"

# XXX connect capacity and expansion variables
function createCapaCns!(part::TechPart,prepTech_dic::Dict{Symbol,NamedTuple},cns_dic::Dict{Symbol,cnsCont})
    for capaVar in filter(x -> occursin("capa",string(x)),keys(prepTech_dic))
        index_arr = intCol(part.var[capaVar])
        join_arr = part.type != :mature ? index_arr : filter(x -> x != :Ts_expSup,collect(index_arr))

        # joins corresponding capacity and expansion variables together
		expVar_sym = Symbol(replace(string(capaVar),"capa" => "exp"))
		if !(expVar_sym in keys(part.var)) continue end
        expVar_df = flatten(part.var[expVar_sym],:Ts_disSup)
        cns_df = rename(innerjoin(part.var[capaVar],combine(groupby(expVar_df,join_arr), :var => (x -> sum(x)) => :exp); on = join_arr),:var => :capa)

        # creates final constraint object
		cns_df[!,:cnsExpr] = map(x -> x.capa - x.capa.constant - x.exp,eachrow(cns_df))
		cns_dic[Symbol(capaVar)] = cnsCont(select(cns_df,Not([:capa,:exp])),:equal)
    end
end

# XXX adds column with JuMP variable to dataframe
function createVar(setData_df::DataFrame,name_str::String,upBd_fl::Union{Float64,Array{Float64,1}},optModel::Model,lock_::ReentrantLock,sets::Dict{Symbol,Tree}; scaFac::Float64 = 1.0, lowBd::Float64 = 0.0)
	# adds an upper bound to all variables if provided within the options
	#if isempty(setData_df) return DataFrame(var = AffExpr[]) end
	arr_boo = typeof(upBd_fl) <: Array
	if arr_boo
		info = VariableInfo.(!isnan(lowBd), lowBd, .!isnan.(upBd_fl), upBd_fl, false, NaN, false, NaN, false, false)
		var_obj = JuMP.build_variable.(error, info)
	else
		info = VariableInfo(!isnan(lowBd), lowBd, !isnan(upBd_fl), upBd_fl, false, NaN, false, NaN, false, false)
		var_obj = JuMP.build_variable(error, info)
	end

	# writes full name for each variable
	setData_df = orderDf(setData_df)
	dim_arr = map(x -> Symbol(split(String(x),"_")[1]), filter(r -> r != :id,intCol(setData_df)))
	dim_int = length(dim_arr)
	setData_df[!,:name] = string.(name_str,"[",map(x -> join(map(y -> sets[dim_arr[y]].nodes[x[y]].val,1:dim_int),", "),eachrow(setData_df)),"]")

	lock(lock_)
	if arr_boo
		setData_df[!,:var] = [AffExpr(0,JuMP.add_variable(optModel, nameItr[1], nameItr[2]) => scaFac) for nameItr in zip(var_obj,setData_df[!,:name])]
	else
		setData_df[!,:var] = [AffExpr(0,JuMP.add_variable(optModel, var_obj, nameItr) => scaFac) for nameItr in setData_df[!,:name]]
	end
	unlock(lock_)

	return setData_df[!,Not(:name)]
end

# XXX scales expressions in the dataframe to be within the range defined within options
function scaleCnsExpr!(cnsExpr_df::DataFrame,coefRng::NamedTuple{(:mat,:rhs),Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}},checkRng_fl::Float64)

	if isempty(cnsExpr_df) return end

	# filters variables with factor zero, because they will lead to scaling errors
	foreach(z -> foreach(y ->  delete!(z.terms,y), filter(x -> z.terms[x] == 0.0 ,collect(keys(z.terms)))), cnsExpr_df[!,:cnsExpr])

	if !all(isnan.(coefRng.rhs))
		# scale expression defining constraint so rhs coefficients are within desired range
		rhs_arr = abs.(getfield.(cnsExpr_df[!,:cnsExpr],:constant))
		findall(rhs_arr .!= 0.0) |> (x -> cnsExpr_df[x,:cnsExpr] = scaleRng(cnsExpr_df[x,:cnsExpr],rhs_arr[x],coefRng.rhs, true))
	end

	if !all(isnan.(coefRng.mat))
		# scale expression defining constraint so matrix coefficients are within desired range
		matRng_arr = map(x -> abs.(values(x.terms)) |> (y -> isempty(y) ? (0.0,0.0) : (minimum(y),maximum(y))), cnsExpr_df[!,:cnsExpr])
		findall(map(x -> x != (0.0,0.0),matRng_arr)) |> (x -> cnsExpr_df[x,:cnsExpr] = scaleRng(cnsExpr_df[x,:cnsExpr],matRng_arr[x],coefRng.mat,false))
	end

	if !isnan(checkRng_fl)
		checkExprRng(cnsExpr_df[:,:cnsExpr],checkRng_fl)
	end
end

# XXX used to perform scaling of expression array based on range of coefficients provided
function scaleRng(expr_arr::Array{AffExpr,1},rng_arr::Array,rng_tup::Tuple{Float64,Float64}, rhs_boo::Bool)
	scaRel_arr = rhs_boo ? union(findall(rng_arr .< rng_tup[1]), findall(rng_arr .> rng_tup[2])) : union(findall(getindex.(rng_arr,1) .< rng_tup[1]), findall(getindex.(rng_arr,2) .> rng_tup[2]))
	if !isempty(scaRel_arr)
		expr_arr[scaRel_arr] = map(x -> x[1] < rng_tup[1] ? rng_tup[1]/x[1] : rng_tup[2]/x[rhs_boo ? 1 : 2], rng_arr[scaRel_arr]) .* expr_arr[scaRel_arr]
	end
	return expr_arr
end

# XXX check range of coefficients in expressions within input array
function checkExprRng(expr_arr::Array{AffExpr,1},rngThres_fl::Float64)
	# obtains range of coefficients for matrix and rhs
	matRng_arr = map(x -> abs.(values(x.terms)) |> (y -> isempty(y) ? (0.0,0.0) : (minimum(y),maximum(y))), expr_arr)
	rhs_arr = abs.(getfield.(expr_arr,:constant))
	both_arr = max.(getindex.(matRng_arr,2),replace(rhs_arr,0.0 => -Inf)) ./ min.(getindex.(matRng_arr,1),replace(rhs_arr,0.0 => Inf))

	# filters rows where reange of coefficients is above threshold
	aboveThres_arr = findall(both_arr .> rngThres_fl)

	for expr in expr_arr[aboveThres_arr]
		println(expr)
	end
end

# XXX creates an actual jump constraint based on the constraint container provided
function createCns(cnsCont_obj::cnsCont,optModel::Model)
	cns_df = cnsCont_obj.data
	if cnsCont_obj.sign == :equal
		cns_df[!,:cns] = map(x -> @constraint(optModel, x.cnsExpr == 0),eachrow(cns_df))
	elseif cnsCont_obj.sign == :greater
		cns_df[!,:cns] = map(x -> @constraint(optModel, x.cnsExpr >= 0),eachrow(cns_df))
	elseif cnsCont_obj.sign == :smaller
		cns_df[!,:cns] = map(x -> @constraint(optModel, x.cnsExpr <= 0),eachrow(cns_df))
	end
	return select!(cns_df,Not(:cnsExpr))
end

# XXX adjusts resolution of var_df according to information in first two tuples
function checkTechReso!(tRes_tup::Tuple{Int,Int},cBalRes_tup::Tuple{Int,Int},var_df::DataFrame,sets_dic::Dict{Symbol,Tree})
	# if dispatch regions for technology were disaggregated, replace the disaggregated with the ones relevant for the carrier
	for (idx,dim) in enumerate([:Ts_dis,:R_dis])
		if cBalRes_tup[idx] != tRes_tup[idx]
			set_sym = Symbol(split(string(dim),"_")[1])
			dim_dic = Dict(x => getAncestors(x,sets_dic[set_sym],:int,cBalRes_tup[idx])[end] for x in unique(var_df[!,dim]))
			var_df[!,dim] = map(x -> dim_dic[x],var_df[!,dim])
		end
	end
end

# </editor-fold>


text:
# import Base.IO

const END_OF_DATA_SECTION = -1
const JULIA_EXCEPTION_THROWN = -2
const TIMING_DATA = -3
const END_OF_STREAM = -4
const NULL = -5
const PAIR_TUPLE = -6
const ARRAY_VALUE = -7
const ARRAY_END = -8
const INTEGER = -9
const STRING_START = -100

readint(io::IO) = ntoh(read(io, Int32))

"""
Read data object from a ioet. Returns code and byte array:
 * if code is negative, it's considered as a special command code
 * if code is positive, it's considered as array length
"""
function readobj(io::IO)
    len = readint(io)
    if (len > 0)
        (len , deserialized(read(io, len)))
    elseif (len == PAIR_TUPLE)
        res = (readobj(io)[2], readobj(io)[2])
        (len , res)
    elseif (len == ARRAY_VALUE)
        arr = Any[]
        while len == ARRAY_VALUE
            push!(arr, readobj(io)[2])
            len = readint(io)
        end
        (len, arr)
    elseif (len == ARRAY_END)
        (len, Any[])
    elseif (len < STRING_START)
        (len , String(read(io, -len + STRING_START)))
    elseif (len == STRING_START)
        (len, "")
    elseif (len == INTEGER)
        (len, ntoh(read(io, Int64)))
    else
        (len, [])
    end
end

writeint(io::IO, x::Int) = write(io, hton(Int32(x)))

"""Write object to stream"""
function writeobj(io::IO, obj::Any)
    sobj = serialized(obj)
    writeint(io, length(sobj))
    write(io, sobj)
end

function writeobj(io::IO, obj::Tuple{Any,Any})
    writeint(io, PAIR_TUPLE)
    writeobj(io, obj[1])
    writeobj(io, obj[2])
end

function writeobj(io::IO, str::AbstractString)
    utf8 = Vector(transcode(UInt8, str))
    writeint(io, STRING_START - length(utf8))
    write(io, utf8)
end

function writeobj(io::IO, x::Integer)
    writeint(io, INTEGER)
    write(io, hton(Int64(x)))
end

function load_stream(io::IO)
    ch = Channel{Any}(1024)
    code, _next = readobj(io)
    @async begin
        while code != END_OF_DATA_SECTION
            put!(ch, _next)
            code, _next = readobj(io)
        end
        close(ch)
    end
    return ch
end

function dump_stream(io::IO, it)
    for v in it
        writeobj(io, v)
    end
end


text:
using KernelMethods.KMap: fftraversal, size_criterion, dnet
import KernelMethods.KMap: centroid!
using SimilaritySearch
using NearNeighborGraph
using TextModel

export codebook, FFTraversal, ApproxKDCentroids

abstract type Clustering end

struct FFTraversal <: Clustering
    numcenters::Int
    k::Int
end

# Approximate k distance centroids
struct ApproxKDCentroids <: Clustering
    numcenters::Int
    k::Int
    maxiter::Int
    tol::Float64
    recall::Float64
end

ApproxKDCentroids(numcenters, k; maxiter=10, tol=0.001, recall=0.9) = ApproxKDCentroids(numcenters, k, maxiter, tol, recall)

function create_distance(::Type{VBOW})
    #cosine_distance
    angle_distance
end

function create_distance(::Type{DenseCosine{Float32}})
    #CosineDistance()
    AngleDistance()
end

function create_distance(::Type{Vector{F}}) where {F <: AbstractFloat}
    L2Distance()
end

"""
Creates an nearest neighbor index to speedup `codebook`
"""
function create_index(db, recall)
    T = typeof(db[1])
    # return Sequential(db, create_distance(T))
    index = LocalSearchIndex(T, create_distance(T), search=BeamSearch(), neighborhood=LogSatNeighborhood(), recall=recall)
    NearNeighborGraph.fit!(index, db)
    index
end

"""
Creates a sampling of the dataset using the farthest first algorithm
(an approximation of the kcenter problem)

"""
function codebook(algo::FFTraversal, X::AbstractVector{T}) where {T <: Union{DenseCosine{Float32},VBOW}}
    n = length(X)
    rlist = Int[]
    # dmax = 0.0
    function callback(pivID, _dmax)
        push!(rlist, pivID)
        # dmax = _dmax
    end

    xcodes = KnnResult[KnnResult(algo.k) for i in 1:n]
    function callbackdist(pivID, objID, d)
        push!(xcodes[objID], pivID, d)
    end

    fftraversal(callback, X, CosineDistance(), size_criterion(algo.numcenters), callbackdist)
    sort!(rlist)

    centers = X[rlist]
    H = Dict(r => i for (i, r) in enumerate(rlist))
    codes = [[H[p.objID] for p in xcodes[i]] for i in 1:n]
    centers, codes, rlist
end

function associate_centroids_and_score(algo, C, X, codes, distances)
    index = create_index(C, algo.recall)

    Threads.@threads for objID in 1:length(X)
        res = search(index, X[objID], KnnResult(algo.k))
        codes[objID] = [p.objID for p in res]
        distances[objID] = last(res).dist
    end

    mean(distances)
end

const _empty_int_array = Int[]
"""
Computes a clustering based on random sampling and an approximate nearest neighbor index, the procedure is iterative like the kmeans algorithm but instead of use the nearest neighbor `codebook` uses `k` nearest neighbors to converge and an approximate nn algorithm.

"""
function codebook(algo::ApproxKDCentroids, X::AbstractVector{T}) where {T <: Union{DenseCosine{Float32},VBOW}}
    n = length(X)
    rlist = rand(1:n, algo.numcenters) |> Set |> collect   # initial selection
    codes = [_empty_int_array for i in 1:n]
    C = X[rlist]
    iter = 0
    distances = zeros(Float64, n)
    scores = [typemax(Float64), associate_centroids_and_score(algo, C, X, codes, distances)]
    @assert !isnan(scores[end]) "ERROR invalid score $scores"
    while iter < algo.maxiter && abs(scores[end-1] - scores[end]) > algo.tol
        iter += 1
        info("*** starting iteration: $iter; scores: $scores ***")
        invindex = [Int[] for i in 1:algo.numcenters]
        for (objID, plist) in enumerate(codes)
            for refID in plist
                push!(invindex[refID], objID)
            end
        end

        info("*** computing centroids ***")
        Threads.@threads for i in 1:length(invindex)
            plist = invindex[i]
            # C[i] can be empty because we could be using approximate search
            if length(plist) > 0
                C[i] = centroid!(X[plist])
            end
        end

        info("*** computing $(algo.k) nearest references ***")
        s = associate_centroids_and_score(algo, C, X, codes, distances)

        push!(scores, s)
        @assert !isnan(scores[end]) "ERROR invalid score $scores"
        info("*** new score with $(algo.k) references: $scores ***")
    end

    info("*** finished computation of $(algo.k) references, scores: $scores ***")
    C, codes, rlist
end


text:
"""
    SARSOPSimulator

simulate a sarsop policy from the xml files
* `sim_len`  number of steps to use in simulation (default 10)
* `sim_num` number of simulations
* `fast` use fast (but very picky) alternate parser for .pomdp files
* `srand` set the random seed for the simulation
* `output_file` where to write the results
* `policy_filename` policy file location
* `pomdp_filename` pomdpx file location
"""
@with_kw struct SARSOPSimulator
    sim_len::Int = 10
    sim_num::Int = 10
    fast::Bool = false
    srand::Union{Int64, Nothing} = nothing
    output_file::AbstractString = ""
    policy_filename::AbstractString = "policy.out"
    pomdp_filename::AbstractString = "model.pomdpx"
end

function get_simulator_options(sim::SARSOPSimulator)
    options = String[]
    push!(options, "--simLen")
    push!(options, string(sim.sim_len))
    push!(options, "--simNum")
    push!(options, string(sim.sim_num))
    if sim.fast
        push!(options, "--fast")
    end
    if sim.srand != nothing
        push!(options, "--srand")
        push!(options, string(sim.srand))
    end
    push!(options, "--output-file")
    push!(options, sim.output_file)
    return options
end

"""
    POMDPs.simulate(sim::SARSOPSimulator)

Simulate a SARSOP policy according to the config specified by `sim::SARSOPSimulator`
"""
function POMDPs.simulate(sim::SARSOPSimulator)
    options_list = get_simulator_options(sim)
    pomdpsim() do pomdpsim_path
        run(`$pomdpsim_path $(sim.pomdp_filename) --policy-file $(sim.policy_filename) $options_list`)
    end
end

"""
    SARSOPEvaluator

simulate a sarsop policy from the xml files
* `sim_len`  number of steps to use in simulation (default 10)
* `sim_num` number of simulations
* `fast` use fast (but very picky) alternate parser for .pomdp files
* `srand` set the random seed for the simulation
* `memory` [MB] No memory limit by default. If memory usage exceeds the specified value, the evaluator will switch back
to a more conservative (and slow) methods.
* `output_file` where to write the results
* `policy_filename` policy file location
* `pomdp_filename` pomdpx file location
"""
@with_kw struct SARSOPEvaluator
    sim_len::Int = 10
    sim_num::Int = 10
    fast::Bool = false
    srand::Union{Int64, Nothing} = nothing
    memory::Union{Float64, Nothing} = nothing
    output_file::AbstractString = ""
    policy_filename::AbstractString = "policy.out"
    pomdp_filename::AbstractString = "model.pomdpx"
end

function get_evaluator_options(ev::SARSOPEvaluator)
    options = String[]
    push!(options, "--simLen")
    push!(options, string(ev.sim_len))
    push!(options, "--simNum")
    push!(options, string(ev.sim_num))
    if ev.fast
        push!(options, "--fast")
    end
    if ev.srand != nothing
        push!(options, "--srand")
        push!(options, string(ev.srand))
    end
    if ev.memory != nothing
        push!(options, "--memory")
        push!(options, string(ev.memory))
    end
    push!(options, "--output-file")
    push!(options, ev.output_file)
    return options
end

"""
    evaluate(ev::SARSOPEvaluator)

simulate a SARSOP policy according to the configuration specified by `ev::SARSOPEvaluator`
"""
function evaluate(ev::SARSOPEvaluator)
    options_list = get_evaluator_options(ev)
    pomdpeval() do pomdpeval_path
        run(`$pomdpeval_path $(ev.pomdp_filename) --policy-file $(ev.policy_filename) $options_list`)
    end
end


text:
# __precompile__(false)
# module GuSTO

using BulletCollision

using Plots
using RigidBodySim, RigidBodyDynamics
using MeshCat, MeshCatMechanisms
using CoordinateTransformations
using Interact, Reactive
using MAT, FileIO, MeshIO
using MechanismGeometries
using ForwardDiff
import GeometryTypes: HyperRectangle, HyperSphere, HomogenousMesh
import ColorTypes: RGB, RGBA

using StaticArrays, DataStructures
using JuMP, Convex
using Ipopt, Mosek, SCS
using MosekTools
using Gurobi
using MathProgBase, MathOptInterface
using NLsolve, DifferentialEquations

using MAT
using GeometryTypes
using FillArrays
using LinearAlgebra

const MOI = MathOptInterface

include("types.jl")
include("goals.jl")

include("traj_opt.jl")
include("scp.jl")
include("shooting.jl")

include("utils/quat_functions.jl")

include("robot.jl")
include("dynamics.jl")
include("environment.jl")

# end

max_stars_repo_name: StanfordASL/GuSTO

max_issues_repo_name: StanfordASL/GuSTO

max_forks_repo_name: StanfordASL/GuSTO


text:
module EMLvalve
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("valve/valve_simplified.jl")
	include("valve/valve_flow.jl")
	include("valve/valve.jl")
end