text:
using Colors
import Base: delete!

export Atom
export Bond
export Ensamble

include("Constants.jl")

struct Atom
    center::Point3f0
    symbol::String
    plot::AbstractPlot
end

struct Bond
    atom1::Atom
    atom2::Atom
    plot::Arrows
    Bond(A1, A2, p) = A1 != A2 ? new(A1, A2, p) : error("Cannot create a bound within a single atom")
end

mutable struct Ensamble
    fig::Figure
    scene::Scene
    mode::AbstractMode
    atoms::Vector{Atom}
    bonds::Vector{Bond}
    selected::Vector{Int}
end

function CreateAtom!(E::Ensamble, S::String, x::T, y::T, z::T) where T <: Number
    r = get_atom_radius(S)
    center = Point3f0(x,y,z)
    sphere = Sphere(center, r)
    plot = mesh!(E.scene, sphere, show_axis=false, color=atom_color[S])
    A = Atom(center, S, plot)
    push!(E.atoms, A)
    return A
end

function CreateBond!(E::Ensamble, i::Int, j::Int)
    A1, A2 = E.atoms[i], E.atoms[j]
    xyz = [A1.center]
    uvw = [A2.center - A1.center]
    plot = arrows!(E.scene, xyz, uvw, linewidth=0.1, arrowsize=0.0, linecolor=stegeman)
    B = Bond(E.atoms[i], E.atoms[j], plot)
    push!(E.bonds, B)
    return B
end

function CreateBond!(E::Ensamble)
    A1, A2 = E.atoms[end-1], E.atoms[end]
    xyz = [A1.center]
    uvw = [A2.center - A1.center]
    plot = arrows!(E.scene, xyz, uvw, linewidth=0.1, arrowsize=0.0, linecolor=stegeman)
    B = Bond(E.atoms[end-1], E.atoms[end], plot)
    push!(E.bonds, B)
    return B
end

function Ensamble()
    fig = Figure()
    s = fig.scene
    m = SelectionMode()
    atoms=Vector{Atom}()
    bonds=Vector{Bond}()
    selected=Vector{Int}()
    Ensamble(fig, s, m, atoms, bonds, selected)
end

function Ensamble(filename::String)
    E = Ensamble()
    read_xyz_to!(E, filename)
    find_bonds!(E)
    return E
end

function read_xyz_to!(E::Ensamble, filename)
    for l = eachline(filename)
        m = match(XYZ_REGEX, l)
        if m !== nothing
            S = String(m.captures[1])
            xyz = m.captures[2:end]
            @assert length(xyz) == 3
            xyz = parse.(Float32, xyz)
            CreateAtom!(E, S, xyz...)
        end
    end
end

function xyz_to!(E::Ensamble, xyz::String)
    for l = split(xyz, "\n")
        m = match(XYZ_REGEX, l)
        if m !== nothing
            S = String(m.captures[1])
            xyz = m.captures[2:end]
            @assert length(xyz) == 3
            xyz = parse.(Float32, xyz)
            CreateAtom!(E, S, xyz...)
        end
    end
end

function find_bonds!(E::Ensamble)

    for i in eachindex(E.atoms)
        r1 = E.atoms[i].center
        for j in (i+1):length(E.atoms)
            r2 = E.atoms[j].center
            d = √sum((r1 .- r2).^2)
            if d < 1.55
                CreateBond!(E, i, j)
            end
        end
    end
end

function delete!(E::Ensamble, A::Atom)

    # Remove atom from ensamble
    filter!(a -> a != A, E.atoms)

    # Delete any bonds involving the atom
    let to_be_deleted = Vector{Bond}()
        # Mark atoms for deletion
        for B in E.bonds
             if A in (B.atom1, B.atom2)
                 push!(to_be_deleted, B)
             end
        end
        # Delete them
        for B in to_be_deleted
            delete!(E, B)
        end
    end

    scene = E.scene
    plot = A.plot
    # Copied from Makie.jl source
    len = length(scene.plots)
    filter!(x -> x !== plot, scene.plots)
    if length(scene.plots) == len
        error("$(typeof(plot)) not in scene!")
    end
    for screen in scene.current_screens
        delete!(screen, scene, plot)
    end
end

function delete!(E::Ensamble, B::Bond)
    # Remove Bond from ensamble
    filter!(b -> b != B, E.bonds)

    scene = E.scene
    plot = B.plot

    # Copied from Makie.jl source
    len = length(scene.plots)
    filter!(x -> x !== plot, scene.plots)
    if length(scene.plots) == len
        error("$(typeof(plot)) not in scene!")
    end
    for screen in scene.current_screens
        delete!(screen, scene, plot)
    end
end

function isselected(E::Ensamble, A::Atom)
    for i = eachindex(E.atoms)
        if A === E.atoms[i]
            if i in E.selected
                return true
            else
                return false
            end
        end
    end
    @error "Atom not found"
end

function isselected(E::Ensamble, B::Bond)
    for i = eachindex(E.bonds)
        if B === E.bonds[i]
            if -i in E.selected
                return true
            else
                return false
            end
        end
    end
    @error "Bond not found"
end

function select!(E::Ensamble, A::Atom)
    if !isselected(E, A)
        for i = eachindex(E.atoms)
            if A === E.atoms[i]
                push!(E.selected, i)
                A.plot.attributes[:color] = (:lightgoldenrod1, 0.90)
                A.plot.attributes[:shininess] = 70.0f0
                A.plot.attributes[:diffuse] = Float32[0.6, 0.6, 0.6]
                println(E.selected)
            end
        end
    else
        for i = eachindex(E.atoms)
            if A === E.atoms[i]
                filter!(x -> x !=i, E.selected)
                A.plot.attributes[:color] = atom_color[A.symbol]
                A.plot.attributes[:diffuse] = Float32[0.4, 0.4, 0.4]
                println(E.selected)
            end
        end
    end
end

function select!(E::Ensamble, B::Bond)
    if !isselected(E, B)
        for i = eachindex(E.bonds)
            if B === E.bonds[i]
                push!(E.selected, -i)
                B.plot.attributes[:color] = (:lightgoldenrod1, 0.90)
                B.plot.attributes[:diffuse] = Float32[0.6, 0.6, 0.6]
            end
        end
    else
        for i = eachindex(E.bonds)
            if B === E.bonds[i]
                filter!(x -> x != -i, E.selected)
                B.plot.attributes[:color] = atom_color[A.symbol]
                B.plot.attributes[:diffuse] = Float32[0.4, 0.4, 0.4]
            end
        end
    end
end

function clear_selection!(E::Ensamble)
    selected_structs = deepcopy(E.selected)
    for i in selected_structs
        if i > 0
            select!(E, E.atoms[i])
        else
            select!(E, E.bonds[-i])
        end
    end
end


text:
module Diagnostics

    using ..Tokens
    using ..Tokens: SourceRange, SourceNode, SourceLocToken, SourceExpr, √
    using ..LineNumbers: SourceFile, compute_line

    abstract AbstractDiagnostic

    immutable Message
        severity::Symbol
        location::SourceRange
        text::AbstractString
    end

    immutable Diagnostic <: AbstractDiagnostic
        elements::Vector{Message}
    end

    immutable Incomplete <: AbstractDiagnostic
        tag::Symbol
        d::Diagnostic
    end

    function after(loc)
        loc = normalize_loc(loc)
        loc == SourceRange() && return SourceRange()
        SourceRange(loc.offset+loc.length,1,loc.file)
    end

    function before(loc)
        loc = normalize_loc(loc)
        loc == SourceRange() && return SourceRange()
        SourceRange(loc.offset-1,1,loc.file)
    end

    function normalize_loc(loc)
        isa(loc, SourceNode) && (loc = loc.loc)
        isa(loc, Union{SourceExpr, SourceLocToken}) && (loc = √loc)
        isa(loc, Token) && (loc = SourceRange())
        loc == nothing && (loc = SourceRange())
        loc
    end

    function diag(loc, message, severity = :error)
        loc = normalize_loc(loc)
        Diagnostic([Message(severity, loc, message)])
    end

    function diag(D::Diagnostic, loc, message, severity = :note)
        loc = normalize_loc(loc)
        loc !== nothing && push!(D.elements,Message(severity, loc, message))
    end
    diag(D::Incomplete, loc, message, severity = :note) = diag(D.d, loc, message, severity)

    function display_diagnostic(io::IO, code, diag; filename = "none")
        file = SourceFile(code)
        for message in diag.elements
            if message.location == SourceRange()
                # Don't show notes if they don't have a location
                if message.severity == :note || message.severity == :fixit
                    continue
                end
                print_with_color(message.severity == :error ? :red : :magenta, io, string(message.severity))
                println(io, ": ", message.text)
                continue
            end
            offset = message.location.offset
            line = compute_line(file, offset)
            str  = String(file[line])
            col  = sum(charwidth, str[1:(offset - file.offsets[line] + 1)])
            if message.severity == :fixit
                print(io, " "^(col-1))
                print_with_color(:green, io, message.text)
                println(io)
            else
                print(io, "$filename:$line:$col " )
                print_with_color(message.severity == :error ? :red : :magenta, io, string(message.severity))
                println(io, ": ", message.text)
                println(io, rstrip(str))
                print(io, " "^(col-1))
                print_with_color(:green, io, string('^',"~"^(max(0,message.location.length-1))))
                println(io)
            end
        end
    end
end

ext: jl


text:
"""
    PreprocessedDataCausalityTest

A causality test where the data is subjected to some preprocessing before
applying the test.
"""
abstract type PreprocessedDataCausalityTest end

export PreprocessedDataCausalityTest


text:
@testset "Renewable data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type RenewableDispatch, consider changing the device models"
    model = DeviceModel(RenewableDispatch, RenewableFullDispatch)
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Renewable, model)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Renewable, model)
end

@testset "Renewable DCPLossLess FullDispatch" begin
    model = DeviceModel(RenewableDispatch, RenewableFullDispatch)
    c_sys5_re = build_system("c_sys5_re")

    #5 Bus testing case
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_re)
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, false, 72, 0, 72, 0, 0, false)

    psi_checkobjfun_test(op_problem, GAEVF)

    # Using Parameters Testing
    op_problem =
        OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_re; use_parameters = true)
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, true, 72, 0, 72, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast - No Parameters Testing
    op_problem = OperationsProblem(
        TestOpProblem,
        DCPPowerModel,
        c_sys5_re;
        use_forecast_data = false,
    )
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, false, 3, 0, 3, 3, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast Testing
    op_problem = OperationsProblem(
        TestOpProblem,
        DCPPowerModel,
        c_sys5_re;
        use_parameters = true,
        use_forecast_data = false,
    )
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, true, 3, 0, 3, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Renewable ACPPower Full Dispatch" begin
    model = DeviceModel(RenewableDispatch, RenewableFullDispatch)
    c_sys5_re = build_system("c_sys5_re")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_re; use_parameters = p)
        construct_device!(op_problem, :Renewable, model)
        if p
            moi_tests(op_problem, p, 144, 0, 144, 72, 0, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        else
            moi_tests(op_problem, p, 144, 0, 144, 72, 0, false)

            psi_checkobjfun_test(op_problem, GAEVF)
        end
    end
    # No Forecast Test
    for p in [true, false]
        op_problem = OperationsProblem(
            TestOpProblem,
            ACPPowerModel,
            c_sys5_re;
            use_forecast_data = false,
            use_parameters = p,
        )
        construct_device!(op_problem, :Renewable, model)
        if p
            moi_tests(op_problem, p, 6, 0, 6, 3, 0, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        else
            moi_tests(op_problem, p, 6, 0, 6, 6, 0, false)

            psi_checkobjfun_test(op_problem, GAEVF)
        end
    end
end

@testset "Renewable DCPLossLess ConstantPowerFactor" begin
    model = DeviceModel(RenewableDispatch, RenewableConstantPowerFactor)
    c_sys5_re = build_system("c_sys5_re")

    #5 Bus testing case
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_re)
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, false, 72, 0, 72, 0, 0, false)

    psi_checkobjfun_test(op_problem, GAEVF)

    # Using Parameters Testing
    op_problem =
        OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_re; use_parameters = true)
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, true, 72, 0, 72, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast - No Parameters Testing
    op_problem = OperationsProblem(
        TestOpProblem,
        DCPPowerModel,
        c_sys5_re;
        use_forecast_data = false,
    )
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, false, 3, 0, 3, 3, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast Testing
    op_problem = OperationsProblem(
        TestOpProblem,
        DCPPowerModel,
        c_sys5_re;
        use_parameters = true,
        use_forecast_data = false,
    )
    construct_device!(op_problem, :Renewable, model)
    moi_tests(op_problem, true, 3, 0, 3, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Renewable ACPPower ConstantPowerFactor" begin
    model = DeviceModel(RenewableDispatch, RenewableConstantPowerFactor)
    c_sys5_re = build_system("c_sys5_re")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_re; use_parameters = p)
        construct_device!(op_problem, :Renewable, model)
        if p
            moi_tests(op_problem, p, 144, 0, 72, 0, 72, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        else
            moi_tests(op_problem, p, 144, 0, 72, 0, 72, false)

            psi_checkobjfun_test(op_problem, GAEVF)
        end
    end
    # No Forecast Test
    for p in [true, false]
        op_problem = OperationsProblem(
            TestOpProblem,
            ACPPowerModel,
            c_sys5_re;
            use_forecast_data = false,
            use_parameters = p,
        )
        construct_device!(op_problem, :Renewable, model)
        if p
            moi_tests(op_problem, p, 6, 0, 3, 0, 3, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        else
            moi_tests(op_problem, p, 6, 0, 3, 3, 3, false)

            psi_checkobjfun_test(op_problem, GAEVF)
        end
    end
end

@testset "Renewable DCPLossLess FixedOutput" begin
    model = DeviceModel(RenewableDispatch, FixedOutput)
    c_sys5_re = build_system("c_sys5_re")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_re; use_parameters = p)
        construct_device!(op_problem, :Renewable, model)
        if p
            moi_tests(op_problem, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        else
            moi_tests(op_problem, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        end
    end
end

@testset "Renewable ACPPowerModel FixedOutput" begin
    model = DeviceModel(RenewableDispatch, FixedOutput)
    c_sys5_re = build_system("c_sys5_re")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_re; use_parameters = p)
        construct_device!(op_problem, :Renewable, model)
        if p
            moi_tests(op_problem, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        else
            moi_tests(op_problem, p, 0, 0, 0, 0, 0, false)
            psi_checkobjfun_test(op_problem, GAEVF)
        end
    end
end


text:
struct WSolver
    solver::WSymbol
end

macro WS_str(s::AbstractString)
    WSolver(WSymbol(s))
end

function wjaggedlist(x::AbstractArray{WElem,N}) where {N}
    if N == 1
        WExpr(W"List", map(y -> y.val, x))
    else
        s = size(x)
        x = reshape(x, s[1], :)
        WExpr(W"List", wjaggedlist.(reshape(x[i, :], s[2:N]) for i in 1:s[1]))
    end
end

function parsewlist(x::WExpr, reccall::Bool = false)
    iswlist(x) || error("")
    if all(y -> iswlist(y), x.args)
        x = parsewlist.(x.args, true)
        x = cat(x...; dims = ndims(first(x)) + 1)
        reccall ? x : permutedims(x, ndims(x):-1:1)
    else
        WElem.(x.args)
    end
end

function (s::WSolver)(a::AbstractArray{WElem}, rissingle::Bool = false; kargs...)
    r = weval((s.solver)(wjaggedlist(a)); kargs...)
    if rissingle
        iswlist(r) ? parsewlist(r) : WElem(r)
    else
        iswlist(r) || error("")
        map(x -> iswlist(x) ? parsewlist(x) : WElem(x), r.args)
    end
end


text:
module AgentsX

using Agents
using Mustache

include("layers.jl")
include("ABM.jl")
include("codegen.jl")

export ParamLayer
export AbstractLayer

export agent_step!
export AxBM
export create_agent_template

end


text:
include("../../utils.jl")
using Random

mutable struct TreeNode
   key::Int
   left::TreeNode
   right::TreeNode
   TreeNode() = new()
   TreeNode(x) = new(x)
   TreeNode(x,y,z) = new(x,y,z)
end

function insert(key, n::TreeNode)
   if key < n.key
      if !isdefined(n,:left)
         n.left = TreeNode(key)
      else
         insert(key, n.left)
      end
   elseif key > n.key
      if !isdefined(n,:right)
         n.right = TreeNode(key)
      else
         insert(key, n.right)
      end
   end
end

function sumTree(n::TreeNode)
   sum = n.key
   if isdefined(n,:left)
      sum += sumTree(n.left)
   end
   if isdefined(n,:right)
       sum += sumTree(n.right)
   end

   return sum
end

# tree_size is the number of elements in mb
function tree(n=4)
    n *= 1024^2
    rng = Xoshiro(12345)
    temp = rand(rng, Int, n)
    root::TreeNode = TreeNode(temp[1])
    for i = 2:n
       insert(temp[i], root)
    end
    return sumTree(root)
end

@gctime tree()

max_stars_repo_name: d-netto/GCBenchmarks

max_issues_repo_name: d-netto/GCBenchmarks

max_forks_repo_name: d-netto/GCBenchmarks


text:
using DataStructures
using Logging

import Base.show
# package code goes here

Logging.configure(level=DEBUG)
typealias TT_TYPE Array{UInt,1}
typealias IDX_TYPE UInt
typealias TV_TYPE Array{Real,1}

function __init__()
	println("loading MJPlayGround.jl")
	Logging.configure(level=DEBUG)
end

function init()
	global vvals = TV_TYPE()
	global pvals = TV_TYPE()
	global vidx = 1
	global pidx = 1
	global tt = TT_TYPE()
end

#export
#	Var, Par, @test, @var, @test1

#the AD types below
const TYPE_V = 1	#variable node
const TYPE_P = 2	#param node
const TYPE_OU = 3	#unary op
const TYPE_OB = 4  	#binary op
abstract Placeholder
immutable AD_V <: Placeholder
	idx::IDX_TYPE  #index on tape
	t::Uint    #type code
end
immutable AD_P <: Placeholder
	idx::IDX_TYPE  #index on tape
	t::Uint    #type code
end
immutable AD_O <: Placeholder
	idx::IDX_TYPE  #index on tape
    t::Uint    #type code
end

function AD_V(tt::TT_TYPE,val=NaN)
	push!(tt,vidx)
	push!(tt,TYPE_V)
	push!(vvals,val)
	global vidx +=1
	this = AD_V(length(tt),TYPE_V)
	return this
end

function AD_P(tt::TT_TYPE,val=NaN)
	push!(tt,pidx)
	push!(tt,TYPE_P)
	push!(pvals,val)
	global pidx +=1
	this = AD_P(length(tt),TYPE_P)
	return this
end

function AD_O(tt::TT_TYPE,oc, lidx::UInt)
	push!(tt,lidx)
	push!(tt,oc)
	push!(tt,TYPE_OU)
	this = AD_O(length(tt),TYPE_OU)
	return this
end

function AD_O(tt::TT_TYPE,oc, lidx::UInt, ridx::UInt)
	push!(tt,ridx)
	push!(tt,lidx)
	push!(tt,oc)
	push!(tt,TYPE_OB)
	this = AD_O(length(tt),TYPE_OB)
	return this
end

function Base.show(io::IO,m::AD_V)
	print("AD_V[",m.idx,"]")
end

function Base.show(io::IO,m::AD_P)
	print("AD_P[",m.idx,"]")
end

function Base.show(io::IO,m::AD_O)
	# print("(",m.l)
	# print(string(getOpSym(m)))
	# if(!isnull(m.r))
	# 	print(m.r.value)
	# end
	# print(")")
	print("AD_O[",m.idx,"]")
end

#Operator overloading function for AD types
const OP = (:+,:-,:*,:/,:sin,:cos)

function getOpSym(oc)
	return OP[oc]
end

const B_OP_START = 1
const B_OP_END = 4
for oc = B_OP_START:1:B_OP_END
	o = OP[oc]
	debug("setup operator ",o)
	eval(quote
			($o)(l::Placeholder,r::Placeholder) =
			begin
				debug("at",OP[$(oc)])
				this = AD_O(tt,$(oc),l.idx,r.idx)
				return this
			end
		end)
end

const U_OP_START = 5
const U_OP_END = 6
for oc = U_OP_START:1:U_OP_END
	o = OP[oc]
	debug("setup operator ",o)
	eval(quote
			($o)(l::Placeholder) =
			begin
				debug("at",$(oc))
				this = AD_O(tt,$(oc),l.idx)
				return this
			end
		end)
end

function evaluate(tt::TT_TYPE,idx::IDX_TYPE)
	println("enter - ",idx)
	ret = NaN
	if(idx != 0)
		ntype = tt[idx] #type node
		idx -= 1
		if(ntype == TYPE_P)
			ret = pvals[tt[idx]]
			idx -=1
		elseif(ntype == TYPE_V)
			ret = vvals[tt[idx]]
			idx -= 1
		elseif(ntype == TYPE_OB)
			oc = tt[idx]
			assert(B_OP_START<= oc <= B_OP_END)
			idx -= 1
			lidx = tt[idx]
			idx -= 1
			ridx = tt[idx]
			idx -= 1
			println("before left - ",lidx)
			(lval,idx) = evaluate(tt,lidx)
			println("before right - ",idx)
			(rval,idx) = evaluate(tt,ridx)
			ret = evaluate(oc,lval,rval)
		else
			oc = tt[idx]
			assert(U_OP_START <= oc <= U_OP_END)
			idx -= 1
			lidx = tt[idx]
			idx -= 1
			println("unary before  - ",idx)
			(lval,idx) = evaluate(tt,lidx)
			ret = evaluate(oc,lval)
		end
		println("now - ",idx)
	end
	assert(!isnan(ret))
	println("exit - ",idx)
	return ret,idx
end

function evaluate(oc::UInt,lval::Real,rval::Real)
	assert(B_OP_START<= oc <= B_OP_END)
	ex = Expr(:call,OP[oc],lval,rval)
	println("evaluate:",ex)
	return eval(ex)
end

function evaluate(oc::UInt,lval::Real)
	assert(U_OP_START<= oc <= U_OP_END)
	ex = Expr(:call,OP[oc],lval)
	println("evaluate:",ex)
	return eval(ex)
end

## Interface method
function feval(tt::TT_TYPE)
	(val,idx) = evaluate(tt,convert(UInt,length(tt)))
	assert(idx == 0)
	return val
end

#macro rewrite(expr)
#	print(expr)
#	debug("rewrite ",expr)
#	return :(4+3)
#end

macro var(x)
	quote
	debug("var - ",$(esc(x)i))
	#v = Var($global_idxi)
	#$(esc(global_idx)) += 1
	global global_idx += 1
	#push!(s,v)
	end
end

macro test(args...)
	println("---",length(args))
	arg1 = args[1]
	arg2 = args[2]
	code = quote
		local a = $(esc(m))
		#local b = $expr;
		local c = $(esc(arg1))
		dump($(esc(arg1)))
	end
	return code
end

ex = :(3+3)

macro test1(ex)
	code =:()
	code = :($code;ex=3
		#dump($(esc(ex)))
		#println($(esc(ex)))
		#a = $(esc(ex))
		#dump($ex)
		#println(eval($(ex)))
	 	#$(ex), $(esc(ex))
	)
	return code
end

macro test2(ex...)
	quote
		x = 1
		$(esc(ex...))
		nothing
	end
end

macro p(n)
    if typeof(n) == Expr
       println(n.args)
    end
    eval(n)
end

macro f(x)
       quote
           s = 4
           (s, $(esc(s)))
       end
end

macro dotimes(n,body)
	quote
		for i=1:$(esc(n))
			$(esc(body))
		end
	end
end


text:
function JPLEphemeridesLoader(arg0::JString, arg1::JPLEphemeridesLoader_EphemerisType)
    return JPLEphemeridesLoader((JString, JPLEphemeridesLoader_EphemerisType), arg0, arg1)
end

function JPLEphemeridesLoader(arg0::JString, arg1::JPLEphemeridesLoader_EphemerisType, arg2::DataProvidersManager, arg3::TimeScales, arg4::Frame)
    return JPLEphemeridesLoader((JString, JPLEphemeridesLoader_EphemerisType, DataProvidersManager, TimeScales, Frame), arg0, arg1, arg2, arg3, arg4)
end

function get_loaded_astronomical_unit(obj::JPLEphemeridesLoader)
    return jcall(obj, "getLoadedAstronomicalUnit", jdouble, ())
end

function get_loaded_constant(obj::JPLEphemeridesLoader, arg0::Vector{JString})
    return jcall(obj, "getLoadedConstant", jdouble, (Vector{JString},), arg0)
end

function get_loaded_earth_moon_mass_ratio(obj::JPLEphemeridesLoader)
    return jcall(obj, "getLoadedEarthMoonMassRatio", jdouble, ())
end

function get_loaded_gravitational_coefficient(obj::JPLEphemeridesLoader, arg0::JPLEphemeridesLoader_EphemerisType)
    return jcall(obj, "getLoadedGravitationalCoefficient", jdouble, (JPLEphemeridesLoader_EphemerisType,), arg0)
end

function get_max_chunks_duration(obj::JPLEphemeridesLoader)
    return jcall(obj, "getMaxChunksDuration", jdouble, ())
end

function load_celestial_body(obj::JPLEphemeridesLoader, arg0::JString)
    return jcall(obj, "loadCelestialBody", CelestialBody, (JString,), arg0)
end


text:
export EliminateGraph, subgraph, refresh, check_validity
export nv0, nv, ne0, ne
export isconnected, unsafe_connect!, unsafe_disconnect!, find_cluster

"""
    EliminateGraph <: AbstractGraph
    EliminateGraph(tbl::AbstractMatrix) -> EliminateGraph

A graph type for algorithms that involve node elimination.
With this type, vertex elimination and recover do not allocate.
`tbl` in the constructor is a boolean table for connection.
"""
mutable struct EliminateGraph{T} <: LightGraphs.AbstractGraph{T}
    tbl::Matrix{Bool}
    vertices::Vector{T}
    ptr::Vector{Int}
    level::Int
    nv::Int
end

EliminateGraph(tbl::AbstractMatrix) = EliminateGraph(Matrix{Bool}(tbl))
function EliminateGraph(tbl::Matrix{Bool})
    N = size(tbl, 1)
    vertices = collect(1:N)
    ptr = zeros(Int,N)
    EliminateGraph(tbl, vertices, ptr, 0, N)
end
function EliminateGraph(nv::Int, pairs::Vector{Pair{Int,Int}})
    tbl = zeros(Bool, nv, nv)
    for (vi, vj) in pairs
        tbl[vi,vj] = true
    end
    return EliminateGraph(tbl .| tbl')
end

Base.copy(eg::EliminateGraph) = EliminateGraph(eg.tbl, eg.vertices |> copy, eg.ptr|>copy, eg.level, eg.nv)
Base.:(==)(eg1::EliminateGraph, eg2::EliminateGraph) = eg1.tbl == eg2.tbl && Set(vertices(eg1)) == Set(vertices(eg2))

"""initial size of a `EliminateGraph`."""
nv0(eg::EliminateGraph) = size(eg.tbl, 1)
"""initial number of edges in a `EliminateGraph`."""
ne0(eg::EliminateGraph) = sum(eg.tbl) ÷ 2
"""current size of a `EliminateGraph`."""
nv(eg::EliminateGraph) = eg.nv
"""current number of edges in a `EliminateGraph`."""
function ne(eg::EliminateGraph)
    res = 0
    vs = vertices(eg)
    for j in 1:nv(eg)
        @inbounds vj = vs[j]
        @inbounds for vi in vertices(eg)[j+1:end]
            res += isconnected(eg, vi, vj)
        end
    end
    return res
end

"""
A eliminate graph is valid.
"""
function check_validity(eg::EliminateGraph)
    mes = ""
    eg.tbl' == eg.tbl || (mes *= "connection table not symmetric\n")
    !any(i->eg.tbl[i,i], 1:nv0(eg)) || (mes *= "diagonal part of table non-zero\n")
    Set(eg.vertices) == Set(1:nv0(eg)) || (mes *= "vertices list incorrect\n")
    0 <= eg.nv <= nv0(eg) || (mes *= "number of vertices not in range\n")
    0 <= eg.level <= nv0(eg) || (mes *= "level not in range\n")
    if eg.level > 0
        eg.ptr[eg.level] == nv0(eg)-eg.nv || (mes *= "ptr points to incorrect position\n")
        nv0(eg) >= eg.ptr[1] > 0 || (mes *= "ptr not in range\n")
        all(diff(eg.ptr[1:eg.level]) .> 0) || (mes *= "non-increasing ptr")
    end
    if mes == ""
        return true
    else
        println(mes)
        return false
    end
end

"""undo elimination for a `EliminateGraph`."""
refresh(eg::EliminateGraph) = EliminateGraph(eg.tbl)

"""
    isconnected(eg::EliminateGraph, vi::Int, vj::Int) -> Bool

Return true if `vi`, `vj` are connected in `eg`.
Note: This function does not check `vi`, `vj` out of bound error!
"""
isconnected(eg::EliminateGraph, vi::Int, vj::Int) = @inbounds eg.tbl[vi,vj]

getid(eg::EliminateGraph, vi::Int) = findfirst(==(vi), vertices(eg))

"""
    find_cluster(eg::EliminateGraph, vi::Int) -> Vector{Int}

Find the cluster connected to `vi` in `eg`.
"""
function find_cluster(eg::EliminateGraph, vi::Int)
    vset = zeros(Bool, nv0(eg))
    find_cluster!(eg, vi, vset)
    return findall(!iszero, vset)
end

function find_cluster!(eg::EliminateGraph, vi::Int, vset::Vector{Bool})
    vset[vi] && return

    vset[vi] = true
    for vj in neighbors(eg, vi)
        find_cluster!(eg, vj, vset)
    end
end

function isconnected(eg::EliminateGraph)
    eg.nv <=1 && return true
    length(find_cluster(eg, vertices(eg)[1])) == nv(eg)
end

subgraph(eg::EliminateGraph, vs) = EliminateGraph(eg.tbl[vs, vs])

"""
    unsafe_connect!(eg::EliminateGraph, vi::Int, vj::Int) -> EliminateGraph

connect two vertices.
"""
function unsafe_connect!(eg::EliminateGraph, vi::Int, vj::Int)
    @inbounds eg.tbl[vi, vj] = true
    @inbounds eg.tbl[vj, vi] = true
    eg
end

"""
    unsafe_disconnect!(eg::EliminateGraph, vi::Int, vj::Int) -> EliminateGraph

connect two vertices.
"""
function unsafe_disconnect!(eg::EliminateGraph, vi::Int, vj::Int)
    @inbounds eg.tbl[vi, vj] = false
    @inbounds eg.tbl[vj, vi] = false
    eg
end

function Base.show(io::IO, eg::EliminateGraph)
    N = nv0(eg)
    println(io, "EliminateGraph")
    vs = vertices(eg)
    for i=1:N
        for j=1:N
            print(io, "  ", (i in vs && j in vs) ? Int(eg.tbl[i,j]) : "⋅")
        end
        println(io)
    end
end

# convert to light graphs
LightGraphs.SimpleGraph(eg::EliminateGraph) = SimpleGraph(ne(eg), [[neighbors(eg, iv)...] for iv=1:nv(eg)])
function EliminateGraph(sg::SimpleGraph)
    N = nv(sg)
    tbl = zeros(Bool, N, N)
    for iv = 1:N
        @inbounds for jv in sg.fadjlist[iv]
            tbl[iv, jv] = true
        end
    end
    EliminateGraph(tbl)
end


text:
# # Example 2: Optimization of a State-to-State Transfer in a Lambda System in the RWA

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`state_to_state_rwa.ipynb`](@__NBVIEWER_ROOT_URL__/examples/state_to_state_rwa.ipynb)

# This example is illustrates the use of complex-valued control fields. This is
# accomplished by rewriting the Hamiltonian as the sum of two independent
# controls (real and imaginary parts). We consider a 3-level system in a
# $\Lambda$ configuration.

const σ̂_z = ComplexF64[1 0; 0 -1]
const σ̂_x = ComplexF64[0 1; 1  0]

"""Two-level-system Hamiltonian."""
function hamiltonian(Ω=1.0, E0=0.2)

    Ĥ₀ = -0.5 * Ω * σ̂_z
    Ĥ₁ = σ̂_x

    return [Ĥ₀, Ĥ₁]

end


text:
# Test the availability of drivers.
# See https://github.com/JuliaGeo/GDAL.jl/issues/65 for an overview.

# not a complete list, but these drivers are tested for availability
available_drivers = [
    # raster drivers
    "AAIGrid",
    "GPKG",
    "GTiff",
    "MEM",
    "PCRaster",
    # libcurl raster drivers
    "EEDAI",
    "PLMosaic",
    "RDA",
    "WCS",
    "WMS",
    "WMTS",
    # vector drivers
    "ARCGEN",
    "GeoJSON",
    "GML",
    "GPKG",
    "KML",
    "OpenFileGDB",
    "ESRI Shapefile",
    "TopoJSON",
    "VRT",
    "SQLite",
    # libcurl vector drivers
    "AmigoCloud",
    "Carto",
    "Cloudant",
    "CouchDB",
    "CSW",
    "EEDA",
    "ElasticSearch",
    "GFT",
    "NGW",
    "PLScenes",
    "WFS",
    "WFS3",
]

for drivername in available_drivers
    @test GDAL.gdalgetdriverbyname(drivername) != C_NULL
end

not_available_drivers = [
    "JP2OpenJPEG",  # up next, see #64
]

for drivername in not_available_drivers
    @test_broken GDAL.gdalgetdriverbyname(drivername) != C_NULL
end


text:
ENV["JULIA_PYTHONCALL_EXE"] = "/home/zhangyong/miniconda3/bin/python"
using Faiss
using PythonCall
using ProgressMeter
using NearestNeighbors
# np = pyimport("numpy")

function test()
    dir_1 = "/mnt/zy_data/data/longhu_1/sorted_2/"
    feats = np.load(joinpath(dir_1, "feats.npy"))
    println(typeof(feats), feats.shape)
    # feats = PyArray{Float32, 2, true, true}(feats)
    feats = pyconvert(Array{Float32, 2}, feats)  # py to julia Matrix
    # D, I = local_rank(vs_query, vs_gallery, k=10, metric="IP", gpus="")

    feat_dim = size(feats, 2)
    idx = Index(feat_dim; str="IDMap2,Flat", metric="L2", gpus="4")  # IDMap2. L2,IP
    Faiss.show(idx)
    k = 10
    @showprogress for i in range(1, 1000)
        vs_gallery = feats[100*i+1:100*(i+1),:]
        # println(typeof(feats), size(feats))
        vs_query = vs_gallery

        # D, I = add_search(idx, vs_query, vs_gallery; k=10, flag=true, metric="cos")
        # D, I = add_search_with_ids(idx, vs_query, vs_gallery; k=10)
        ids = collect(range(100*i+1, 100*(i+1))) .+ 100
        println(typeof(ids), size(ids))
        add_with_ids(idx, vs_gallery, ids)
        D, I = search(idx, vs_query, k)
        # println(typeof(D), size(D))
        # println(typeof(I), size(I))
        if i == 2
            println(D[1:2, 1:5])
            println(I[1:2, 1:5])
            break
        end
    end
end

function faiss_test_1()
    feats = rand(10^6, 128);
    top_k = 100
    vs_gallery = feats;
    vs_query = feats[1:10^4, :];

    D, I = local_rank(vs_query, vs_gallery, k=top_k, metric="IP", gpus="")

    # println(typeof(D), size(D))
    # println(D[1:5, :])
end

function nn_test()
    data = rand(128, 10^6)
    k = 100
    query = rand(128, 10^4)

    kdtree = KDTree(data)
    idxs, dists = knn(kdtree, query, k, true)
    # println(typeof(idxs), size(idxs))
end

# @time test()
@time faiss_test_1()
# @time nn_test()

#=
julia --project=/home/zhangyong/codes/Faiss.jl/Project.toml "/home/zhangyong/codes/Faiss.jl/test/test_1.jl"

cpu: bs=100,k=10
237.4seconds=4min (6.41 M allocations: 800.244 MiB, 0.13% gc time, 1.39% compilation time)
gpu: bs=100,k=10
25.010865 seconds (6.38 M allocations: 797.798 MiB, 1.01% gc time, 14.27% compilation time)

=#


text:
mutable struct CUHK_SYSU
    w # weight
end

function CUHK_SYSU(weight::Int, scale=0.01)
    w = weight
    return CUHK_SYSU(w)
end

function (l::CUHK_SYSU)(x)
    return l.w[:, x]
end

function load_queries(l::CUHK_SYSU)
    queries = []
    # TODO replicate from seqnet code
    return queries
end

function load_split_img_names(l::CUHK_SYSU)
    images = []
    # TODO replicate from seqnet code
    return images
end

function load_annotations(l::CUHK_SYSU)
    annotations = []
    # TODO replicate from seqnet code
    return annotations
end

max_stars_repo_name: ucatbas/seqnet-julia

max_issues_repo_name: ucatbas/seqnet-julia

max_forks_repo_name: ucatbas/seqnet-julia


text: Dict{String, Any}("contributors" => Any[Dict{String, Any}("value" => "Foo Bar <foo@example.com>", "type" => "string"), Dict{String, Any}("name" => Dict{String, Any}("value" => "Baz Qux", "type" => "string"), "url" => Dict{String, Any}("value" => "https://example.com/bazqux", "type" => "string"), "email" => Dict{String, Any}("value" => "bazqux@example.com", "type" => "string"))])


text:
abstract type AbstractUnconstrainedSolution <: AbstractSolution end
abstract type AbstractConstrainedSolution   <: AbstractSolution end
# multi-objective solution are constrained by default
abstract type AbstractMultiObjectiveSolution <: AbstractConstrainedSolution end

mutable struct xf_indiv <: AbstractUnconstrainedSolution # Single Objective
    x::Vector{Float64}
    f::Float64
end

mutable struct xfgh_indiv <: AbstractConstrainedSolution # Single Objective Constraied
    x::Vector{Float64}
    f::Float64
    g::Vector{Float64}
    h::Vector{Float64}
    sum_violations::Float64 # ∑ max(0,g) + ∑|h|
    is_feasible::Bool

end

function xfgh_indiv(
    x::Vector{Float64},
    f::Float64,
    g::Vector{Float64},
    h::Vector{Float64};
    sum_violations = 0.0,
    ε = 0.0
)
    if sum_violations <= 0.0
        sum_violations = violationsSum(g, h; ε=ε)
    end

    xfgh_indiv(x, f, g, h, sum_violations, sum_violations == 0.0)
end

mutable struct xFgh_indiv <: AbstractMultiObjectiveSolution# Single Objective Constraied
    x::Vector{Float64}
    f::Vector{Float64}
    g::Vector{Float64}
    h::Vector{Float64}
    rank::Int
    crowding::Float64
    sum_violations::Float64 # ∑ max(0,g) + ∑|h|
    is_feasible::Bool
end

function xFgh_indiv(
    x::Vector{Float64},
    f::Vector{Float64},
    g::Vector{Float64},
    h::Vector{Float64};
    rank = 0,
    crowding = 0.0,
    sum_violations = 0.0,
    ε = 0.0
)

    if sum_violations <= 0
        sum_violations = violationsSum(g, h;ε=ε)
    end
    xFgh_indiv(x, f, g, h, Int(rank), crowding, sum_violations, sum_violations == 0.0)
end

const Solution = Union{xf_indiv, xfgh_indiv, xFgh_indiv}
const Population = Array{Solution, 1}

############################################################
# Generate solutions depending on the objective function
# output
#############################################################

# single objective
"""
    Metaheuristics.create_child(x, fx)

Constructor for a solution depending on the result of `fx`.

## Example

```
julia> import Metaheuristics

julia> Metaheuristics.create_child(rand(3), 1.0)
| f(x) = 1
| solution.x = [0.2700437125780806, 0.5233263210622989, 0.12871108215859772]

julia> Metaheuristics.create_child(rand(3), (1.0, [2.0, 0.2], [3.0, 0.3]))
| f(x) = 1
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.9881102595664819, 0.4816273348099591, 0.7742585077942159]

julia> Metaheuristics.create_child(rand(3), ([-1, -2.0], [2.0, 0.2], [3.0, 0.3]))
| f(x) = [-1.0, -2.0]
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.23983577719146854, 0.3611544510766811, 0.7998754930109109]

julia> population = [ Metaheuristics.create_child(rand(2), (randn(2),  randn(2), rand(2))) for i = 1:100  ]
                           F space
          ┌────────────────────────────────────────┐
        2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠂⠀⠀⠀⠀⠀⡇⠈⡀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠘⠀⡇⠀⠀⠘⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⢀⠠⠐⠀⡇⠄⠁⠀⠀⠀⡀⠀⢁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠂⢈⠀⠈⡇⠀⡐⠃⠀⠄⠄⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠄⢐⠠⠀⡄⠀⠀⡇⠀⠂⠈⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠉⠉⠉⠉⠋⠉⠉⠉⠉⠉⠉⠙⢉⠉⠙⠉⠉⡏⠉⠉⠩⠋⠉⠩⠉⠉⠉⡉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
   f_2    │⠀⠀⠀⠀⠀⡀⠀⠀⠀⠄⠀⠀⡀⠀⠀⠂⠀⡇⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠐⡇⠠⠀⠀⠀⠈⢀⠄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠄⠀⡀⠀⠂⡇⠐⠘⠈⠂⠀⠈⡀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⡇⠀⠈⢀⠐⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⢁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          -3                                       4
                             f_1
```

"""
function create_child(x::Vector{Float64}, fResult::Float64;ε=0.0)
    return xf_indiv(x, fResult)
end

# constrained single objective
function create_child(x::Vector{Float64},
        fResult::Tuple{Float64,Array{Float64,1},Array{Float64,1}};
        ε=0.0
    )
    f, g, h = fResult
    isempty(g) && isempty(h) && (return xf_indiv(x, f))

    return xfgh_indiv(x, f, g, h; ε=ε)
end

# constrained multi-objective
function create_child(x::Vector{Float64},
        fResult::Tuple{Array{Float64,1},Array{Float64,1},Array{Float64,1}};
        ε = 0.0
    )
    f, g, h = fResult
    return xFgh_indiv(x, f, g, h;ε=ε)
end

##########################################################3
# parallel evaluation
##########################################################3
# single objective
function create_child(X::AbstractMatrix, fResult::AbstractVector;ε=0.0)
    size(X,1) != length(fResult) && error("Error in parallel evaluation: size(X,2) != length(f(X))")
    return [xf_indiv(X[i,:], fResult[i]) for i in 1:size(X,1)]
end

# unconstrained multi objective
function create_child(X::AbstractMatrix, fResult::AbstractMatrix;ε=0.0)
    size(fResult,2) > 1 && (return create_child(X, (fResult, zeros(0,0), zeros(0,0))))

    size(X,1) != size(fResult,1) && error("Error in parallel evaluation: size(X,2) != length(f(X))")

    return [xf_indiv(X[i,:], fResult[i]) for i in 1:size(X,1)]
end

# constrained single objective
function create_child(X::AbstractMatrix,
        fResult::Tuple{AbstractVector, T, T};
        ε=0.0
    ) where T <: AbstractMatrix{Float64}

    F, G, H = fResult

    size_ok = size(X,1) == length(F) #&& size(X,1) == size(G,1) && size(X,1) == size(H,1)
    !size_ok && error("Error in parallel evaluation: Size of X, F(X),G(X) or H(X) differs.")

    @assert isempty(G) || !isempty(G) && size(X,1) == size(G,1)
    @assert isempty(H) || !isempty(H) && size(X,1) == size(H,1)

    population = xfgh_indiv[]

    for i in 1:size(X,1)
        g = isempty(G) ? zeros(0) : G[i,:]
        h = isempty(H) ? zeros(0) : H[i,:]
        push!(population, xfgh_indiv(X[i,:], F[i], g, h;ε=ε))
    end
    population
end

# constrained multi-objective
function create_child(X::AbstractMatrix,
        fResult::Tuple{T, T, T};
        ε=0.0
    ) where T <: AbstractMatrix{Float64}

    F, G, H = fResult

    # single-objective
    size(F,2) == 1 && (return create_child(X, (F[:,1], G, H)))

    @assert size(X,1) == size(F,1)
    @assert isempty(G) || !isempty(G) && size(X,1) == size(G,1)
    @assert isempty(H) || !isempty(H) && size(X,1) == size(H,1)
    # !size_ok && error("Error in parallel evaluation: size(X,1) != size(F(X),1).")
    #&& error("Error in parallel evaluation: size(X,1) != size(G(X),1).")
    #&& error("Error in parallel evaluation: size(X,1) != size(H(X),1).")

    population = xFgh_indiv[]

    for i in 1:size(X,1)
        g = isempty(G) ? zeros(0) : G[i,:]
        h = isempty(H) ? zeros(0) : H[i,:]
        push!(population, xFgh_indiv(X[i,:], F[i,:], g, h;ε=ε))
    end
    population
end

##########################################################3

function generate_population(N::Int, problem;ε=0.0, parallel_evaluation = false)
    a = view(problem.bounds, 1, :)'
    b = view(problem.bounds, 2, :)'
    D = length(a)

    X = a .+ (b - a) .* rand(N, D)

    if problem.parallel_evaluation
        return create_solutions(X, problem; ε=ε)
    end

    population = [ create_solution(X[i,:], problem; ε=ε) for i in 1:N]

    return population
end

"""
    Metaheuristics.create_child(x, fx)

Constructor for a solution depending on the result of `fx`.

## Example

```
julia> import Metaheuristics

julia> Metaheuristics.create_child(rand(3), 1.0)
| f(x) = 1
| solution.x = [0.2700437125780806, 0.5233263210622989, 0.12871108215859772]

julia> Metaheuristics.create_child(rand(3), (1.0, [2.0, 0.2], [3.0, 0.3]))
| f(x) = 1
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.9881102595664819, 0.4816273348099591, 0.7742585077942159]

julia> Metaheuristics.create_child(rand(3), ([-1, -2.0], [2.0, 0.2], [3.0, 0.3]))
| f(x) = [-1.0, -2.0]
| g(x) = [2.0, 0.2]
| h(x) = [3.0, 0.3]
| x = [0.23983577719146854, 0.3611544510766811, 0.7998754930109109]

julia> population = [ Metaheuristics.create_child(rand(2), (randn(2),  randn(2), rand(2))) for i = 1:100  ]
                           F space
          ┌────────────────────────────────────────┐
        2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠂⠀⠀⠀⠀⠀⡇⠈⡀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠘⠀⡇⠀⠀⠘⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⢀⠠⠐⠀⡇⠄⠁⠀⠀⠀⡀⠀⢁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠂⢈⠀⠈⡇⠀⡐⠃⠀⠄⠄⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠄⢐⠠⠀⡄⠀⠀⡇⠀⠂⠈⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠉⠉⠉⠉⠋⠉⠉⠉⠉⠉⠉⠙⢉⠉⠙⠉⠉⡏⠉⠉⠩⠋⠉⠩⠉⠉⠉⡉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
   f_2    │⠀⠀⠀⠀⠀⡀⠀⠀⠀⠄⠀⠀⡀⠀⠀⠂⠀⡇⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠐⡇⠠⠀⠀⠀⠈⢀⠄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠄⠀⡀⠀⠂⡇⠐⠘⠈⠂⠀⠈⡀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⡇⠀⠈⢀⠐⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⢁⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       -3 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          └────────────────────────────────────────┘
          -3                                       4
                             f_1
```

"""
generateChild(x, fx) = create_child(x, fx)

function create_solution(x::AbstractVector, problem::Problem; ε=0.0)
    problem.f_calls += 1
    return create_child(x, problem.f(x), ε=ε)
end

function create_solutions(X::AbstractMatrix, problem::Problem; ε=0.0)
    problem.f_calls += size(X,1)
    if problem.parallel_evaluation
        return create_child(X, problem.f(X), ε=ε)
    end
    [create_child(X[i,:], problem.f(X[i,:]), ε=ε) for i in 1:size(X,1)]
end

# getters for the above structures

"""
    get_position(solution)

Get the position vector.
"""
get_position(solution::Solution) = solution.x

function positions(population::AbstractArray)
    if isempty(population)
        zeros(0,0)
    end

    D = length(get_position(population[1]))
    [ get_position(sol)[i] for sol in population, i in 1:D]
end
"""
    fval(solution)

Get the objective function value (fitness) of a solution.
"""
fval(solution::Solution) = solution.f
sum_violations(solution::T) where T <: AbstractConstrainedSolution = solution.sum_violations

"""
    gval(solution)

Get the inequality constraints of a solution.
"""
gval(solution::T) where T <: AbstractConstrainedSolution = solution.g

"""
    hval(solution)

Get the equality constraints of a solution.
"""
hval(solution::T) where T <: AbstractConstrainedSolution = solution.h

function fvals(population::AbstractArray{T}) where T <: AbstractMultiObjectiveSolution
    isempty(population) && (return zeros(0, 0))

    M = length(fval(population[1]))
    [ fval(sol)[i] for sol in population, i in 1:M]
end

fvals(population::AbstractArray{T}) where T <: AbstractSolution = fval.(population)

function gvals(population::AbstractArray{T}) where T <: AbstractConstrainedSolution
    if isempty(population)
        return zeros(0, 0)
    end

    M = length(gval(population[1]))
    [ gval(sol)[i] for sol in population, i in 1:M]
end

function hvals(population::AbstractArray{T}) where T <: AbstractConstrainedSolution
    if isempty(population)
        return zeros(0, 0)
    end

    M = length(hval(population[1]))
    [ hval(sol)[i] for sol in population, i in 1:M]
end

"""
    is_feasible(solution)
Returns `true` if solution is feasible, otherwise returns `false`.
"""
is_feasible(solution::T) where T <: AbstractConstrainedSolution = solution.is_feasible
is_feasible(solution::T) where T <: AbstractSolution = true

"""
    pareto_front(state::State)
Returns the non-dominated solutions in state.population.
"""
pareto_front(st::State) = pareto_front(st.population)

"""
    pareto_front(population::Array)
Returns non-dominated solutions.
"""
function pareto_front(population::AbstractArray{T}) where T <: AbstractMultiObjectiveSolution
    return fvals(get_non_dominated_solutions(population))
end


text:
# important: set d=1 in BridgeLandmarks.jl
using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD2

workdir = @__DIR__
println(workdir)
cd(workdir)

n = 3
nshapes = 1
T = 1.0

a = 2.0     # Hamiltonian kernel parameter (the larger, the stronger landmarks behave similarly)
c = 0.1     # multiplicative factor in kernel
γ = 1.0     # Noise level

Ptrue = MarslandShardlow(a, c, γ, 0.0, n)

q0 = [PointF(.5), PointF(1.0), PointF(1.3)]
p0 = zeros(PointF,n)
x0 = State(q0, p0)
xobs0 = x0.q
xobsT = [exp.(xobs0[i])  for i in 1:n ]

JLD2.@save "data_exp1-1D.jld2" xobs0 xobsT n x0 nshapes


text:
"""
Maximal ratio of circumcircle radius and length of shortest edge before  a triangle
is considered bad.
"""
beta = sqrt(2.)

minArea = 1e-5

"""
    compute_circumcenter_radius(mesh::Mesh, p::VertexIndex, r::VertexIndex, q::VertexIndex)

Computes the position and radius of the triangle circumcenter.

# Returns
* cx, cy, rs with ``c_x``, ``c_y`` the position and ``r_s`` the radius square of the circum center.
"""
function compute_circumcenter_radius(mesh::Mesh, p::VertexIndex, r::VertexIndex, q::VertexIndex)
  local pr = mesh.tesselation.vertices[p]
  local qr = mesh.tesselation.vertices[q]
  local rr = mesh.tesselation.vertices[r]
  local ax = pr[1]
  local ay = pr[2]
  local bx = qr[1]
  local by = qr[2]
  local cx = rr[1]
  local cy = rr[2]
  local sx = (ax^2.0*(by-cy) + ay^2.0*(by-cy) + bx^2.0*(cy-ay) + by^2.0*(cy-ay) + cx^2.0*(ay-by) + cy^2.0*(ay-by))/2.0
  local sy = (ax^2.0*(cx-bx) + ay^2.0*(cx-bx) + bx^2.0*(ax-cx) + by^2.0*(ax-cx) + cx^2.0*(bx-ax) + cy^2.0*(bx-ax))/2.0
  local a = ax*(by - cy) + bx * (cy-ay) + cx*(ay-by)
  local b = ax^2.0 * (bx*cy - by * cx) + ay^2.0 * (bx*cy - by*cx) + bx^2.0*(ay*cx-ax*cy) + by^2.0*(ay*cx-ax*cy) + cx^2.0*(ax*by-ay*bx) + cy^2.0*(ax*by-ay*bx)
  local rs = ((b/a + (sx*sx + sy*sy)/(a*a)))
  return sx/a, sy/a, rs
end

"""
    compute_offcenter(mesh:Mesh, p::VertexIndex, q::VertexIndex, r::VertexIndex)

Computes position of offcenter[^1] for given triangle ``p - q - r``.

# Remarks
* ``p - q`` must be the shortest edge of the triangle, i.e. the smallest angle is located at ``r``.

# References
[^1] Alper Üngör (2009). Off-centers: A new type of Steiner points for computing size-optimal quality-guaranteed Delaunay triangulations, Computational Geometry, *42* (2)

![SteinerPoint](SteinerPoint.svg)
"""
function compute_offcenter(mesh::Mesh, p::VertexIndex, q::VertexIndex, r::VertexIndex)
  local pr = mesh.tesselation.vertices[p]
  local qr = mesh.tesselation.vertices[q]
  local rr = mesh.tesselation.vertices[r]

  # position of bisection start of pq and circum center
  local bsr = 0.5 * (pr + qr)
  ccx, ccy, ccrs = compute_circumcenter_radius(mesh, p, q, r)

  # length of pq and radius of offcenter
  local lpq = norm(pr - qr)
  local ocr = beta * lpq * 0.999999

  # distance edge midpoint to center of offcircle ``c_2``
  local dobsoc = sqrt(ocr*ocr - (lpq/2.)*(lpq/2.) )

  # distance edge midpoint to center of outcircle
  local dbscc = norm(bsr - [ccx, ccy])
  if (dbscc <= dobsoc)
    # we take the out circle center
    return [ccx, ccy]
  else
    # we take off center

    # bisection: pq-midpoint
    local c2 = bsr + dobsoc * normalize([ccx, ccy] - bsr)
    local offcenter = c2 + ocr * normalize([ccx, ccy] - bsr)
    return offcenter
  end
end

"""
    vertex_encroaches_segment(mesh::Mesh, ei::Int, vertex::Point)

Checks if a potential vertex encroaches the segment denoted by `ei`.

# Remarks
* If edge `ei` is not a segment, i.e. not marked as a boundary, it always
  returns `false`
"""
function vertex_encroaches_segment(mesh::Mesh, ei::Int, vertex::Point)
    if (!isconstraint(mesh.tesselation, ei))
      return false
    end

    vio = org(mesh.tesselation, ei)
    vid = dest(mesh.tesselation, ei)
    v1 = mesh.tesselation.vertices[vio]
    v2 = mesh.tesselation.vertices[vid]
    mp = 0.5*(v1 + v2)
    dr = vertex - mp # from midpoint to new vertex
    dx = v1 - mp # from midpoint to end point

    return (dr[1]^2 + dr[2]^2) <= (dx[1]^2 + dx[2]^2)
  end

""""
    get_shortest_edge(mesh::Mesh, ai::VertexIndex, bi::VertexIndex, ci::VertexIndex)

Given a triangle with CCW-ordered vertices `ai` - `bi` - `ci`, it computes the
edge length and returns a sorted tuple `p` - `q` - `r` such that the triangle
``p - q - r`` has its shortest edge between `p` and `q`.
"""
function get_shortest_edge(mesh::Mesh, ai::VertexIndex, bi::VertexIndex, ci::VertexIndex)
  # nomenclature as in Bronstein et al..
  xa = mesh.tesselation.vertices[ai]
  xb = mesh.tesselation.vertices[bi]
  xc = mesh.tesselation.vertices[ci]
  vab = (xb - xa) # edge a
  vac = (xc - xa) # edge b
  vbc = (xb - xc) # edge a
  c = norm(vab)
  b = norm(vac)
  a = norm(vbc)
  #cosa = (b'*b + c'*c - a'*a)/(2. * b' * c)
  #cosb = (a'*a + c'*c - b'*b)/(2. * c' * a)
  #cosc = (a'*a + b'*b - c'*c)/(2. * a' * b)
  ordered = [(a, ai, bi, ci), (b, bi, ci, ai), (c, ci, ai, bi)]

  # order in place
  if (ordered[1][1] > ordered[2][1])
    ordered[1], ordered[2] = ordered[2], ordered[1]
  end
  if (ordered[2][1] > ordered[3][1])
    ordered[2], ordered[3] = ordered[3], ordered[2]
  end
  if (ordered[1][1] > ordered[2][1])
    ordered[1], ordered[2] = ordered[2], ordered[1]
  end

  return (ordered[1][2], ordered[1][3], ordered[1][4])
end

 """
    function refine_triangle(mesh::Mesh, tri::Triangle)

Refines the triangle located to the left of edge `ei` by adding a new vertex
at the offcenter or the midpoint of encroached edges.

# Remarks
* If a new vertex at the off-center of the triangle does not encroach any segments
  it is inserted at the off-center.
* Otherwise, a new vertex is inserted at the midpoint of each edge that the
  off-center would encroach on.
  try
"""
function refine_triangle(mesh::Mesh, tri::Triangle)

  # get ordered vertices
  p, q, r = get_shortest_edge(mesh, tri.vertices[1], tri.vertices[2], tri.vertices[3])

  # get offcenter
  offcenter = compute_offcenter(mesh, p, q, r)

  # get indices of all encroached segments
  indices = [i for i in 1:4:length(mesh.tesselation.edges)]
  encroachedIndices = filter(i->vertex_encroaches_segment(mesh, i, offcenter), indices)

  if isempty(encroachedIndices)
    # insert off-cente  r
    # we need to push the *unscaled point* so we push it straight to tesselation
    push_scaled!(mesh, offcenter')
  else
    # insert mid points
    for ei in encroachedIndices
      # get midpoint and insert it
      oi = org(mesh.tesselation, ei)
      di = dest(mesh.tesselation, ei)
      mp = 0.5 * (mesh.tesselation.vertices[oi] + mesh.tesselation.vertices[di])

      # we need to push the *unscaled point* so we push it straight to tesselation
      newVertexIndex = length(mesh.tesselation.vertices) + 1
      push_scaled!(mesh, mp')

      # we need to fix the face locations for the new faces.. we do that by
      # iterating over the star of the new vertex, starting at the segment
      startEdge = mesh.tesselation.vertexCache[newVertexIndex]
      while (!isconstraint(mesh.tesselation, startEdge))
        startEdge = onext(mesh.tesselation, startEdge)
      end
      faceIsInside = mesh.tesselation.edges[startEdge].constraintType == OnRightSide
      curEdge = startEdge
      while (true)
        # mark left face as inside/outside
        mesh.faceLocation[lface(mesh.tesselation, curEdge)] = faceIsInside
        curEdge = onext(mesh.tesselation, curEdge)
        if (curEdge == startEdge)
          break
        elseif (isconstraint(mesh.tesselation, curEdge))
          faceIsInside = !faceIsInside
        end
      end
    end
  end
end

"""
    check_triangle(mesh::Mesh, tri::Triangle)

Checks if the triangle satisfies the quality measure ``r/\bar{pq} \le \beta``,
with ``r`` the radius of the circumcircle and ``\bar{pq}`` the length of its
shortest edge.

# Returns
* `true`, if the triangle satisfies the quality measure.
"""
function check_triangle(mesh::Mesh, tri::Triangle)
  # TODO: we really need to get rid of the double computation of offcenter etc.
  #       maybe put a nullable in Triangle?

  # check if the area is below limit
  ra = mesh.tesselation.vertices[tri.vertices[1]]
  rb = mesh.tesselation.vertices[tri.vertices[2]]
  rc = mesh.tesselation.vertices[tri.vertices[3]]

  xa, ya = ra[1], ra[2]
  xb, yb = rb[1], rb[2]
  xc, yc = rc[1], rc[2]

  area = 0.5*abs((xa-xc)*(yb-ya)-(xa-xb)*(yc-ya))
  if area < minArea
    return true
  end

  # get ordered vertices and circumcircle
  p, q, r = get_shortest_edge(mesh, tri.vertices[1], tri.vertices[2], tri.vertices[3])
  sx, sy, rs = DelaunayMeshes.compute_circumcenter_radius(mesh, 1, 2, 3)

  lpq = norm(mesh.tesselation.vertices[p] - mesh.tesselation.vertices[q])
  return sqrt(rs)/lpq <= DelaunayMeshes.beta
end

"""
    refine_valid_triangles!(mesh::

Gets the first triangle that does not satisfy the quality measure or `null`.
"""
function refine_valid_triangles!(mesh::Mesh)
  foundBadTriangle = true
  counter = 0
  while foundBadTriangle && counter < 10000
    counter = counter + 1
    foundBadTriangle = false

    # get the first valid but bad triangle
    itState = start(mesh.tesselation)
    while (!done(mesh.tesselation, itState))
      tri, itState = next(mesh.tesselation, itState)
      if (mesh.faceLocation[tri.face] && !check_triangle(mesh, tri))

        # the triangle is bad and valid - we need to refine it
        refine_triangle(mesh, tri)
        foundBadTriangle = true

        # break out to outer loop
        break
      end
    end
  end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/division.jl
# Contains definitions of division.
#############################################################################

#=
@inline function div_MV(::Diff, x::MCNoGrad, y::MCNoGrad, z::Interval{Float64})
    if (0.0 < y.Intv.lo)
        cv, cv_grad = div_diffcv(x, y)
        cc, cc_grad = div_diffcv(-x, y)
        return MC{N,Diff}(cv, -cc, z, cv_grad, -cc_grad, x.cnst && y.cnst)
    end
    cv, cv_grad = div_diffcv(-x, -y)
    cc, cc_grad = div_diffcv(-x, y)
    return MC{N,Diff}(cv, -cc, z, cv_grad, -cc_grad, x.cnst && y.cnst)
end

@inline function div_kernel(t::Union{NS,MV}, x::MCNoGrad, y::MCNoGrad, z::Interval{Float64})
    if (x === y)
        zMC = one(MCNoGrad)
    else
        zMC = mult_kernel(t, x, inv(y), z)
    end
    return zMC
end

@inline function div_kernel(t::ANYRELAX, x::MCNoGrad, y::MCNoGrad, z::Interval{Float64})
    degen1 = (x.Intv.hi - x.Intv.lo == 0.0)
    degen2 = (y.Intv.hi - y.Intv.lo == 0.0)
    if x === y
        zMC = one(MCNoGrad)
    elseif  !degen1 || !degen2
        zMC = div_MV(t, x, y, z)
    else
        zMC = mult_kernel(t, x, inv(y), z)
    end
    return zMC
end
=#

@inline function /(t::ANYRELAX, x::MCNoGrad, y::MCNoGrad)
    if 0.0 ∉ y
        z, _ =  div_kernel(t, x, y, x.Intv/y.Intv)
        return z
    end
    return nan(MCNoGrad)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using CImGui
using CImGui: ImVec2

"""
    ShowDemoWindowColumns()
"""
function ShowDemoWindowColumns()
    CImGui.CollapsingHeader("Columns") || return

    CImGui.PushID("Columns")

    # basic columns
    if CImGui.TreeNode("Basic")
        CImGui.Text("Without border:")
        CImGui.Columns(3, "mycolumns3", false) # 3-ways, no border
        CImGui.Separator()
        for n = 0:13
            if CImGui.Selectable("Item $n")
                @info "Trigger Item $n | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            # if (CImGui.Button(label, ImVec2(-1,0))) {}
            CImGui.NextColumn()
        end
        CImGui.Columns(1)
        CImGui.Separator()

        CImGui.Text("With border:")
        CImGui.Columns(4, "mycolumns") # 4-ways, with border
        CImGui.Separator()
        CImGui.Text("ID"); CImGui.NextColumn()
        CImGui.Text("Name"); CImGui.NextColumn()
        CImGui.Text("Path"); CImGui.NextColumn()
        CImGui.Text("Hovered"); CImGui.NextColumn()
        CImGui.Separator()
        names = ["One", "Two", "Three"]
        paths = ["/path/one", "/path/two", "/path/three"]
        @cstatic selected=Cint(-1) begin
            for i = 0:2
                label = @sprintf "%04d" i
                if CImGui.Selectable(label, selected == i, CImGui.ImGuiSelectableFlags_SpanAllColumns)
                    selected = i
                end
                hovered = CImGui.IsItemHovered()
                CImGui.NextColumn()
                CImGui.Text(names[i+1]); CImGui.NextColumn()
                CImGui.Text(paths[i+1]); CImGui.NextColumn()
                CImGui.Text("hovered"); CImGui.NextColumn()
            end
        end
        CImGui.Columns(1)
        CImGui.Separator()
        CImGui.TreePop()
    end

    # create multiple items in a same cell before switching to next column
    if CImGui.TreeNode("Mixed items")
        CImGui.Columns(3, "mixed")
        CImGui.Separator()

        CImGui.Text("Hello")
        CImGui.Button("Banana")
        CImGui.NextColumn()

        CImGui.Text("ImGui")
        CImGui.Button("Apple")
        @cstatic foo=Cfloat(1.0) @c CImGui.InputFloat("red", &foo, 0.05, 0, "%.3f")
        CImGui.Text("An extra line here.")
        CImGui.NextColumn()

        CImGui.Text("Sailor")
        CImGui.Button("Corniflower")
        @cstatic bar=Cfloat(1.0) @c CImGui.InputFloat("blue", &bar, 0.05, 0, "%.3f")
        CImGui.NextColumn()

        CImGui.CollapsingHeader("Category A") && CImGui.Text("Blah blah blah")
        CImGui.NextColumn()
        CImGui.CollapsingHeader("Category B") && CImGui.Text("Blah blah blah")
        CImGui.NextColumn()
        CImGui.CollapsingHeader("Category C") && CImGui.Text("Blah blah blah")
        CImGui.NextColumn()
        CImGui.Columns(1)
        CImGui.Separator()
        CImGui.TreePop()
    end

    # word wrapping
    if CImGui.TreeNode("Word-wrapping")
        CImGui.Columns(2, "word-wrapping")
        CImGui.Separator()
        CImGui.TextWrapped("The quick brown fox jumps over the lazy dog.")
        CImGui.TextWrapped("Hello Left")
        CImGui.NextColumn()
        CImGui.TextWrapped("The quick brown fox jumps over the lazy dog.")
        CImGui.TextWrapped("Hello Right")
        CImGui.Columns(1)
        CImGui.Separator()
        CImGui.TreePop()
    end

    if CImGui.TreeNode("Borders")
        # NB: Future columns API should allow automatic horizontal borders.
        @cstatic h_borders=true v_borders=true begin
            @c CImGui.Checkbox("horizontal", &h_borders)
            CImGui.SameLine()
            @c CImGui.Checkbox("vertical", &v_borders)
            CImGui.Columns(4, C_NULL, v_borders)
            for i = 0:4*3-1
                (h_borders && CImGui.GetColumnIndex() == 0) && CImGui.Separator()
                CImGui.Text(@sprintf("%c%c%c", 'a'+i, 'a'+i, 'a'+i))
                CImGui.Text(@sprintf("Width %.2f\nOffset %.2f", CImGui.GetColumnWidth(), CImGui.GetColumnOffset()))
                CImGui.NextColumn()
            end
            CImGui.Columns(1)
            h_borders && CImGui.Separator()
        end
        CImGui.TreePop()
    end

    # scrolling columns
    # if CImGui.TreeNode("Vertical Scrolling")
    #     CImGui.BeginChild("##header", ImVec2(0, CImGui.GetTextLineHeightWithSpacing()+CImGui.GetStyle().ItemSpacing.y))
    #     CImGui.Columns(3)
    #     CImGui.Text("ID"); CImGui.NextColumn()
    #     CImGui.Text("Name"); CImGui.NextColumn()
    #     CImGui.Text("Path"); CImGui.NextColumn()
    #     CImGui.Columns(1)
    #     CImGui.Separator()
    #     CImGui.EndChild()
    #     CImGui.BeginChild("##scrollingregion", (0, 60))
    #     CImGui.Columns(3)
    #     for i = 0:9
    #         CImGui.Text(@sprintf("%04d", i))
    #         CImGui.NextColumn()
    #         CImGui.Text("Foobar")
    #         CImGui.NextColumn()
    #         CImGui.Text(@sprintf("/path/foobar/%04d/", i))
    #         CImGui.NextColumn()
    #     end
    #     CImGui.Columns(1)
    #     CImGui.EndChild()
    #     CImGui.TreePop()
    # end

    if CImGui.TreeNode("Horizontal Scrolling")
        CImGui.SetNextWindowContentSize((1500.0, 0.0))
        CImGui.BeginChild("##ScrollingRegion", ImVec2(0, CImGui.GetFontSize() * 20), false, CImGui.ImGuiWindowFlags_HorizontalScrollbar)
        CImGui.Columns(10)
        ITEMS_COUNT = 2000

        clipper = CImGui.Clipper(ITEMS_COUNT) # also demonstrate using the clipper for large list
        while CImGui.Step(clipper)
            s = CImGui.Get(clipper, :DisplayStart)
            e = CImGui.Get(clipper, :DisplayEnd)-1
            for i = s:e, j = 0:9
                CImGui.Text("Line $i Column $j...")
                CImGui.NextColumn()
            end
        end
        CImGui.Destroy(clipper)

        CImGui.Columns(1)
        CImGui.EndChild()
        CImGui.TreePop()
    end

    node_open = CImGui.TreeNode("Tree within single cell")
    CImGui.SameLine()
    ShowHelpMarker("NB: Tree node must be poped before ending the cell. There's no storage of state per-cell.")
    if node_open
        CImGui.Columns(2, "tree items")
        CImGui.Separator()
        if CImGui.TreeNode("Hello")
            CImGui.BulletText("Sailor")
            CImGui.TreePop()
        end
        CImGui.NextColumn()
        if CImGui.TreeNode("Bonjour")
            CImGui.BulletText("Marin")
            CImGui.TreePop()
        end
        CImGui.NextColumn()
        CImGui.Columns(1)
        CImGui.Separator()
        CImGui.TreePop()
    end
    CImGui.PopID()
end


text:
include("expression_iterator.jl")
using Compat
function lintpkg(pkg::AbstractString)
    if occursin("/", pkg) # pkg is a file path
        return LintResult(lintpkgforfile(pkg))
    end

    try
        p = Base.find_package(pkg)
        LintResult(lintpkgforfile(p))
    catch
        throw("cannot find package: " * pkg)
    end
end

"""
Lint the package for a file.
If file is in src, lint Package.jl and its includes.
If file is in test, lint runtests.jl and its includes.
If file is in base lint all files in base dir.
"""
function lintpkgforfile(path::AbstractString, ctx::LintContext=LintContext())
    path = abspath(path)
    if ispath(ctx.path)
        if Sys.iswindows()
            len = count(x -> x == '\\', path)
        else
            len = count(x -> x == '/', path) - 1
        end
        for _ = 1:len
            path, folder = splitdir(path)
            if folder == "src"
                file = joinpath(path, folder, basename(path) * ".jl")
                if ispath(file)
                    lintinclude(ctx, file)
                end
                break
            elseif folder == "test"
                file = joinpath(path, folder, "runtests.jl")
                if ispath(file)
                    lintinclude(ctx, file)
                end
                break
            elseif folder == "base"
                lintdir(joinpath(path, folder), ctx)
                break
            end
        end
    end
    finish(ctx)
    ctx.messages
end

function lintfile(f::AbstractString)
    if !isfile(f)
        throw("no such file exists")
    end
    lintfile(f, read(f, String))
end

function lintfile(f::AbstractString, code::AbstractString)
    ctx = LintContext(f)

    _lintstr(code, ctx)
    finish(ctx)
    msgs = ctx.messages

    # If we have an undeclared symbol, lint the package to try and resolve
    for message in msgs
        if message.code == :E321 # undeclared symbol
            ctx = LintContext(f)
            message = lintpkgforfile(file(ctx), ctx)
            break
        end
    end

    filter!(msg -> file(msg) == f, msgs)
    sort!(msgs)

    LintResult(msgs)
end

"Lint over each expression in each line.

Calls `lintexpr` over each parseable-parsed expression.
Each parse is called over each line."
function _lintstr(str::AbstractString, ctx::LintContext, lineoffset = 0)
    non_empty_lines=split(str, "\n", limit=0, keepempty=false)
    offset_where_last_expression_ends=nothing
    try
        for (ex, line_begin, line_end) in ExpressionIterator.each_expression_and_offset(str)
            # inform context of current line
            ctx.line = ctx.lineabs = line_end + lineoffset +1
            lintexpr(ex, ctx)
        end
    catch y
        # report an unexpected error
        # end-of-input and parsing errors are expected
        if typeof(y) != Meta.ParseError || y.msg != "end of input"
            msg(ctx, :E111, string(y))
        end
    end

end

"""
    lintstr(s::AbstractString)

Unlike other lint functions, like `lintpkg` and `lintfile`, `lintstr` does not
ignore any messages by default.
"""
function lintstr(str::AbstractString,
                 ctx::LintContext = (c = LintContext(); c.ignore = []; c))
    _lintstr(str, ctx)
    finish(ctx)
    ctx.messages
end

"""
Lint all .jl ending files at a given directory.
Will ignore LintContext file and already included files.
"""
function lintdir(dir::AbstractString, ctx::LintContext=LintContext())
    for file in readdir(dir)
        if endswith(file, ".jl")
            file = joinpath(dir, file)
            if !isdir(file)
                lintinclude(ctx, file)
            end
        end
    end
    finish(ctx)
    ctx.messages
end


text:
using LinearAlgebra
using StatsBase
using Distributions
cd(@__DIR__)
include("../examples/lutkepohl_example.jl")

h = 8

U = residuals(v)
T, K = size(U)
p = v.lags
nparam = size(v.B, 1)
Z = v.Z
ΣUml = ((T - nparam) / T) * v.ΣU

U2 = vcat(zeros(p + (h - p), K), U)
Ulag = VectorAutoregressions.lag_matrix(U2, h)
ZUlag = [Z Ulag]
AD = (ZUlag' * ZUlag) \ (ZUlag' * U)

E = U - ZUlag * AD
ΣEml = (E' * E) / T

# Test stat from Lütkepohl 2006, p.173
s = sqrt((K^4 * h^2 - 4) / (K^2 + (K^2 * h^2) - 5))
N = T - (K*p) - 1 - (K*h) - (K - K*h +1)/2
df1 = h * K^2
df2 = (N * s) - ((K^2 * h) / 2) + 1
Qss = ((det(ΣUml) / det(ΣEml))^(1 / s) - 1) * (df2 / df1)
pvalss = ccdf(FDist(df1, df2), Qss)

# Test stat from Killian and Lütkepohl 2017, p.54
Q = T * (K - tr(ΣUml \ ΣEml))
df = h * K^2
pval = ccdf(Chisq(df), Q)


text:
importall LightGraphs.SimpleGraphs
import LightGraphs.SimpleGraphs: fadj, badj, adj
struct DummySimpleGraph <: AbstractSimpleGraph end
@testset "SimpleGraphs" begin
    adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
    adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
    # specific concrete generators - no need for loop
    @test @inferred(eltype(SimpleGraph())) == Int
    @test @inferred(eltype(SimpleGraph(adjmx1))) == Int
    @test_throws ErrorException SimpleGraph(adjmx2)

    @test_throws ErrorException badj(DummySimpleGraph())

    @test @inferred(ne(SimpleGraph(PathDiGraph(5)))) == 4
    @test @inferred(!is_directed(SimpleGraph))

    @test @inferred(eltype(SimpleDiGraph())) == Int
    @test @inferred(eltype(SimpleDiGraph(adjmx2))) == Int
    @test @inferred(ne(SimpleDiGraph(PathGraph(5)))) == 8
    @test @inferred(is_directed(SimpleDiGraph))

    for gbig in [Graph(0xff), DiGraph(0xff)]
        @test @inferred(!add_vertex!(gbig))    # overflow
        @test @inferred(!add_vertices!(gbig, 10))
    end

    gdx = PathDiGraph(4)
    gx = SimpleGraph()
    for g in testgraphs(gx)
        T = eltype(g)
        @test sprint(show, g) == "empty undirected simple $T graph"
        @inferred(add_vertices!(g, 5))
        @test sprint(show, g) == "{5, 0} undirected simple $T graph"
    end
    gx = SimpleDiGraph()
    for g in testdigraphs(gx)
        T = eltype(g)
        @test sprint(show, g) == "empty directed simple $T graph"
        @inferred(add_vertices!(g, 5))
        @test sprint(show, g) == "{5, 0} directed simple $T graph"
    end

    gx = PathGraph(4)
    for g in testgraphs(gx)
        @test @inferred(vertices(g)) == 1:4
        @test Edge(2,3) in edges(g)
        @test @inferred(nv(g)) == 4
        @test @inferred(fadj(g)) == badj(g) == adj(g) == g.fadjlist
        @test @inferred(fadj(g,2)) == badj(g,2) == adj(g,2) == g.fadjlist[2]

        @test @inferred(has_edge(g, 2, 3))
        @test @inferred(has_edge(g, 3, 2))

        gc = copy(g)
        @test @inferred(add_edge!(gc, 4=>1)) && gc == CycleGraph(4)
        @test @inferred(has_edge(gc, 4=>1)) && has_edge(gc, 0x04=>0x01)
        gc = copy(g)
        @test @inferred(add_edge!(gc, (4,1))) && gc == CycleGraph(4)
        @test @inferred(has_edge(gc, (4,1))) && has_edge(gc, (0x04, 0x01))
        gc = copy(g)
        @test add_edge!(gc, 4, 1) && gc == CycleGraph(4)

        @test @inferred(in_neighbors(g, 2)) == @inferred(out_neighbors(g, 2)) == @inferred(neighbors(g,2)) == [1,3]
        @test @inferred(add_vertex!(gc))   # out of order, but we want it for issubset
        @test @inferred(g ⊆ gc)
        @test @inferred(has_vertex(gc, 5))

        @test @inferred(ne(g)) == 3

        @test @inferred(rem_edge!(gc, 1, 2)) && @inferred(!has_edge(gc, 1, 2))
        ga = @inferred(copy(g))
        @test @inferred(rem_vertex!(ga, 2)) && ne(ga) == 1
        @test @inferred(!rem_vertex!(ga, 10))

        @test @inferred(zero(g)) == SimpleGraph{eltype(g)}()

        # concrete tests below

        @test @inferred(eltype(g)) == eltype(fadj(g,1)) == eltype(nv(g))
        T = @inferred(eltype(g))
        @test @inferred(nv(SimpleGraph{T}(6))) == 6

        @test @inferred(eltype(SimpleGraph(T))) == T
        @test @inferred(eltype(SimpleGraph{T}(adjmx1))) == T

        ga = SimpleGraph(10)
        @test @inferred(eltype(SimpleGraph{T}(ga))) == T

        for gd in testdigraphs(gdx)
            U = eltype(gd)
            @test @inferred(eltype(SimpleGraph(gd))) == U
        end

        @test @inferred(edgetype(g)) == SimpleGraphEdge{T}
        @test @inferred(copy(g)) == g
        @test @inferred(!is_directed(g))

        e = first(edges(g))
        @test @inferred(has_edge(g, e))
    end

    gdx = PathDiGraph(4)
    for g in testdigraphs(gdx)
        @test @inferred(vertices(g)) == 1:4
        @test Edge(2,3) in edges(g)
        @test !(Edge(3,2) in edges(g))
        @test @inferred(nv(g)) == 4
        @test @inferred(fadj(g)[2]) == fadj(g, 2) == [3]
        @test @inferred(badj(g)[2]) == badj(g, 2) == [1]
        @test_throws MethodError adj(g)

        @test @inferred(has_edge(g, 2, 3))
        @test @inferred(!has_edge(g, 3, 2))

        gc = copy(g)
        @test @inferred(add_edge!(gc, 4=>1)) && gc == CycleDiGraph(4)
        @test @inferred(has_edge(gc, 4=>1)) && has_edge(gc, 0x04=>0x01)
        gc = copy(g)
        @test @inferred(add_edge!(gc, (4,1))) && gc == CycleDiGraph(4)
        @test @inferred(has_edge(gc, (4,1))) && has_edge(gc, (0x04, 0x01))
        gc = @inferred(copy(g))
        @test @inferred(add_edge!(gc, 4, 1)) && gc == CycleDiGraph(4)

        @test @inferred(in_neighbors(g, 2)) == [1]
        @test @inferred(out_neighbors(g, 2)) == @inferred(neighbors(g,2)) == [3]
        @test @inferred(add_vertex!(gc))   # out of order, but we want it for issubset
        @test @inferred(g ⊆ gc)
        @test @inferred(has_vertex(gc, 5))

        @test @inferred(ne(g)) == 3

        @test @inferred(!rem_edge!(gc, 2, 1))
        @test @inferred(rem_edge!(gc, 1, 2)) && @inferred(!has_edge(gc, 1, 2))
        ga = @inferred(copy(g))
        @test @inferred(rem_vertex!(ga, 2)) && ne(ga) == 1
        @test @inferred(!rem_vertex!(ga, 10))

        @test @inferred(zero(g)) == SimpleDiGraph{eltype(g)}()

        # concrete tests below

        @test @inferred(eltype(g)) == eltype(@inferred(fadj(g,1))) == eltype(nv(g))
        T = @inferred(eltype(g))
        @test @inferred(nv(SimpleDiGraph{T}(6))) == 6

        @test @inferred(eltype(SimpleDiGraph(T))) == T
        @test @inferred(eltype(SimpleDiGraph{T}(adjmx2))) == T

        ga = SimpleDiGraph(10)
        @test @inferred(eltype(SimpleDiGraph{T}(ga))) == T

        for gu in testgraphs(gx)
            U = @inferred(eltype(gu))
            @test @inferred(eltype(SimpleDiGraph(gu))) == U
        end

        @test @inferred(edgetype(g)) == SimpleDiGraphEdge{T}
        @test @inferred(copy(g)) == g
        @test @inferred(is_directed(g))

        e = first(@inferred(edges(g)))
        @test @inferred(has_edge(g, e))
    end

    gdx = CompleteDiGraph(4)
    for g in testdigraphs(gdx)
        @test rem_vertex!(g, 2)
        @test nv(g) == 3 && ne(g) == 6
    end

end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Test
using MARY_fRG.Triangulated
using MARY_fRG.Fermi
using MARY_fRG.Fermi.Surface
using MARY_fRG.Drawers

@testset "正方格子的费米面" begin
    qs = common_square_lattice(0.)
    ltris, ladjs = split_square(qs.brillouin, 100)
    sf = const_energy_line(ltris, ladjs, 0., qs.dispersion[1])
    plt = draw_lines(qs.brillouin.edges)
    draw_lines!(plt, sf; lcidx=ones(Int64, length(sf)))
    #draw_polygon!(plt, ltris)
    sf2 = const_energy_line(ltris, ladjs, -0.2, qs.dispersion[1])
    draw_lines!(plt, sf2; lcidx=2ones(Int64, length(sf2)))
    @savepng plt "surface1"
end

@testset "三角格子的费米面" begin
    ts = common_triangle_lattice(0.)
    ltris, ladjs = split_hexagon(ts.brillouin, 100)
    sf = const_energy_line(ltris, ladjs, 0., ts.dispersion[1])
    plt = draw_lines(ts.brillouin.edges)
    draw_lines!(plt, sf; lcidx=ones(Int64, length(sf)))
    sf2 = const_energy_line(ltris, ladjs, 2.0, ts.dispersion[1])
    draw_lines!(plt, sf2; lcidx=2ones(Int64, length(sf2)))
    @savepng plt "surface2"
end


text:
# Helper types for implementing `push`

struct SingletonVector{T} <: AbstractVector{T}
    value::T
end

Base.size(::SingletonVector) = (1,)

function Base.getindex(v::SingletonVector, i::Integer)
    @boundscheck i == 1 || throw(BoundsError(v, i))
    return v.value
end

struct SingletonDict{K, V} <: AbstractDict{K, V}
    key::K
    value::V
end

Base.iterate(d::SingletonDict) = (d.key => d.value, nothing)
Base.iterate(d::SingletonDict, ::Nothing) = nothing

function Base.getindex(d::SingletonDict{K}, key::K) where K
    @boundscheck d.key == key || throw(BoundsError(d, key))
    return d.value
end


text:
@testset "radial map" begin
    γ = 1.0
    λ = 0.2
    δ = 1e-4
    κ = 5.0

    M = RadialMap(4,2, γ=γ, λ=λ, δ=δ, κ=κ)
    @test size(M) == (4,2)

    @test M.γ==1.0
    @test M.λ==0.2
    @test M.δ==1e-4
    @test M.κ==5.0

    for i=1:4
        @test size(M.C[i]) ==(i,2)
    end

    M = RadialMap(2,2)
    M.C[1].coeff[1] .= ones(2+3)
    M.C[2].coeff[1] .= 5*ones(2+1)
    M.C[2].coeff[2] .= ones(2+3)

    x = randn()
    y = randn(2)

    C11 = RadialMapComponent(1,2, [zeros(4)], [ones(4)], [ones(2+3)])
    C22 = RadialMapComponent(2,2, [zeros(2),zeros(4)], [ones(2), ones(4)], [5*ones(2+1), ones(2+3)])

    @test isapprox(M([x; x]), [C11([x]); C22([x; x])], atol = 1e-10)
    @test isapprox(M(y), [C11([y[1]]); C22(y)], atol = 1e-10)

    C = RadialMapComponent(2,2)
    A = randn(8)
    modifycoeff!(C, A)

    Cprime = ForwardDiff.gradient(x->C(x), y)

    @test isapprox(Cprime[1], C.coeff[1][1] + C.coeff[1][2]*rbf′(0.0, 1.0)(y[1]) + C.coeff[1][3]*rbf′(0.0, 1.0)(y[1]), atol = 1e-10)
    @test isapprox(Cprime[2], 0.0*C.coeff[2][1] + C.coeff[2][2]*ψ₀′(0.0, 1.0)(y[2]) + C.coeff[2][3]*rbf(0.0, 1.0)(y[2]) + C.coeff[2][4]*rbf(0.0, 1.0)(y[2]) + C.coeff[2][5]*ψpp1′(0.0, 1.0)(y[2]), atol = 1e-10)
end

@testset "Sparse radial map I" begin
    γ = 1.0
    λ = 0.2
    δ = 1e-4
    κ = 5.0
    # p = [[2];[-1; 1]; []

    M = SparseRadialMap(4,fill(2,4), γ=γ, λ=λ, δ=δ, κ=κ)
    @test size(M) == (4,[fill(2,i) for i=1:4])

    @test M.γ==1.0
    @test M.λ==0.2
    @test M.δ==1e-4
    @test M.κ==5.0

    for i=1:4
        @test size(M.C[i]) ==(i,fill(2,i))
    end
end

@testset "Sparse radial map II" begin
    p = [Int64[-1], [0; -1], [2; 0; 1], [-1 ;2;-1; 0]]
    M = SparseRadialMap(4,p)

    for i=1:4
        @test M.p[i]==p[i]
    end

    z =randn(4)

    @test M(z)[1]==z[1]

    M.C[2].coeff[1] =randn(1)

    M.C[3].ξ[1] = randn(2)
    M.C[3].σ[1] = rand(2)
    M.C[3].coeff[1] = randn(3)

    M(z)
    M.C[3].coeff[2] = randn(1)

    M.C[3].ξ[3] = randn(3)
    M.C[3].σ[3] = rand(3)
    M.C[3].coeff[3] = randn(4)

    M(z)

    M.C[4].ξ[2] = randn(2)
    M.C[4].σ[2] = rand(2)
    M.C[4].coeff[2] = randn(3)

    M.C[4].coeff[4] = randn(2)

    @test norm(M(z)[2] - ( ui(p[2][1],Float64[],Float64[],M.C[2].coeff[1])(z[1]) + z[2]))<1e-12

    @test norm(M(z)[3] -( ui(p[3][1],M.C[3].ξ[1],M.C[3].σ[1],M.C[3].coeff[1])(z[1])+
                     ui(p[3][2],M.C[3].ξ[2],M.C[3].σ[2],M.C[3].coeff[2])(z[2])+
                     uk(p[3][3],M.C[3].ξ[3],M.C[3].σ[3],M.C[3].coeff[3])(z[3])))<1e-12

    @test norm(M(z)[4] -( ui(p[4][2],M.C[4].ξ[2],M.C[4].σ[2],M.C[4].coeff[2])(z[2])+
                    uk(p[4][4],M.C[4].ξ[4],M.C[4].σ[4],M.C[4].coeff[4])(z[4])))<1e-12

end


text:
using Documenter

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

using TestContexts

makedocs(
    sitename = "TestContexts.jl",
    modules = [TestContexts],
    authors = "Oren Ben-Kiki",
    clean = true,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = ["index.md"],
)


text:
module ReactionEngine
# Write your package code here.
include("Inspect.jl")
include("ThermoProbe.jl")
include("Plug.jl")
include("Cstr.jl")
include("Batch.jl")
include("Equil.jl")
include("Transport.jl")
include("Sens.jl")
export inspect, thermoprobe, cstr, plug, batch, equilibrate, transport_properties, global_sensitivity

"""
inspect(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
inspect(input_file::AbstractString, lib_dir::AbstractString) = Inspect.run(input_file, lib_dir)

"""
thermoprobe(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
thermoprobe(input_file::AbstractString, lib_dir::AbstractString ) = ThermoProbe.run(input_file,lib_dir)

"""
cstr(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
cstr(input_file::AbstractString, lib_dir::AbstractString) = Cstr.cstr(input_file,lib_dir)

"""
plug(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
plug(input_file::AbstractString, lib_dir::AbstractString) = Plug.plug(input_file,lib_dir)

"""
batch(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
batch(input_file::AbstractString, lib_dir::AbstractString) = Batch.batch(input_file, lib_dir)

"""
equilibrate(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
equilibrate(input_file::AbstractString, lib_dir::AbstractString) = Equil.equil(input_file, lib_dir)

"""
transport_properties(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The  directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
transport_properties(input_file::AbstractString, lib_dir::AbstractString) = Transport.run(input_file,lib_dir)

"""
global_sensitivity(input_file::AbstractString, lib_dir::AbstractString)
- input_file : Input xml file (String)
- lib_dir : The directory where the data files are present.
            In the present version, the data files are present in the lib directory within the test folder.
            The user must specify the relative path to the lib directory. The name of the lib directory need not be
            "lib"
"""
global_sensitivity(input_file::AbstractString, lib_dir::AbstractString) = Sens.gsa_xml(input_file,lib_dir)

end


text:
module PGA2D

using CliffordAlgebras

import Base.join, Base.show

export pga2d, PGA2DT, PGA2DMV
export e0,e1,e2, e01, e12, e20, e012, 𝐈
export direction, point, line
export point_coordinates
export direction_coordinates
export line_coordinates
export is_point, is_direction
export is_line, is_motor
export normalize
export join_pp
export meet_ll
export l_ortho_lp
export p_proj_lp
export l_para_lp
export d_ortho_l
export dist_pp
export dist_ll
export angle_ll
export angle_dl
export dist_orient_lp
export l_bisect_ll
export l_bisect_pp
export p_bisect_pp
export reflect_l
export area_ppp
export area_loop
export length_loop
export motor_pa
export motor_ll
export move_pa
export move_ll
export move_m

const pga2d = CliffordAlgebra(:PGA2D)
const PGA2DT = typeof(pga2d)
const PGA2DMV = MultiVector{PGA2DT}

const e0 = pga2d.e0
const e1 = pga2d.e1
const e2 = pga2d.e2

const e01 = pga2d.e0e1
const e12 = pga2d.e1e2
const e20 = pga2d.e2e0

const e012 = pga2d.e1e2e0
const 𝐈 = e012

"""
    point(x::Real, y::Real)

Constructs a MultiVector that encodes a point (x,y).
"""
point(x::Real, y::Real) = direction(x, y) + e12

"""
    direction(x::Real, y::Real)

Constructs a MultiVector that encodes a direction (x,y).
"""
direction(x::Real, y::Real) = x * e20 + y * e01

"""
    line(a::Real, b::Real, c::Real)

Constructs a MultiVector that encodes a line ax + bx + c = 0.
"""
line(a::Real, b::Real, c::Real) = a * e1 + b * e2 + c * e0

"""
    point_coordinates(P::PGA2DMV)

Extracts the point coordinates (x,y) from a MultiVector P. If the MultiVector does not encode a point a DomainError is thrown.
"""
function point_coordinates(P::PGA2DMV)
    w = P.e1e2
    if isgrade(P,2) && !iszero(w)
        wi = inv(w)
        return wi.*(P.e2e0, P.e0e1)
    else
        throw(DomainError(P, " is not a point."))
    end
end

"""
    direction_coordinates(d::PGA2DMV)

Extracts the direction coordinates (x,y) from a MultiVector P. If the MultiVector does not encode a direction a DomainError is thrown.
"""
function direction_coordinates(d::PGA2DMV)
    w = d.e1e2
    if isgrade(d,2) && iszero(w)
        return (d.e2e0, d.e0e1)
    else
        throw(DomainError(d, " is not a direction."))
    end
end

"""
    line_coordinates(l::PGA2DMV)

Extracts the line coordinates (a,b,c) from a MultiVector l. If the MultiVector does not encode a line a DomainError is thrown.
"""
function line_coordinates(l::PGA2DMV)
    if isgrade(l,1) && (!iszero(l.e1) || !iszero(l.e2))
        return (l.e1, l.e2, l.e0)
    else
        throw(DomainError(l, " is not a line."))
    end
end

"""
    is_point(P::PGA2DMV)

Returns true iff the MultiVector P encodes a point.
"""
function is_point(P::PGA2DMV)
    isgrade(P,2) && !iszero(P.e1e2)
end

"""
    is_direction(d::PGA2DMV)

Returns true iff the MultiVector d encodes a direction.
"""
function is_direction(d::PGA2DMV)
    isgrade(d,2) && iszero(d.e1e2) && !iszero(d)
end

"""
    is_line(l::PGA2DMV)

Returns true iff the MultiVector l encodes a line.
"""
function is_line(l::PGA2DMV)
    isgrade(l,1) && (!iszero(l.e1) || !iszero(l.e2))
end

"""
    is_motor(m::PGA2DMV)

Returns true iff the MultiVector m encodes a motor.
"""
function is_motor(m::PGA2DMV)
    iszero(m.e1) && iszero(m.e2) && iszero(m.e0) && iszero(m.e1e2e0) && !iszero(m)
end

function show(io::IO, x::PGA2DMV)
    if is_point(x)
        print(io, "point", point_coordinates(x), " ≜ ")
    elseif is_direction(x)
        print(io, "direction", direction_coordinates(x), " ≜ ")
    elseif is_line(x)
        print(io, "line", line_coordinates(normalize(x)), " ≜ ")
    elseif is_motor(x)
        m = normalize(x)
        P = grade(m,2)
        if is_point(P)
            orientation = sign(P.e1e2)
            α = 2 * orientation * atan(norm(P), scalar(m))
            print(io, "motor: rotation $(α/π)π around $P", " ≜ ")
        elseif is_direction(P)
            d = 2 * sqrt(P.e0e1^2 + P.e2e0^2)
            d_ortho = (P.e2e0,-P.e0e1).*(2/d)
            print(io, "motor: translation $d direction $(d_ortho)", " ≜ ")
        end
    end
    show_multivector(io, x)
end

"""
    normalize(x::PGA2DMV)

Scales the MultiVector x so that it has unit norm.
"""
normalize(x::PGA2DMV) = x / norm(x)

"""
    meet_ll(l1::PGA2DMV, l2::PGA2DMV)

Calculates the point of intersection of the two line-encoding MultiVectors l1 and l2.
If the lines are parallel the result is a direction parallel to the lines. If the lines coincides, the result is 0.
"""
meet_ll(l1::PGA2DMV, l2::PGA2DMV) = l1 ∧ l2

"""
    join_pp(P1::PGA2DMV, P2::PGA2DMV)

Calculates the line through the points P1 and P2. If P1 and P2 coincide, the result is 0.
"""
join_pp(P1::PGA2DMV, P2::PGA2DMV) = P1 ∨ P2

"""
    l_ortho_lp(l::PGA2DMV, P::PGA2DMV)

Calculates the line that is orthogonal to l and passes through P.
"""
l_ortho_lp(l::PGA2DMV, P::PGA2DMV) = l ⋅ P

"""
    p_proj_lp(l::PGA2DMV, P::PGA2DMV)

Orthogonally projects the point P onto the line l.
"""
p_proj_lp(l::PGA2DMV, P::PGA2DMV) = (l ⋅ P) * l

"""
    l_para_lp(l::PGA2DMV, P::PGA2DMV)

Constructs a parallel line to the line l in the point P.
"""
l_para_lp(l::PGA2DMV, P::PGA2DMV) = (l ⋅ P ) * P

"""
    d_ortho_l(l::PGA2DMV)

Calculates the direction orthogonal to l.
"""
d_ortho_l(l::PGA2DMV) = l * 𝐈

"""
    dist_pp(P1::PGA2DMV, P2::PGA2DMV)

Calculates the Euclidean distance between P1 and P2.
"""
dist_pp(P1::PGA2DMV, P2::PGA2DMV) = norm(P1 ∨ P2) / (norm(P1) * norm(P2))

"""
    dist_ll(l1::PGA2DMV, l2::PGA2DMV)

Calculates the orthogonal Euclidean distance between parallel lines l1 and l2.
"""
dist_ll(l1::PGA2DMV, l2::PGA2DMV) = norm(dual(l1 ∧ l2)) / (norm(l1) * norm(l2))

"""
    angle_ll(l1::PGA2DMV, l2::PGA2DMV)

Calculates the angle between two lines l1 and l2.
"""
angle_ll(l1::PGA2DMV, l2::PGA2DMV) = asin( norm(l1 ∧ l2) / (norm(l1) * norm(l2)))

"""
    angle_dl(d::PGA2DMV, l2::PGA2DMV)

Calculates the angle between the direction d and the line l.
"""
angle_dl(d::PGA2DMV, l::PGA2DMV) = asin( norm(dual(d ∧ l)) / (norm(d) * norm(l)))

"""
    dist_orient_lp(l::PGA2DMV, P::PGA2DMV)

Calculates the oriented Euclidean distance between the line l and the point P.
"""
dist_orient_lp(l::PGA2DMV, P::PGA2DMV) = scalar(P ∨ l) / (norm(P) * norm(l))

"""
    l_bisect_ll(l1::PGA2DMV, l2::PGA2DMV)

Calculates both angle-bisecting lines between l1 and l2 and return them in a tuple of MultiVectors.
"""
function l_bisect_ll(l1::PGA2DMV, l2::PGA2DMV)
    l1n = normalize(l1)
    l2n = normalize(l2)
    (l1n + l2n, l1n - l2n)
end

"""
    l_bisect_pp(P1::PGA2DMV, P2::PGA2DMV)

Calculates the orthogonal bisecting line between the points P1 and P2.
"""
function l_bisect_pp(P1::PGA2DMV, P2::PGA2DMV)
    P1n = normalize(P1)
    P2n = normalize(P2)
    (P1n + P2n) ×₋ (P1n ∨ P2n)
end

"""
    p_bisect_pp(P1::PGA2DMV, P2::PGA2DMV)

Calculates the midpoint between P1 and P2.
"""
p_bisect_pp(P1::PGA2DMV, P2::PGA2DMV) = normalize(P1) + normalize(P2)

"""
    reflect_l(x::PGA2DMV, l::PGA2DMV)

Reflects any MultiVector encoded object x on the line l.
"""
reflect_l(x::PGA2DMV, l::PGA2DMV) = l * x * l

"""
    area_ppp(P1::PGA2DMV, P2::PGA2DMV, P3::PGA2DMV)

Calculates the oriented area of the triangle P1,P2,P3.
"""
area_ppp(P1::PGA2DMV, P2::PGA2DMV, P3::PGA2DMV) = scalar((P1 ∨ P2 ∨ P3) / (2 * norm(P1) * norm(P2) * norm(P3)))

"""
    area_loop(points)

Calculates the oriented area of the loop defined by the iteratable points.
"""
function area_loop(points)
    area = 0
    for k = 2:(length(points)-1)
        area += triangle_area(points[1],points[k],points[k+1])
    end
    area
end

"""
    length_loop(points)

Calculates the length of the loop defined by the iteratable points.
"""
function length_loop(points)
    Pp = last(points)
    Ppn = normalize(Pp)
    l = zero(eltype(Ppn))
    for P in points
        Pn = P / norm(P)
        l += norm( Ppn ∨ Pn )
        Ppn = Pn
    end
    l
end

"""
    motor_pa(P::PGA2DMV, α::Real)

Calculates the motor for a rotation by the angle α around P. If P is a direction, then the motor describes a translation by α orthogonal to the direction.
"""
function motor_pa(P::PGA2DMV, α::Real)
    if is_direction(P)
        1 + α/(2*sqrt(P.e0e1^2+P.e2e0^2)) * P
    else
        exp(α/2*normalize(P))
    end
end

"""
    motor_ll(l1::PGA2DMV, l2::PGA2DMV)

Calculates the motor that maps the lines l1 to l2.
"""
function motor_ll(l1::PGA2DMV, l2::PGA2DMV)
    l1n = normalize(l1)
    l2n = normalize(l2)
    l2l1 = l2n * l1n
    l2l1n = normalize(l2l1)
    1 + l2l1n
end

"""
    move_pa(x::PGA2DMV, P::PGA2DMV, α::Real)

Moves x with motor_pa(x,P,α).
"""
move_pa(x::PGA2DMV, P::PGA2DMV, α::Real) = motor_pa(P,α) ≀ x

"""
    move_ll(x::PGA2DMV, l1::PGA2DMV, l2::PGA2DMV)

Moves x with motor_ll(l1,l2).
"""
move_ll(x::PGA2DMV, l1::PGA2DMV, l2::PGA2DMV) = motor_ll(l1,l2) ≀ x

"""
    move_m(x::PGA2DMV, m::PGA2DMV)

Moves x with the motor m.
"""
move_m(x::PGA2DMV, m::PGA2DMV) = m ≀ x

using Plots

@recipe function multivector_plot_recipe(mv::PGA2DMV)
    if is_point(mv)
        seriestype := :path
        markershape --> :circle
        markersize --> 5
        (x,y) = point_coordinates(mv)
        [x],[y]
    elseif is_line(mv)
        seriestype := :straightline
        (a,b,c) = line_coordinates(mv)
        if abs(a) > abs(b)
            [-c/a, -c/a+b],[0,-a]
        else
            [0,-b],[-c/b, -c/b+a]
        end
    else
        nothing
    end
end

end


text:
module Jobs

using Distributed
using OrderedCollections

export launch_workers,
    current_results, current_times, current_time_info, run_on_workers

capacity = 10000
const jobs = RemoteChannel(()->Channel{Tuple{Any, Function, Any}}(capacity))
const results = RemoteChannel(()->Channel{Pair{Tuple{Any, Any}, Tuple}}(
                    capacity))
const collected = OrderedDict()

function kill_workers(; waitfor=30)
    ws = workers()
    (length(ws) == 1 && 1 in ws) && return
    rmprocs(ws, waitfor=waitfor)
end

function clear_jobs()
    while isready(jobs)
        take!(jobs)
    end
    nothing
end

function take_results()
    while isready(results)
        (job_id, job_args), (
            result, t, bytes, gctime, memallocs) = take!(results)
        collected[job_id][job_args] = (result, t, bytes, gctime, memallocs)
    end
    nothing
end

function launch_workers(n::Integer)
    num_new = n - length(setdiff!(Set(workers()), [1]))
    # Create worker processes
    new_workers = addprocs(num_new; exeflags="--project")
    length(workers()) == n || error("didn't start the right number of workers")
    1 in workers() && error("main process is listed as a worker")
    nothing
end

function current_results(job_id)
    take_results()
    OrderedDict(
         k => val
         for (k, (val,)) in collected[job_id]
    )
end

function current_times(job_id)
    take_results()
    OrderedDict(
         k => time
         for (k, (val, time)) in collected[job_id]
    )
end

function current_time_info(job_id)
    take_results()
    OrderedDict(
         k => v[2:end]
         for (k, v) in collected[job_id]
    )
end

function fetch_results_internal(job_id, plot_i, default::Tuple{Float64,Float64,Float64,Float64}=(1.0,1.0,1.0,1.0))
    x_arr = x_axis_for_plot(plot_i)
    y_arr_dict = Dict{Int, Vector{Tuple{Float64,Float64,Float64,Float64}}}()
    results = current_results(job_id)
    for ((res_plot_i, _, d, e, th,  _, _), val) in results
        res_plot_i == plot_i || continue
        if !(d in keys(y_arr_dict))
            y_arr_dict[d] = repeat([default], length(x_arr))
        end
        if plot_i == 1
            i = indexin([e], x_arr)[1]
        elseif 2 <= plot_i <= 3
            i = indexin([th], x_arr)[1]
        elseif plot_i == 4
            i = indexin([th], x_arr)[1]
        end
        i === nothing || (y_arr_dict[d][i] = val)
    end
    dists = sort!(collect(keys(y_arr_dict)))
    y_arrs = [y_arr_dict[d] for d in dists]
    x_arr, dists, y_arrs
end

function run_on_workers(job_id, job_f, job_args_list)
    job_id in keys(collected) && error("job_id <$job_id> has already been used")

    1 in workers() && error("run launch_workers(n) first")

    # Set up workerpool
    wp = CachingPool(workers())

    # Ensure do_work is defined
    @everywhere function _do_work(job_id, job_f, job_args, results)
        print("Worker $(myid()) begins: ", job_args, "\n")
        flush(stdout)
        val, t, bytes, gctime, memallocs = @timed job_f(job_args...)
        put!(results,
             (job_id, job_args) => (val, t, bytes, gctime, memallocs))

        #while isready(jobs)
        #    job_id, job_f, job_args = take!(jobs)
        #    # Run the job
        #    val, t, bytes, gctime, memallocs = @timed job_f(job_args...)
        #    put!(results,
        #         (job_id, job_args) => (val, t, bytes, gctime, memallocs))
        #end
        print("Worker $(myid()) is done: ", job_args, "\n")
        flush(stdout)
    end

    # Make jobs
    collected[job_id] = OrderedDict()
    #function make_jobs()
    #    for job_args in job_args_list
    #        put!(jobs, (job_id, job_f, job_args))
    #    end
    #end
    ##@async
    #make_jobs()

    ## Start worker tasks
    #@sync for pid in workers()
    #    @async begin
    #        remotecall_wait(Main._do_work, pid, jobs, results)
    #    end
    #end
    # Start worker tasks
    @sync for job_args in job_args_list
        @async remotecall_wait(Main._do_work, wp, job_id, job_f, job_args, results)
    end

    # Wait and Collect results
    time_taken = reduce((x,y)-> x+y,
             let tt = v
               tt
             end
             for (_,v) in current_times(job_id)
    )
    println(time_taken, " seconds")
    current_results(job_id)

end

end  # module

max_stars_repo_name: yongshanding/VLQ

max_issues_repo_name: yongshanding/VLQ

max_forks_repo_name: yongshanding/VLQ