text:
@testset "mssql" begin

    # Microsoft drivers (mssql) don't support reading usernames and passwords from .odbc.ini
    # files so you have to specify those explicitly in the DSN setup call:
    # dsn = ODBC.DSN("MSSQL-test", "yourusername", "yourpassword")
    dsn = ODBC.DSN("MSSQL-test", "SA", "YourStrong!Passw0rd")

    @testset "basic queries" begin
        dbs = ODBC.query(dsn, "select name from sys.databases")
        data = ODBC.query(dsn, "select table_name from master.information_schema.tables")
        data = ODBC.query(dsn, "select table_name from tempdb.information_schema.tables")
        data = ODBC.query(dsn, "select table_name from model.information_schema.tables")
        data = ODBC.query(dsn, "select table_name from msdb.information_schema.tables")
    end

    @testset "create testdb" begin
        ODBC.execute!(dsn, "drop database if exists testdb")
        ODBC.execute!(dsn, "create database testdb")
        ODBC.execute!(dsn, "use testdb")
    end

    @testset "test1" begin
        ODBC.execute!(dsn, "drop table if exists test1")
        ODBC.execute!(dsn, "create table test1
                            (test_bigint bigint,
                             test_bit bit,
                             test_decimal decimal,
                             test_int int,
                             test_money money,
                             test_numeric numeric,
                             test_smallint smallint,
                             test_smallmoney smallmoney,
                             test_tiny_int tinyint,
                             test_float float,
                             test_real real,
                             test_date date,
                             test_datetime2 datetime2,
                             test_datetime datetime,
                             test_datetimeoffset datetimeoffset,
                             test_smalldatetime smalldatetime,
                             test_time time,
                             test_char char(1),
                             test_varchar varchar(16),
                             test_nchar nchar(1),
                             test_nvarchar nvarchar(16),
                             test_binary binary(2),
                             test_varbinary varbinary(16)
                            )")
        data = DataFrame(ODBC.Query(dsn, "select * from information_schema.columns where table_name = 'test1'"))
        ODBC.execute!(dsn, "insert test1 VALUES
                            (1, -- bigint
                             1, -- bit
                             1.0, -- decimal
                             1, -- int
                             1.0, -- money
                             1.0, -- numeric
                             1, -- smallint
                             1.0, -- smallmoney
                             1, -- tinyint
                             1.2, -- float
                             1.2, -- real
                             '2016-01-01', -- date
                             '2016-01-01 01:01:01', -- datetime2
                             '2016-01-01 01:01:01', -- datetime
                             '2016-01-01 01:01:01-05:00', -- datetimeoffset
                             '2016-01-01 01:01:01', -- smalldatetime
                             '01:01:01', -- time
                             'A', -- char(1)
                             'hey there sailor', -- varchar
                             'B', -- nchar(1)
                             'hey there sally', -- nvarchar
                             cast(123456 as binary(2)), -- binary
                             cast(123456 as varbinary(16)) -- varbinary
                            )")
        data = DataFrame(ODBC.Query(dsn, "select * from test1"))

        @test size(data) == (1,23)
        @test Tables.schema(data).types == (
            Union{Int64, Missing},
            Union{Int8, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int32, Missing},
            Union{DecFP.Dec64, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int16, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int8, Missing},
            Union{Float64, Missing},
            Union{Float32, Missing},
            Union{ODBC.API.SQLDate, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTime, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{Array{UInt8,1}, Missing},
        )

        @test data[1][1] === Int64(1)
        @test data[2][1] === Int8(1)
        @test data[3][1] === DecFP.Dec64(1)
        @test data[4][1] === Int32(1)
        @test data[5][1] == DecFP.Dec64(1)
        @test data[6][1] === DecFP.Dec64(1)
        @test data[7][1] === Int16(1)
        @test data[8][1] == DecFP.Dec64(1)
        @test data[9][1] === Int8(1)
        @test data[10][1] === Float64(1.2)
        @test data[11][1] === Float32(1.2)
        @test data[12][1] === ODBC.API.SQLDate(2016,1,1)
        @test data[13][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0)
        @test data[14][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0)
        @test_broken data[15][1] === ODBC.API.SQLTimestamp(2016,1,1,6,1,1,0)
        @test data[16][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,0,0)
        @test data[17][1] === ODBC.API.SQLTime(1,1,1)
        @test data[18][1] == "A"
        @test data[19][1] == "hey there sailor"
        @test data[20][1] == "B"
        @test data[21][1] == "hey there sally"
        @test data[22][1] == UInt8[0xe2, 0x40]
        @test data[23][1] == UInt8[0x00, 0x01, 0xe2, 0x40]

        ODBC.execute!(dsn, "insert test1 VALUES
                            (2, -- bigint
                             1, -- bit
                             2.0, -- decimal
                             2, -- int
                             2.0, -- money
                             2.0, -- numeric
                             2, -- smallint
                             2.0, -- smallmoney
                             2, -- tinyint
                             2.2, -- float
                             2.2, -- real
                             '2016-01-01', -- date
                             '2016-01-01 01:01:01', -- datetime2
                             '2016-01-01 01:01:01', -- datetime
                             '2016-01-01 01:01:01-05:00', -- datetimeoffset
                             '2016-01-01 01:01:01', -- smalldatetime
                             '01:01:01', -- time
                             'A', -- char(1)
                             'hey there sailor', -- varchar
                             'B', -- nchar(1)
                             'hey there sally', -- nvarchar
                             cast(123456 as binary(2)), -- binary
                             cast(123456 as varbinary(16)) -- varbinary
                            )")
        data = DataFrame(ODBC.Query(dsn, "select * from test1"))
        @test size(data) == (2,23)
        @test data[1][1] === Int64(1)
        @test data[1][2] === Int64(2)

        ODBC.execute!(dsn, "drop table if exists test1")
    end

    @testset "employee" begin
        ODBC.execute!(dsn, "drop table if exists employee")

        ODBC.execute!(dsn, """
        CREATE TABLE Employee
        (
            ID INT NOT NULL IDENTITY(1,1) PRIMARY KEY,
            Name VARCHAR(255),
            Salary FLOAT,
            JoinDate DATE,
            LastLogin DATETIME,
            LunchTime TIME,
            OfficeNo TINYINT,
            Senior BIT,
            empno SMALLINT
        );""")
        ODBC.execute!(dsn, """
        INSERT INTO Employee
        (Name, Salary, JoinDate, LastLogin, LunchTime, OfficeNo, Senior, empno)
         VALUES ('John', 10000.50, '2015-8-3', '2015-9-5 12:31:30', '12:00:00', 1, 1, 1301),
         ('Tom', 20000.25, '2015-8-4', '2015-10-12 13:12:14', '13:00:00', 12, 1, 1422),
         ('Jim', 30000.00, '2015-6-2', '2015-9-5 10:05:10', '12:30:00', 45, 0, 1567),
         ('Tim', 15000.50, '2015-7-25', '2015-10-10 12:12:25', '12:30:00', 56, 1, 3200);
        """)
        data = DataFrame(ODBC.Query(dsn, "select * from employee"))
        @test size(data) == (4,9)
        ODBC.execute!(dsn, "drop table if exists employee")
    end

    ODBC.disconnect!(dsn)
end


text:
#=
surface_test:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-19
=#
# using Revise
using Test
using Logging
using SparseArrays
using Distributed
if nprocs() == 1
    addprocs(3)
end

@everywhere using LarSurf
# using Plasm
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation

@testset "Extract surface grid" begin
    segmentation = zeros(Int8, 5, 6, 7)

    segmentation[2:5,2:5,2:6] .= 1
    obj_sz = [4, 4, 5]
    # Plasm.view(Plasm.numbering(.6)((V,[VV, EV, filteredFV])))

    # reducedLARmodel, Flin, (V, model) = LarSurf.get_surface_grid_old(segmentation; return_all=true)
    reducedLARmodel, Flin, (V, model) = LarSurf.get_surface_grid(segmentation; return_all=true)

    (VV, EV, FV, CV) = model
    filteredFV = reducedLARmodel[2][1]
    larmodel1 = (V,[VV, EV, filteredFV])

    @test LarSurf.check_LARmodel(larmodel1)
    # using Plasm
    # Plasm.View(larmodel1)
end

@testset "Extract surface grid and check euler" begin

    segmentation = LarSurf.data234()
    # reducedLARmodel, Flin, (V, model) = LarSurf.get_surface_grid_old(segmentation; return_all=true)
    reducedLARmodel, Flin, (V, model) = LarSurf.get_surface_grid(segmentation; return_all=true)

    (VV, EV, FV, CV) = model
    filteredFV = reducedLARmodel[2][1]
    larmodel1 = (V,[VV, EV, filteredFV])

    # Plasm.View(larmodel1)
    display(filteredFV)

    @test LarSurf.check_surface_euler(filteredFV)

end

@testset "Extract surface grid per block full" begin

    # segmentation = zeros(Int8, 2, 3, 4)
    # segmentation[1:2,2:3,3:4] .= 1
    segmentation = LarSurf.data234()
    # Plasm.view(Plasm.numbering(.6)((V,[VV, EV, filteredFV])))
    # just for visualization
    data_size = LarSurf.size_as_array(size(segmentation))
    V, (VV, EV, FV, CV) = LarSurf.Lar.cuboidGrid(data_size, true)
    Fchar = LarSurf.__grid_get_surface_Fchar_per_block(segmentation, [2,2,2])

    V1, topology1 = LarSurf.get_surface_grid_per_block_full(segmentation, [2,2,2])
    # (VV, EV, FV, CV) = model
    larmodel1 = (V,[VV, EV, topology1[1]])
    @test LarSurf.check_LARmodel(larmodel1)
    @test LarSurf.check_surface_euler(topology1[1])
    # Plasm.View(larmodel1)

    # second implementation
    V2, FV2 = LarSurf.grid_Fchar_to_V_FVreduced(Fchar, size(segmentation))
    larmodel2 = (V2,[VV, EV, FV2])
    @test LarSurf.check_LARmodel(larmodel2)
    @test LarSurf.check_surface_euler(FV2)
    # Plasm.view( Plasm.numbering(.6)(larmodel2) )

    # third implementation
    V3, FV3 = LarSurf.grid_Fchar_to_Vreduced_FVreduced_old(Fchar, data_size)
    larmodel3 = (V3, [FV3])
    @test LarSurf.check_LARmodel(larmodel3)
    @test LarSurf.check_surface_euler(FV3)
    # Plasm.view(larmodel3)

    V4, FV4 = LarSurf.grid_Fchar_to_Vreduced_FVreduced(Fchar, data_size)
    larmodel4 = (V4, [FV4])
    @test LarSurf.check_LARmodel(larmodel4)
    @test LarSurf.check_surface_euler(FV4)
    # VV3 = [[i] for i=1:size(V3,2)]
    # Plasm.view((V3, [VV3, FV3]))
    # Plasm.view(larmodel4)

end

@testset "Extract surface grid per block reduced FV" begin

    segmentation = LarSurf.data234()
    # segmentation = zeros(Int8, 2, 3, 4)
    # segmentation[1:2,2:3,3:4] .= 1
    # obj_sz = [2, 2, 2]
    # Plasm.view(Plasm.numbering(.6)((V,[VV, EV, filteredFV])))

    larmodel1 = LarSurf.get_surface_grid_per_block_FVreduced(segmentation, [2,2,2])
    @test LarSurf.check_LARmodel(larmodel1)
    @test LarSurf.check_surface_euler(larmodel1[2][1])
    # Plasm.View(larmodel1)

end

@testset "Extract surface grid per block reduced FV and V" begin

    segmentation = zeros(Int8, 5, 6, 7)
    segmentation[2:5,2:5,2:6] .= 1
    obj_sz = [4, 4, 5]

    segmentation = zeros(Int8, 2, 3, 4)

    segmentation[1:2,2:3,3:4] .= 1
    obj_sz = [2, 2, 2]
    # Plasm.view(Plasm.numbering(.6)((V,[VV, EV, filteredFV])))

    larmodel1 = LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_old(segmentation, [2,2,2])
    @test LarSurf.check_LARmodel(larmodel1)
    @test LarSurf.check_surface_euler(larmodel1[2][1])
    # Plasm.View(larmodel1)

    larmodel1 = LarSurf.get_surface_grid_per_block_Vreduced_FVreduced(segmentation, [2,2,2])
    @test LarSurf.check_LARmodel(larmodel1)
    @test LarSurf.check_surface_euler(larmodel1[2][1])
    # Plasm.View(larmodel1)

end

@testset "Get Flin test with fixed block" begin

    segmentation = LarSurf.data234()
    # segmentation = zeros(Int8, 2, 3, 4)
    # segmentation[1:2,2:3,3:4] .= 1
    # obj_sz = [2, 2, 2]
    block_size = [2,2,2]

    Flin = LarSurf.__grid_get_surface_Fchar_per_block(segmentation, block_size)
    # Slin, oneS, b3 = LarSurf.grid_get_surf_Fvec_larmodel_loc_fixed_block_size(segmentation, [2,2,2])
    Flin_loc, offsets, blocks_number_axis, larmodel1 = LarSurf.grid_get_surface_Bchar_loc_fixed_block_size(segmentation, block_size)
    Flin, new_data_size = LarSurf.__grid_get_surface_Fchar_per_fixed_block_size(segmentation, block_size)
    larmodel1 = LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_fixed_block_size(segmentation, block_size)
    @test LarSurf.check_LARmodel(larmodel1)
    @test LarSurf.check_surface_euler(larmodel1[2][1])
    # Plasm.View(larmodel1)
end

@testset "Get Flin test parallel" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]

    # Flin = LarSurf.__grid_get_surface_Fchar_per_block(segmentation, block_size)
    # Slin, oneS, b3 = LarSurf.grid_get_surf_Fvec_larmodel_loc_fixed_block_size(segmentation, [2,2,2])
    # Flin_loc, offsets, blocks_number_axis, larmodel1 = LarSurf.grid_get_surface_Bchar_loc_fixed_block_size(segmentation, block_size)
    # Flin, new_data_size = LarSurf.__grid_get_surface_Fchar_per_fixed_block_size(segmentation, block_size)
    larmodel1 = LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_parallel(segmentation, block_size)
    @test LarSurf.check_LARmodel(larmodel1)
    @test LarSurf.check_surface_euler(larmodel1[2][1])
    # Plasm.View(larmodel1)
end

@testset "Get Flin test in one block" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]
    number_of_blocks, bgetter = LarSurf.block_getter(segmentation, block_size)
    LarSurf.get_block(2, bgetter...)
    @test number_of_blocks ==4
    fids = LarSurf.__grid_get_surface_get_Fids_used_in_block(3, bgetter...)
    @test length(fids) > 0
    # TODO why is there only 16 faces? There should be 20 of them.
    # @test length(fids) == 20
    # segmentation = zeros(Int8, 2, 3, 4)
    # segmentation[1:2,2:3,3:4] .= 1
    # obj_sz = [2, 2, 2]
end

@testset "Get Fchar test" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]

    # Flin1 = LarSurf.__grid_get_surface_Fchar(segmentation, block_size)
    # println("get fchar test 1")
    Flin0, lmodel0 = LarSurf.grid_get_surf_Fvec_larmodel(segmentation)
    @test nnz(Flin0) == 6*4
    # println("get fchar test 2")
    Flin1 = LarSurf.__grid_get_surface_Fchar_per_block(segmentation, block_size)
    @test nnz(Flin1) == 6*4

    # println("get fchar test 3")
    Flin2 = LarSurf.__grid_get_surface_Fchar_per_block_parallel_pmap(segmentation, block_size)
    @test nnz(Flin2) == 6*4
    # println("get fchar test 4")
    # TODO concider to make channel implementation work
    # Flin2 = LarSurf.__grid_get_surface_Fchar_per_block_parallel_channel(segmentation, block_size)
    # @test nnz(Flin2) == 6*4
    # println("get fchar test 5")

    # Slin, oneS, b3 = LarSurf.grid_get_surf_Fvec_larmodel_loc_fixed_block_size(segmentation, [2,2,2])
    # Flin_loc, offsets, blocks_number_axis, larmodel1 = LarSurf.grid_get_surface_Bchar_loc_fixed_block_size(segmentation, block_size)
    # Flin, new_data_size = LarSurf.__grid_get_surface_Fchar_per_fixed_block_size(segmentation, block_size)
    # larmodel1 = LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_fixed_block_size(segmentation, block_size)
    # @test LarSurf.check_LARmodel(larmodel1)
    # @test LarSurf.check_surface_euler(larmodel1[2][1])
    # Plasm.View(larmodel1)
end

@testset "surface extraction all functions" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]

    # whole at once functions
    fns_one = [
        LarSurf.get_surface_grid,
        LarSurf.get_surface_grid_old,
    ]
    for fun in fns_one
        larmodel1 = fun(segmentation)
        @test LarSurf.check_LARmodel(larmodel1)
        @test LarSurf.check_surface_euler(larmodel1[2][1])
    end

    # per block functions

    fns = [
        LarSurf.get_surface_grid_per_block,
        LarSurf.get_surface_grid_per_block_full,
        LarSurf.get_surface_grid_per_block_FVreduced,
        LarSurf.get_surface_grid_per_block_Vreduced_FVreduced,
        LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_fixed_block_size,
        LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_old,
        LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_parallel
        # LarSurf.get_surface_grid_per_block,
        # LarSurf.get_surface_grid_per_bloc,
    ]
    for fun in fns
        larmodel1 = fun(segmentation, block_size)
        @test LarSurf.check_LARmodel(larmodel1)
        @test LarSurf.check_surface_euler(larmodel1[2][1])
    end
end

@testset "Get Fchar parallel test" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]
    Flin2 = LarSurf.__grid_get_surface_Fchar_per_block_parallel_pmap(segmentation, block_size)
    println(Flin2)
    @test nnz(Flin2) == 6*4
end

ext: jl


text:
module PointsOnASphere

export Point3D,Point2D,SphericalPoint

struct Point3D{R<:Real,T<:Real,P<:Real}
	r :: R
	θ :: T
	ϕ :: P

	function Point3D(r::R,θ::T,ϕ::P) where {R,T,P}
		r >= 0 || throw(ArgumentError("r must be non-negative"))
		0 <= θ <= π || throw(ArgumentError("θ should lie in [0,π]"))
		0 <= ϕ < 2π || throw(ArgumentError("ϕ should lie in [0,2π)"))
		new{R,T,P}(r,θ,ϕ)
	end
end

struct Point2D{T<:Real,P<:Real}
	θ :: T
	ϕ :: P

	function Point2D(θ::T,ϕ::P) where {T,P}
		0 <= θ <= π || throw(ArgumentError("θ should lie in [0,π]"))
		0 <= ϕ < 2π || throw(ArgumentError("ϕ should lie in [0,2π)"))
		new{T,P}(θ,ϕ)
	end
end

Point3D(r::Real,n::Point2D) = Point3D(r,n.θ,n.ϕ)
Point3D(x::Point3D) = x
Point3D(n::Point2D) = Point3D(1,n)
Point3D(n::Tuple{Real,Real}) = Point3D(1,n...)
Point3D(r::Real,n::Tuple{Real,Real}) = Point3D(r,n...)

Point2D(x::Point3D) = Point2D(x.θ,x.ϕ)
Point2D(x::Point2D) = x

Point2D(n::Tuple{Real,Real}) = Point2D(n...)
Point3D(x::Tuple{Real,Real,Real}) = Point3D(x...)

const SphericalPoint = Union{Point2D,Point3D}

# display
Base.show(io::IO, pt::Point3D) = print(io,"(r=$(round(pt.r,sigdigits=3)), θ=$(round(pt.θ,sigdigits=3)), ϕ=$(round(pt.ϕ,sigdigits=3)))")
Base.show(io::IO, pt::Point2D) = print(io,"(θ=$(round(pt.θ,sigdigits=3)), ϕ=$(round(pt.ϕ,sigdigits=3)))")

Base.show(io::IO, ::MIME"text/plain", pt::Point3D) = print(io,"(r=$(pt.r), θ=$(pt.θ), ϕ=$(pt.ϕ))")
Base.show(io::IO, ::MIME"text/plain", pt::Point2D) = print(io,"(θ=$(pt.θ), ϕ=$(pt.ϕ))")

end # module


text:
using Cbc

include(joinpath(Pkg.dir("MathProgBase"),"test","mixintprog.jl"))
include(joinpath(dirname(@__FILE__),"..","examples","knapsack.jl"))

mixintprogtest(CbcSolver(logLevel=1))

include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
solver = CbcSolver(logLevel=1)
MathProgBase.setparameters!(solver, Silent=true, TimeLimit=100.0)
println("Cbc should not display any output from these tests")
coniclineartest(solver)

if isdir(Pkg.dir("JuMP"))
    include(joinpath(Pkg.dir("JuMP"),"test","runtests.jl"))
end


text:
load("utils/req.jl")
req("utils/failexpect.jl")
req("dag/pshow_dag.jl")
req("tangle.jl")
req("midsection.jl")
req("julia_backend.jl")

code = :( let
    dest[...] = scatter(A.*B + C)
    dest2[...] = A[scatter(B)]
end )

rawdag = tangle(code.args[1])
println("Raw DAG:")
pprintln(rawdag)

sdag = scatter_propagated(rawdag)
println("\nScattered DAG:")
pprintln(sdag)

edag = expand_ellipsis_indexing(sdag, [:i,:j])
println("\n[...] expanded:")
pprintln(edag)

println("\nuntangled:")
print_list(untangle(edag)[2])

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
if VERSION>=v"1.4" && Sys.isapple() && !(haskey(ENV, "DOCUMENTER_KEY"))
    error("""Your Julia version is ≥1.4, and your operation system is MacOSX.
Currently, there is a compatibility issue for this combination.
Please downgrade your Julia version.""")
end

if haskey(ENV, "MANUAL") && ENV["MANUAL"]=="1"
    error("""****** You indicated you want to build ADCME package manually.
To this end, you need to create a dependency file
$(joinpath(@__DIR__, "deps.jl"))
and populate it with appropriate binary locations.
--------------------------------------------------------------------------------------------
BINDIR = ""
LIBDIR = ""
TF_INC = ""
TF_ABI = ""
EIGEN_INC = ""
CC = ""
CXX = ""
CMAKE = ""
MAKE = ""
GIT = ""
PYTHON = ""
TF_LIB_FILE = ""
LIBCUDA = ""
CUDA_INC = ""
__STR__ = join([BINDIR,LIBDIR,TF_INC,TF_ABI,EIGEN_INC,CC,CXX,CMAKE,MAKE,GIT,PYTHON,TF_LIB_FILE,LIBCUDA,CUDA_INC], ";")
--------------------------------------------------------------------------------------------
""")
end

LIBCUDA = ""
CUDA_INC = ""
if isfile(joinpath(@__DIR__, "deps.jl"))
    include(joinpath(@__DIR__, "deps.jl"))
end

if length(LIBCUDA)>0
    @warn("We detect once you have installed ADCME-GPU, we will keep using the corresponding CUDA and CUDA library path as default.")
end

push!(LOAD_PATH, "@stdlib")
using Pkg
using Conda
using CMake
using LibGit2

ENVDIR = abspath(Conda.ROOTENV)

VER = haskey(Pkg.installed(),"ADCME")  ? Pkg.installed()["ADCME"] : "NOT_INSTALLED"
@info """Your Julia version is $VERSION, current ADCME version is $VER, ADCME dependencies installation path: $ENVDIR"""

@info " --------------- (1/6) Install Tensorflow Dependencies  --------------- "

if haskey(ENV, "FORCE_REINSTALL_ADCME") && ENV["FORCE_REINSTALL_ADCME"] in [1, "1"]
    @info " --------------- Install ADCME by force  --------------- "
    include("install_conda.jl")
end

if !("adcme" in Conda._installed_packages())
    try
        Conda.add("adcme", channel="kailaix")
    catch
        @info " --------------- Retry with a fresh new environment  --------------- "
        include("install_conda.jl")
        Conda.add("adcme", channel="kailaix")
    end
end

BINDIR = Sys.iswindows() ? abspath("$ENVDIR/Scripts") : abspath("$ENVDIR/bin")

GIT = "LibGit2"
PYTHON = joinpath(BINDIR, "python")

if Sys.iswindows()
    PYTHON = abspath(joinpath(ENVDIR, "python.exe"))
end
@info " --------------- (2/6) Check Python Version  --------------- "

!haskey(Pkg.installed(), "PyCall") && Pkg.add("PyCall")
ENV["PYTHON"]=PYTHON
Pkg.build("PyCall")
using PyCall
@info """
PyCall Python version: $(PyCall.python)
Conda Python version: $PYTHON
"""

@info " --------------- (3/6) Looking for TensorFlow Dynamic Libraries --------------- "
tf = pyimport("tensorflow")
core_path = abspath(joinpath(tf.sysconfig.get_compile_flags()[1][3:end], ".."))
lib = readdir(core_path)
if Sys.iswindows()
    global TF_LIB_FILE = abspath(joinpath(core_path, "python/_pywrap_tensorflow_internal.lib"))
else
    global TF_LIB_FILE = joinpath(core_path,lib[findall(occursin.("libtensorflow_framework", lib))[end]])
end
TF_INC = tf.sysconfig.get_compile_flags()[1][3:end]
TF_ABI = tf.sysconfig.get_compile_flags()[2][end:end]

if Sys.iswindows()
    if !isdir(joinpath(TF_INC, "tensorflow"))
        @info " --------------- (Windows) Downloading Include Files for Custom Operators --------------- "
        run(`cmd /c rmdir /s /q $TF_INC`)
        LibGit2.clone("https://github.com/kailaix/tensorflow-1.15-include", TF_INC)
    end
end

@info " --------------- (4/6) Preparing Custom Operator Environment --------------- "
LIBDIR = abspath("$ENVDIR/lib/Libraries")

if !isdir(LIBDIR)
    @info "Downloading dependencies to $LIBDIR..."
    mkdir(LIBDIR)
end

if !isfile("$LIBDIR/eigen.zip")
    download("http://bitbucket.org/eigen/eigen/get/3.3.7.zip","$LIBDIR/eigen.zip")
end

if !isdir("$LIBDIR/eigen3")
    UNZIP =  joinpath(BINDIR, "unzip")
    if Sys.iswindows()
        if !isfile("$LIBDIR/unzip.exe")
            download("http://stahlworks.com/dev/unzip.exe", joinpath(LIBDIR, "unzip.exe"))
        end
        UNZIP =  joinpath(LIBDIR, "unzip.exe")
    end
    run(`$UNZIP -qq $LIBDIR/eigen.zip -d $LIBDIR`)
    mv("$LIBDIR/eigen-eigen-323c052e1731", "$LIBDIR/eigen3", force=true)
end

# If the system has `nvcc` but "GPU" is not specified, warn the users to build with
# ENV["GPU"] = 1
if !haskey(ENV, "GPU")
    try
        if Sys.iswindows()
            run(`cmd /c nvcc --version`)
        else
            run(`which nvcc`)
        end
        @warn("""We detected that you have `nvcc` installed but ENV[\"GPU\"] is not set.
>>> If you want to install ADCME with GPU capabiity enabled, please set `ENV[\"GPU\"]=1`.""")
    catch
    end
end

if haskey(ENV, "GPU") && ENV["GPU"]=="1" && !(Sys.isapple())
    @info " --------------- (5/6) Installing GPU Dependencies --------------- "
    try
        if Sys.iswindows()
            run(`cmd /c nvcc --version`)
        else
            run(`which nvcc`)
        end
    catch
        error("""You specified ENV["GPU"]=1 but nvcc cannot be found.
Make sure `nvcc` is available.""")
    end
    s = Sys.iswindows() ? join(readlines(pipeline(`cmd /c nvcc --version`)), " ") : join(readlines(pipeline(`nvcc --version`)), " ")
    ver = match(r"V(\d+\.\d)", s)[1]
    if ver[1:2]!="10"
        error("TensorFlow backend of ADCME requires CUDA 10.0. But you have CUDA $ver")
    end
    if ver[1:4]!="10.0"
        @warn("TensorFlow is compiled using CUDA 10.0, but you have CUDA $ver. This might cause some problems.")
    end

    if !("adcme-gpu" in Conda._installed_packages())
        Conda.add("adcme-gpu", channel="kailaix")
    end

    pkg_dir = joinpath(Conda.ROOTENV, "pkgs/")
    files = readdir(pkg_dir)
    libpath = filter(x->startswith(x, "cudatoolkit") && isdir(joinpath(pkg_dir,x)), files)
    if length(libpath)==0
        @warn "cudatoolkit* not found in $pkg_dir"
    elseif length(libpath)>1
        @warn "more than 1 cudatoolkit found, use $(libpath[1]) by default"
    end

    if length(libpath)>=1
        LIBCUDA = abspath(joinpath(pkg_dir, libpath[1], "lib"))
    end

    libpath = filter(x->startswith(x, "cudnn") && isdir(joinpath(pkg_dir,x)), files)
    if length(libpath)==0
        @warn "cudnn* not found in $pkg_dir"
    elseif length(libpath)>1
        @warn "more than 1 cudatoolkit found, use $(libpath[1]) by default"
    end

    if length(libpath)>=1
        LIBCUDA = LIBCUDA* (Sys.iswindows() ? ";" : ":") *abspath(joinpath(pkg_dir, libpath[1], "lib"))
        @info " --------------- CUDA include headers  --------------- "
        cudnn = Sys.iswindows() ?
            joinpath(pkg_dir, libpath[1], "Library/include", "cudnn.h") :
            joinpath(pkg_dir, libpath[1], "include", "cudnn.h")
        cp(cudnn, joinpath(TF_INC, "cudnn.h"), force=true)
    end

    NVCC = Sys.iswindows() ?
            strip(String(read(`cmd /c where nvcc`))) :
             readlines(pipeline(`which nvcc`))[1]
    CUDA_INC = joinpath(splitdir(splitdir(NVCC)[1])[1], "include")

else
    @info " --------------- (5/6) Skipped: Installing GPU Dependencies  --------------- "
end

@info """ -------------- (6/6) Write Dependency Files  --------------- """

s = ""
t = []
function adding(k, v)
    global s
    if Sys.iswindows()
        v = replace(v, "\\"=>"\\\\")
    end
    s *= "$k = \"$v\"\n"
    push!(t, "$k")
end
adding("BINDIR", BINDIR)
if Sys.iswindows()
    D = abspath(joinpath(LIBDIR, "lib"))
    !isdir(D) && mkdir(D)
    adding("LIBDIR", D)
else
    adding("LIBDIR", abspath(joinpath(ENVDIR, "lib")))
end
adding("TF_INC", TF_INC)
adding("TF_ABI", TF_ABI)
adding("PREFIXDIR", LIBDIR)
if Sys.isapple()
    adding("CC", joinpath(BINDIR, "clang"))
    adding("CXX", joinpath(BINDIR, "clang++"))
elseif Sys.islinux()
    adding("CC", joinpath(BINDIR, "x86_64-conda_cos6-linux-gnu-gcc"))
    adding("CXX", joinpath(BINDIR, "x86_64-conda_cos6-linux-gnu-g++"))
else
    adding("CC", "")
    adding("CXX", "")
end
adding("CMAKE", cmake)
if Sys.iswindows()
    adding("MAKE", "C:\\Program Files (x86)\\Microsoft Visual Studio\\2017\\Community\\MSBuild\\15.0\\Bin\\MSBuild")
else
    adding("MAKE", joinpath(BINDIR, "make"))
end
adding("GIT", GIT)
adding("PYTHON", PyCall.python)
adding("TF_LIB_FILE", TF_LIB_FILE)
adding("LIBCUDA", LIBCUDA)
adding("CUDA_INC", CUDA_INC)
if Sys.iswindows()
    adding("NINJA", "")
else
    adding("NINJA", joinpath(BINDIR, "ninja"))
end
adding("INCDIR", abspath(joinpath(BINDIR, "..", "include")))

t = "join(["*join(t, ",")*"], \";\")"
s *= "__STR__ = $t"
open("deps.jl", "w") do io
    write(io, s)
end

@info """ --------------- Finished: $(abspath("deps.jl"))  --------------- """


text:
Base.Broadcast.promote_containertype(::Type{DataValueArray}, ::Type{DataValueArray}) = DataValueArray
Base.Broadcast.promote_containertype(::Type{Array}, ::Type{DataValueArray}) = DataValueArray
Base.Broadcast.promote_containertype(::Type{DataValueArray}, ::Type{Array}) = DataValueArray
Base.Broadcast.promote_containertype(::Type{DataValueArray}, _) = DataValueArray
Base.Broadcast.promote_containertype(_, ::Type{DataValueArray}) = DataValueArray

Base.Broadcast._containertype(::Type{<:DataValueArray}) = DataValueArray

Base.Broadcast.broadcast_indices(::Type{DataValueArray}, A) = indices(A)

@inline function Base.Broadcast.broadcast_t(f, ::Type{DataValue{T}}, shape, iter, A, Bs::Vararg{Any,N}) where {N,T}
    C = similar(DataValueArray{T}, shape)
    keeps, Idefaults = Base.Broadcast.map_newindexer(shape, A, Bs)
    Base.Broadcast._broadcast!(f, C, keeps, Idefaults, A, Bs, Val{N}, iter)
    return C
end

# broadcast methods that dispatch on the type of the final container
@inline function Base.Broadcast.broadcast_c(f, ::Type{DataValueArray}, A, Bs...)
    T = Base.Broadcast._broadcast_eltype(f, A, Bs...)
    shape = Base.Broadcast.broadcast_indices(A, Bs...)
    iter = CartesianRange(shape)
    if isleaftype(T)
        return Base.Broadcast.broadcast_t(f, T, shape, iter, A, Bs...)
    end
    if isempty(iter)
        return similar(Array{T}, shape)
    end
    return Base.Broadcast.broadcast_t(f, Any, shape, iter, A, Bs...)
end

# This one is much faster than normal broadcasting but the method won't get called
# in fusing operations like (!).(isna.(x))
Base.broadcast(::typeof(isna), data::DataValueArray) = copy(data.isna)


text:
@doc raw"""
    Circle3

An object `c` of type [`Circle3`](@ref) is a circle in the three-dimensional
Euclidean space ``Ε^3``.

Note that the circle can be degenerated, i.e. the squared radius may be zero.
"""
Circle3

@doc raw"""
    Circle3(center::Point3, r²::Real, plane::Plane3)

Introduces a variable `c` of type [`Circle3`](@ref).

It is initalized to the circle of center `center` and squared radius `r²` in
plane `plane`.

!!! info "Precondition"

    `center` lies in `plane` and `r²` ``\ge`` 0.
"""
Circle3(center::Point3, r²::Real, plane::Plane3)

@cxxdereference Circle3(center::Point3, r²::Real, plane::Plane3) =
    Circle3(center, convert(FT, r²), plane)

@doc raw"""
    Circle3(center::Point3, r²::Real, n::Vector3)

Introduces a variable `c` of type [`Circle3`](@ref).

It is initalized to the circle of center `center` and squared radius `r²` in a
plane normal to the vector `n`.

!!! info "Precondition"

    `r²` ``\ge`` 0.
"""
Circle3(center::Point3, r²::Real, n::Vector3)

@cxxdereference Circle3(center::Point3, r²::Real, n::Vector3) =
    Circle3(center, convert(FT, r²), n)

@doc raw"""
    Circle3(p::Point3, q::Point3, r::Point3)

Introduces a variable `c` of type [`Circle3`](@ref).

It is initialized to the circle passing through the three points.

!!! info "Precondition"

    The three points are not collinear.
"""
Circle3(p::Point3, q::Point3, r::Point3)

@doc raw"""
    Circle3(sphere₁::Sphere3, sphere₂::Sphere3)

Introduces a variable `c` of type [`Circle3`](@ref).

It is initialized to the circle along which the two spheres intersect.

!!! info "Precondition"

    The two spheres intersect along a circle.
"""
Circle3(sphere₁::Sphere3, sphere₂::Sphere3)

@doc raw"""
    Circle3(sphere::Sphere3, plane::Plane3)

Introduces a variable `c` of type [`Circle3`](@ref).

It is initialized to the circle along which the sphere and the plane intersect.

!!! info "Precondition"

    The sphere and the plane intersect along a circle.
"""
Circle3(sphere::Sphere3, plane::Plane3)

@doc raw"""
    Circle3(plane::Plane3, sphere::Sphere3)

Introduces a variable `c` of type [`Circle3`](@ref).

It is initialized to the circle along which the sphere and the plane intersect.

!!! info "Precondition"

    The sphere and the plane intersect along a circle.
"""
Circle3(plane::Plane3, sphere::Sphere3)

"""
    center(c::Circle3)

Returns the center of `c`.
"""
center(c::Circle3)

"""
    squared_radius(c::Circle3)

Returns the squared radius of `c`.
"""
squared_radius(c::Circle3)

"""
    supporting_plane(c::Circle3)

Returns the supporting plane of `c`.
"""
supporting_plane(c::Circle3)

"""
    diametral_sphere(c::Circle3)

Returns the diametral sphere of `c`.
"""
diametral_sphere(c::Circle3)

"""
    area_divided_by_pi(c::Circle3)

Returns the area of `c`, divided by ``π``.
"""
area_divided_by_pi(c::Circle3)

"""
    approximate_area(c::Circle3)

Returns an approximation of the area of `c`.
"""
approximate_area(c::Circle3)

"""
    squared_length_divided_by_pi_square(c::Circle3)

Returns the squared length of `c`, divided by ``π^2``.
"""
squared_length_divided_by_pi_square(c::Circle3)

"""
    approximate_squared_length(c::Circle3)

Returns an approximation of the squared length (i.e. perimeter) of `c`.
"""
approximate_squared_length(c::Circle3)

"""
    ==(c::Circle3, circle3::Circle3)

Returns `true`, iff `c` and `circle3` are equal, i.e. if they have the same
center, same squared radius and same supporting plane.
"""
==(c::Circle3, circle3::Circle3)

"""
    bbox(c::Circle3)

Returns a bounding box containing `c`.
"""
bbox(c::Circle3)


text:
## Load Packages
using Omega
using StatsBase
using Random
using Plots
using Distributions

function one_simulation_prey(rng, n, transitions)

    """
    Simulates one step of gillespie for prey.
    """

    #t = t_list[n](rng)
    hazard_result = hazards_list[n](rng)
    prey_val = prey_list[n](rng)
    labels = ["spawn_prey", "pred_dies","prey2pred"]
    transition = transitions[labels[hazard_result]]
    new_prey = prey_val + transition[1]

    # Enforce only positive integers
    max(1, new_prey)

end

function one_simulation_pred(rng, n, transitions)

    """
    Simulates one step of gillespie for pred.
    """

    #t = t_list[n](rng)
    hazard_result = hazards_list[n](rng)
    pred_val = pred_list[n](rng)
    labels = ["spawn_prey", "pred_dies","prey2pred"]
    transition = transitions[labels[hazard_result]]
    new_pred = pred_val + transition[2]

    # Enforce only positive integers
    max(1, new_pred)

end

function get_hazards(rng, n)
    """
    Compute the hazard function given the current states.  Note that in this
        model the function depends only on state, not time, though this is not
        the case in general. "spawn_prey" represents the event of a prey being born,
        "prey2pred" represents a predator consuming a new prey and consequently spawning
        a new predator, "pred_dies" represents the death of a predator.

    args:
        ecology(dict): a dictionary containing species counts.
        theta(dict): A dictionary where keys represent events (reactions) and
            values are the location hyperparameters for rate parameter constants
            corresponding to those events.
    """

    ecology = Dict("prey" => prey_list[n](rng), "pred" => pred_list[n](rng))

    hazards = Dict(
        "spawn_prey" => theta(rng)["spawn_prey"] * ecology["prey"],
        "prey2pred" => theta(rng)["prey2pred"] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta(rng)["pred_dies"] * ecology["pred"]
        )

    vals = collect(values(hazards))
    sum_vals = sum(vals)
    prob_vals = vals/sum_vals
    categorical(rng, prob_vals)
end

function generate_rates(rng)

    Dict("spawn_prey" => spawn_prey(rng),
         "prey2pred" => prey2pred(rng),
         "pred_dies" => pred_dies(rng))

end

## Initialize paramters
Pre = [[1, 0], [1, 1], [0, 1]]
Post = [[2, 0], [0, 2], [0, 0]]
transition_mat = Post - Pre
transitions = Dict("spawn_prey" => transition_mat[1,],
                    "prey2pred" => transition_mat[2,],
                    "pred_dies" => transition_mat[3,])
#t0=0.0

# Initiate starting values randomely
prey_init = normal(10., .001)
pred_init = normal(10., .001)

# Random variables for rates
spawn_prey = normal(1.5, .01)
prey2pred = normal(1.0, .0001)
pred_dies = normal(3.0, .0075)

# prey_init = normal(50., .001)
# pred_init = normal(100., .001)
#
# # Random variables for rates
# spawn_prey = normal(.9, .01)
# prey2pred = normal(.004, .0001)
# pred_dies = normal(.4, .0075)

## Generate Model
hazards_list = Any[]
prey_list = Any[]
pred_list = Any[]
Random.seed!(1234)
theta = ciid(generate_rates)
push!(prey_list, prey_init)
push!(pred_list, pred_init)
N = 1000

for f in 2:N
    last = f - 1
    hazards_temp = ciid(get_hazards, last) # individual step
    prey_temp = ciid(one_simulation_prey, last, transitions) # individual step
    pred_temp = ciid(one_simulation_pred, last, transitions) # individual step
    push!(hazards_list, hazards_temp)
    push!(prey_list, prey_temp)
    push!(pred_list, pred_temp)
end

random_var_tuple = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list)...,
                Tuple(x for x in pred_list)...,
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)

## Sample
Random.seed!(1234)
samples = rand(random_var_tuple,
                1, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,samples[1][N+x])
    push!(pred_vals,samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Omega Simulation",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Predators"],
        lw = 1.25)

a = shuffle(prey_vals)
b = sortperm(a)
condition_index = b[490:499]

## Sample with conditional
hazards_list = Any[]
prey_list = Any[]
pred_list = Any[]
theta = ciid(generate_rates)
push!(prey_list, prey_init)
push!(pred_list, pred_init)
N = 500
for f in 2:N
    last = f - 1
    hazards_temp = ciid(get_hazards, last) # individual step
    if f in 150:151
        prey_temp_ = ciid(one_simulation_prey, last, transitions) # individual step
        prey_temp = cond(prey_temp_, prey_temp_ ==ₛ 5)
    else
        prey_temp = ciid(one_simulation_prey, last, transitions) # individual step
    end
    #if f in 500:750
    #    pred_temp_ = ciid(one_simulation_pred, last, transitions) # individual step
    #    pred_temp = cond(pred_temp_, pred_temp_ < 120)
    #else
    pred_temp = ciid(one_simulation_pred, last, transitions) # individual step
    #end
    push!(hazards_list, hazards_temp)
    push!(prey_list, prey_temp)
    push!(pred_list, pred_temp)
end

random_var_tuple = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list)...,
                Tuple(x for x in pred_list)...,
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)

Random.seed!(1234)
samples = rand(random_var_tuple,
                1, alg = SSMH)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,samples[1][N+x])
    push!(pred_vals,samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Prey <55 at T=2500 Simulation",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Pred"],
        lw = 1.25)

# Histogram to test that conditional worked
compile = []
for x in 1:50
    push!(compile, samples[x][7500])
end
histogram(compile, bins = 10,
        title = "Prey at T=2500 (<55 Conditional)")
## Counterfactual
Random.seed!(1234)
samples = rand(random_var_tuple, 1, alg = RejectionSample)
# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,samples[1][N+x])
    push!(pred_vals,samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Simulation",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Pred"],
        lw = 1.25)

replace_index = 16
new_prey = replace(prey_list[replace_index],
                    prey_list[replace_index - 1] => 10.0)
prey_list_replace = prey_list
prey_list_replace[replace_index] = new_prey

random_var_tuple_int = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list_replace)...,
                Tuple(x for x in pred_list)...,
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)
Random.seed!(1234)
adj_samples = rand(random_var_tuple_int, 1, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,adj_samples[1][N+x])
    push!(pred_vals,adj_samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Intervention Initial Prey=150",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Predators"],
        lw = 1.25)

## Int fix test
replace_index = 16
test = replace(prey_list[replace_index],
                    prey_list[replace_index - 1] => normal(30,1))

Random.seed!(1234)
samples = rand(random_var_tuple, 1, alg = RejectionSample)
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,samples[1][N+x])
    push!(pred_vals,samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Intervention Prey T10=300",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Pred"],
        lw = 1.25)

## Intervene on initial value
prey_init_new = replace(prey_init, prey_init => 150.)
prey_list_replace = prey_list
prey_list_replace[1] = prey_init_new
random_var_tuple_int = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list_replace)...,
                Tuple(x for x in pred_list)...,
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)
Random.seed!(1234)
adj_samples = rand(random_var_tuple_int, 1, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,adj_samples[1][N+x])
    push!(pred_vals,adj_samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Intervention Initial Prey=150",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Predators"],
        lw = 1.25)

## Recreate Plots
hazards_list = Any[]
prey_list = Any[]
pred_list = Any[]
Random.seed!(2)
theta = ciid(generate_rates)
push!(prey_list, prey_init)
push!(pred_list, pred_init)
N = 200

for f in 2:N
    last = f - 1
    hazards_temp = ciid(get_hazards, last) # individual step
    prey_temp = ciid(one_simulation_prey, last, transitions) # individual step
    pred_temp = ciid(one_simulation_pred, last, transitions) # individual step
    push!(hazards_list, hazards_temp)
    push!(prey_list, prey_temp)
    push!(pred_list, pred_temp)
end

random_var_tuple = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list)...,
                Tuple(x for x in pred_list)...,
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)

## A
Random.seed!(22)
samples = rand(random_var_tuple,
                1, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,samples[1][N+x])
    push!(pred_vals,samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Omega Simulation",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Predators"],
        lw = 1.25)

## B
condition_index = [27, 75, 122, 187] #27

hazards_list = Any[]
prey_list = Any[]
pred_list = Any[]
Random.seed!(2)
theta = ciid(generate_rates)
push!(prey_list, prey_init)
push!(pred_list, pred_init)
N = 200
for f in 2:N
    last = f - 1
    hazards_temp = ciid(get_hazards, last) # individual step
    if f in condition_index
        prey_temp_ = ciid(one_simulation_prey, last, transitions) # individual step
        prey_temp = cond(prey_temp_, prey_temp_ >= 6.)
    else
        prey_temp = ciid(one_simulation_prey, last, transitions) # individual step
    end
    pred_temp = ciid(one_simulation_pred, last, transitions) # individual step
    push!(hazards_list, hazards_temp)
    push!(prey_list, prey_temp)
    push!(pred_list, pred_temp)
end

random_var_tuple = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list)...,
                Tuple(x for x in pred_list)...,
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)

Random.seed!(22)
samples = rand(random_var_tuple,
                1, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals, samples[1][N+x])
    push!(pred_vals, samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Condition On Too Many Prey",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Pred"],
        lw = 1.25)

## C/E
replace_index = 15
new_prey = replace(prey_list[replace_index],
                    prey_list[replace_index - 1] => 10.0)
prey_list_replace = prey_list
prey_list_replace[replace_index] = new_prey
random_var_tuple_int = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list_replace)...,
                Tuple(x for x in pred_list)...,#pred_list_replace
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,adj_samples[1][N+x])
    push!(pred_vals,adj_samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Pred=6 at t=15",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Predators"],
        lw = 1.25)

function increase_prey()

    Random.seed!(rand(1:10000))
    samples = rand(random_var_tuple, 1, alg = RejectionSample)
    adj_samples = rand(random_var_tuple_int, 1, alg = RejectionSample)

    # extract run results and plot
    prey_norm_vals = []
    prey_adj_vals = []
    for x in 1:(N-1)
        push!(prey_norm_vals,samples[1][N+x])
        push!(prey_adj_vals,adj_samples[1][(N*2)+x])
    end

    return sum(prey_adj_vals) - sum(prey_norm_vals)
end

check = [increase_prey() for x=1:100]

histogram(check, bins = 15,
        title = "blah")

##D/F
replace_index = 15
new_pred = replace(pred_list[replace_index],
                    pred_list[replace_index - 1] => 12.0)
pred_list_replace = pred_list
pred_list_replace[replace_index] = new_pred
random_var_tuple_int = (Tuple(x for x in hazards_list)...,
                Tuple(x for x in prey_list)...,
                Tuple(x for x in pred_list_replace)...,#pred_list_replace
                Tuple(Any[spawn_prey,prey2pred,pred_dies,theta])...)

Random.seed!(22)
adj_samples = rand(random_var_tuple_int, 1, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:(N-1)
    push!(prey_vals,adj_samples[1][N+x])
    push!(pred_vals,adj_samples[1][(N*2)+x])
end

plot(hcat(prey_vals,pred_vals),
        title = "Pred=6 at t=15",
        xlabel = "Time",
        ylabel = "Quantity",
        label = ["Prey" "Predators"],
        lw = 1.25)

function increase_pred()

    samples = rand(random_var_tuple, 1, alg = RejectionSample)
    adj_samples = rand(random_var_tuple_int, 1, alg = RejectionSample)

    # extract run results and plot
    prey_norm_vals = []
    prey_adj_vals = []
    for x in 1:(N-1)
        push!(prey_norm_vals,samples[1][N+x])
        push!(prey_adj_vals,adj_samples[1][(N*2)+x])
    end

    return sum(prey_adj_vals) - sum(prey_norm_vals)
end

check = [increase_pred() for x=1:10]

histogram(check,
        title = "blah")


text:
using GAP_pkg_datastructures
using Test

import GAP_pkg_datastructures
println("TODO")


text:
using OptionalUnits
using Documenter

DocMeta.setdocmeta!(OptionalUnits, :DocTestSetup, :(using OptionalUnits, Unitful); recursive=true)

makedocs(;
    modules=[OptionalUnits],
    authors="Justin Ackers <justin.ackers@imte.fraunhofer.de> and contributors",
    repo="https://github.com/jusack/OptionalUnits.jl/blob/{commit}{path}#{line}",
    sitename="OptionalUnits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jusack.github.io/OptionalUnits.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Open Issues" => "todo.md"
    ],
)

deploydocs(;
    repo="github.com/jusack/OptionalUnits.jl",
    devbranch = "main",
)


text:
## NPL
"""
    ParametersNPL <: AbstractParametersBinary

# Description

Contains info about the difficulty of an item under the N-parameter logistic model.
"""
mutable struct ParametersNPL <: AbstractParametersBinary
    a::Vector{Float64}
    bounds_a::Vector{Vector{Float64}}
    b::Float64
    bounds_b::Vector{Float64}
    prior::Distributions.MultivariateDistribution
    posterior::Distributions.MultivariateDistribution
    chain::Vector{Vector{Float64}}
    expected_information::Matrix{Float64}
    calibrated::Bool

    ParametersNPL(
        a,
        bounds_a,
        b,
        bounds_b,
        prior,
        posterior,
        chain,
        expected_information,
        calibrated,
    ) = new(
        a,
        bounds_a,
        b,
        bounds_b,
        prior,
        posterior,
        chain,
        expected_information,
        calibrated,
    )

    # Random Initializers

    function ParametersNPL(
        N_variate_dist::Distributions.MultivariateDistribution,
        bounds::Vector{Vector{Float64}},
    )
        N = size(bounds, 1)
        pars = truncate_rand(N_variate_dist, bounds)
        new(
            string.("par_", collect(1:N)),
            pars[2:end],
            bounds[2:end],
            pars[1],
            bounds[1],
            N_variate_dist,
            N_variate_dist,
            [pars[1]],
            LinearAlgebra.I(N),
            true,
        )
    end

end


text:
module problem120

description = """
Square remainders
Let r be the remainder when (a−1)^n + (a+1)^n is divided by a^2.
For example, if a = 7 and n = 3, then r = 42: 63 + 83 = 728 ≡ 42 mod 49. And as n varies, so too will r, but for a = 7 it turns out that rmax = 42.

For 3 ≤ a ≤ 1000, find ∑ rmax.
"""

# compute the maximum value for r = ((a-1)^n + (a+1)^n) % (a^2)
# For increasing values of n, this is a repeating sequence.
# Stop searching for max r when the first value is repeated.
# Perform the power operations using modulo arithmetic to prevent overflow, reusing the results of the previous iterations.
function rmax(a :: Int)
  as1 = a - 1
  ap1 = a + 1
  asq = a * a
  r1 = (as1 + ap1) % asq
  result = r1
  while true
    as1 = (as1 * (a-1)) % asq
    ap1 = (ap1 * (a+1)) % asq
    r = (as1 + ap1) % asq
    if r == r1
      return result
    end
    result = max(result, r)
  end
end

function sumrmax(low=3, high=1000)
  return sum(map(rmax, low:high))
end

using Base.Test
@test rmax(7) == 42

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using CUDA, Plots, Printf
# pyplot()
do_viz = false

function compute_flux!(qx, qy, T, λ, dx, dy, nx, ny)
    ix = (blockIdx().x-1) * blockDim().x + threadIdx().x # thread ID, dimension x
    iy = (blockIdx().y-1) * blockDim().y + threadIdx().y # thread ID, dimension y

    if (ix<=nx-1 && 2<=iy<=ny-1)  qx[ix,iy-1] = -λ*(T[ix+1,iy]-T[ix,iy])/dx  end
    if (2<=ix<=nx-1 && iy<=ny-1)  qy[ix-1,iy] = -λ*(T[ix,iy+1]-T[ix,iy])/dy  end
    return
end

function update_T!(T, qx, qy, dt, ρCp, dx, dy, nx, ny)
    ix = (blockIdx().x-1) * blockDim().x + threadIdx().x # thread ID, dimension x
    iy = (blockIdx().y-1) * blockDim().y + threadIdx().y # thread ID, dimension y

    if (2<=ix<=nx-1 && 2<=iy<=ny-1)  T[ix,iy] = T[ix,iy] - dt/ρCp*((qx[ix,iy-1]-qx[ix-1,iy-1])/dx + (qy[ix-1,iy]-qy[ix-1,iy-1])/dy)  end
    return
end

@views function heat_2D_gpu()
    # Physics
    lx, ly = 10.0, 10.0
    λ      = 1.0
    ρCp    = 1.0
    nt     = 200
    # Numerics
    BLOCKX = 16
    BLOCKY = 16
    GRIDX  = 8
    GRIDY  = 8
    nx, ny = BLOCKX*GRIDX, BLOCKY*GRIDY
    nout   = 10
    # Derived numerics
    dx, dy = lx/nx, ly/ny
    dt     = min(dx^2,dy^2)/ρCp/λ/4.1
    xc     = LinRange(dx/2, lx-dx/2, nx)
    yc     = LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qx     = CUDA.zeros(Float64, nx-1,ny-2)
    qy     = CUDA.zeros(Float64, nx-2,ny-1)
    # Initial condition
    T      = CuArray( exp.(.-(xc.-lx./2.0).^2 .-(yc.-ly./2.0)'.^2) )
    cuthreads = (BLOCKX, BLOCKY, 1)
    cublocks  = (GRIDX,  GRIDY,  1)
    # Time loop
    for it = 1:nt
        if (it==11) global t0 = Base.time() end
        @cuda blocks=cublocks threads=cuthreads compute_flux!(qx, qy, T, λ, dx, dy, nx, ny)
        synchronize()
        @cuda blocks=cublocks threads=cuthreads update_T!(T, qx, qy, dt, ρCp, dx, dy, nx, ny)
        synchronize()
        # Visualise
        if mod(it,nout)==0 && do_viz
            display(heatmap(xc, yc, Array(T)', xlabel="lx", ylabel="ly", title="heat diffusion, it=$it", clims=(0.,1.)))
            # sleep(.01)
        end
    end
    time_s = (Base.time()-t0)
    @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", time_s, round((2/1e9*nx*ny*sizeof(lx))/(time_s/(nt-10)), sigdigits=2))
    return
end

@time heat_2D_gpu()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
# Window functions

# without α
for ( julia_fun, ippf_prefix ) in   [   (   :bartlett,      "ippsWinBartlett"   ),
                                        (   :blackman,      "ippsWinBlackmanStd"),
                                        (   :hamming,       "ippsWinHamming"    ),
                                        (   :hann,          "ippsWinHann"       ) ]

    julia_fun! = symbol(string(julia_fun, '!'))

    for ( T ) in [  :IPP16f,
                    :IPP32f,
                    :IPP64f,
                    :IPP16fc,
                    :IPP32fc,
                    :IPP64fc  ]

        ippf  = string( ippf_prefix, IPPSuffix( T ) )

        @eval begin
            function $(julia_fun!)( y::Array{$T}, x::Array{$T}  )
                n = length( y )
                length(x) == n || throw( DimensionMismatch( "Inconsistent array lengths." ))
                if n > 0
                    @ippscall( $ippf, ( Ptr{$T},    Ptr{$T},    IPPInt  ),
                                        x,          y,          n       )
                end

                return y
            end

            $(julia_fun!)( x::Array{$T} ) = $(julia_fun!)( x, x )
            $(julia_fun)(  x::Array{$T} ) = $(julia_fun!)( similar(x), x )

            function $(julia_fun!)( y::Array{$T}, n::Integer )
                ny = length(y)
                ny == n || throw( DimensionMismatch( "Trying to fill a $ny element buffer with $n elements." ))
                x  = ones( $T, n )
                $(julia_fun!)( y, x )
            end

            $(julia_fun)( T::Type, n::Integer ) = $(julia_fun!)( Array( T, n ), n )
            $(julia_fun)( n::Integer ) = $(julia_fun!)( Array( IPP64f, n ), n )
        end
    end
end

# with α
for ( julia_fun, ippf_prefix ) in   [   (   :blackman,      "ippsWinBlackman"   ),
                                        (   :kaiser,        "ippsWinKaiser"     ) ]

    julia_fun! = symbol(string(julia_fun, '!'))

    for ( T ) in [  :IPP16f,
                    :IPP32f,
                    :IPP64f,
                    :IPP16fc,
                    :IPP32fc,
                    :IPP64fc  ]

        ippf  = string( ippf_prefix, IPPSuffix( T ) )

        @eval begin
            function $(julia_fun!)( y::Array{$T}, x::Array{$T}, α::FloatingPoint  )
                n = length( y )
                length(x) == n || throw( DimensionMismatch( "Inconsistent array lengths." ))
                if n > 0
                    @ippscall( $ippf, ( Ptr{$T},    Ptr{$T},    IPPInt, IPP64f  ),
                                        x,          y,          n,      α       )
                end

                return y
            end

            $(julia_fun!)( x::Array{$T}, α::FloatingPoint ) = $(julia_fun!)( x, x, α )
            $(julia_fun)(  x::Array{$T}, α::FloatingPoint ) = $(julia_fun!)( similar(x), x, α )

            function $(julia_fun!)( y::Array{$T}, n::Integer, α::FloatingPoint )
                ny = length(y)
                ny == n || throw( DimensionMismatch( "Trying to fill a $ny element buffer with $n elements." ))
                x  = ones( $T, n )
                $(julia_fun!)( y, x, α )
            end

            $(julia_fun)( T::Type, n::Integer, α::FloatingPoint ) = $(julia_fun!)( Array( T, n ), n, α )
            $(julia_fun)( n::Integer, α::FloatingPoint )          = $(julia_fun!)( Array( IPP64f, n ), n, α )
        end
    end
end


text:
#!/usr/bin/env julia

using EvaluationOfRLAlgs, Test
println("Starting tests")
println("Finished tests")

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

using NLopt
const solver = MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LD_SLSQP)

optimizer = MOI.instantiate(solver)

const config = MOIT.TestConfig(atol=1e-2, rtol=1e-2, duals=false,
                               optimal_status=MOI.LOCALLY_SOLVED)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "NLopt"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(optimizer, false)
    @test !MOIU.supports_default_copy_to(optimizer, true)
end

@testset "Non-Linear tests" begin
    MOIT.nlptest(optimizer, config)
end

@testset "Testing getters" begin
    MOI.Test.copytest(MOI.instantiate(solver, with_bridge_type=Float64), MOIU.Model{Float64}())
end

@testset "Bounds set twice" begin
    MOI.Test.set_lower_bound_twice(optimizer, Float64)
    MOI.Test.set_upper_bound_twice(optimizer, Float64)
end

MOI.empty!(optimizer)
bridged = MOIB.full_bridge_optimizer(MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), optimizer), Float64)

@testset "Continuous Linear tests" begin
    MOIT.contlineartest(bridged, config, [
        # Infeasibility and unboundedness not detected by NLopt
        "linear8a", "linear8b", "linear8c", "linear12",
        # Terminates with `:ROUNDOFF_LIMITED` due to the bad scaling of input problem
        "linear9",
        # FIXME invalid NLopt arguments: too many equality constraints
        "linear15"])
end

@testset "Continuous Quadratic tests" begin
    MOIT.qptest(bridged, config)
end

MOIB.add_bridge(bridged, MOIB.Constraint.SOCtoNonConvexQuadBridge{Float64})

@testset "Continuous Conic tests" begin
    MOIT.lintest(bridged, config, [
        # Infeasibility and unboundedness not detected by NLopt
        "lin3", "lin4"
    ])
    MOIT.soctest(bridged, config, [
        # Infeasibility and unboundedness not detected by NLopt
        "soc3"
    ])
    MOIT.rsoctest(bridged, config, [
        # Infeasibility and unboundedness not detected by NLopt
        "rotatedsoc2",
        # FIXME invalid NLopt arguments: bounds 4 fail 1 <= 0 <= 1
        "rotatedsoc3",
    ])
    MOIT.geomeantest(bridged, config, [
        # FIXME Status is ok but solution is not, e.g.
        #   Expression: ≈(MOI.get(model, MOI.ObjectiveValue()), 1.0, atol = atol, rtol = rtol)
        #   Evaluated: 0.006271435980170771 ≈ 1.0 (atol=0.01, rtol=0.01)
        "geomean2v", "geomean2f"
    ])
    MOIT.norminftest(bridged, config, [
        # Infeasibility and unboundedness not detected by NLopt
        "norminf2",
    ])
    MOIT.normonetest(bridged, config, [
        # Infeasibility and unboundedness not detected by NLopt
        "normone2",
    ])
end

@testset "Unit" begin
    exclude = [
        # Integer variables not supported
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3",
        "solve_integer_edge_cases",
        # ObjectiveBound not supported.
        "solve_objbound_edge_cases",
        # NumberOfThreads not supported
        "number_threads",
        # Infeasibility and unboundedness not detected by NLopt
        "solve_unbounded_model"
    ]
    MOIT.unittest(bridged, config, exclude)
end


text:
using ChemEquations
using Test

for str in ("compound", "chemequation", "balance")
    @testset "$str.jl" begin
        include("$str.jl")
    end
end


text:
using Lorenz96
using LogFixPoint16s
using PyPlot
using Statistics, StatsBase
using BFloat16s
# using ColorSchemes
# cmap = ColorMap(ColorSchemes.oslo.colors)

LogFixPoint16s.set_nfrac(7)

# parameters
Nens = 100
N = 10
n = 200
Δt = 0.1
F = 8.0

# coordinates
t = collect(0:N)*Δt
nvec = collect(1:n)

# friction parameter space
αs = collect(1:1e-3:2.0)
nαs = length(αs)

# # preallocate
# f16 = fill(0.0,nαs,Nens)
# lfx16 = fill(0.0,nαs,Nens)

# for (iα,α) = enumerate(αs)

#     # control run
#     X0 = L96(;n,Δt,F,α)[:,1000:end]

#     for ie in 1:Nens
#         # pick random start date from control run
#         X = X0[:,rand(1:size(X0)[2])]

#         # experiments
#         XF64 = L96(Float64;X,N,n,Δt,α,F)[:,end]
#         XF16 = L96(BFloat16;X,N,n,Δt,α,F)[:,end]
#         XLFX16 = L96(LogFixPoint16;X,N,n,Δt,α,F)[:,end]

#         f16[iα,ie] = sqrt(mean((XF64-XF16).^2))
#         lfx16[iα,ie] = sqrt(mean((XF64-XLFX16).^2))
#     end
# end

# E = 5.13    # climatological error

# # monte carlo sampling for mean
# m = Nens
# f16m = fill(0.0,nαs,m)
# lfx16m = fill(0.0,nαs,m)

# for i in 1:m
#     # calculate median from a random half of the ensemble members
#     f16m[:,i] = mean(f16[:,rand(1:Nens,Nens÷2)],dims=2)
#     lfx16m[:,i] = mean(lfx16[:,rand(1:Nens,Nens÷2)],dims=2)
# end

# confidence interval from those
f16lo = [percentile(f16m[i,:],5) for i in 1:nαs]
f16hi = [percentile(f16m[i,:],95) for i in 1:nαs]

lfx16lo = [percentile(lfx16m[i,:],5) for i in 1:nαs]
lfx16hi = [percentile(lfx16m[i,:],95) for i in 1:nαs]

fig,ax = subplots(figsize=(10,3))

ax.semilogy(αs,mean(f16,dims=2)/E,label="Float16")
ax.plot(αs,mean(lfx16,dims=2)/E,label="LogFixPoint16")

ax.fill_between(αs,f16lo/E,f16hi/E,alpha=0.3)
ax.fill_between(αs,lfx16lo/E,lfx16hi/E,alpha=0.3)

ax.legend()

ax.set_ylabel("Error after $(N÷10) mtu [%]")
ax.set_xlabel("Friction parameter α")
ax.set_xlim(1.0,2.0)
tight_layout()


text:
@testset "UnitCell tests" begin
    @test typeof(UnitCell((1.0,))) == UnitCell{1}
    @test typeof(UnitCell((1.0,0.0),(0.0,1.0))) == UnitCell{2}
    @test typeof(UnitCell((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))) == UnitCell{3}
    @test typeof(UnitCell((1.0,0.0,0.0,0.0),(0.0,1.0,0.0,0.0),(0.0,0.0,1.0,0.0),(0.0,0.0,0.0,1.0))) == UnitCell{4}

    uc = UnitCell((1.0,0.0),(0.0,1.0))
    b1 = addBasisSite!(uc, (0.0, 0.0))
    b2 = addBasisSite!(uc, (0.5, 0.5))
    @test length(uc.basis) == 2
    @test length(uc.interactionsOnsite) == 2
    @test length(uc.interactionsField) == 2

    @test_throws ErrorException addInteraction!(uc, b1, b1, ones(3,3), (0,0))
    @test_throws ErrorException addInteraction!(uc, b1, b2, ones(4,3))
    addInteraction!(uc, b1, b1, ones(3,3), (1,0))
    addInteraction!(uc, b1, b1, ones(3,3), (0,1))
    @test length(uc.interactions) == 2

    @test_throws ErrorException setInteractionOnsite!(uc, b1, ones(4,3))
    setInteractionOnsite!(uc, b1, ones(3,3))
    setInteractionOnsite!(uc, b2, 2.0*ones(3,3))
    @test uc.interactionsOnsite[b1] == ones(3,3)
    @test uc.interactionsOnsite[b2] == 2.0*ones(3,3)

    @test_throws ErrorException setField!(uc, b1, ones(4))
    setField!(uc, b1, ones(3))
    setField!(uc, b2, 2.0*ones(3))
    @test uc.interactionsField[b1] == ones(3)
    @test uc.interactionsField[b2] == 2.0*ones(3)
end


text:
##  a bare bones Test-Suite for some critical functions
using Test, Distributions, Random, LinearAlgebra
using DelimitedFiles

include("../src/JS_SAA_main.jl")
include("consts.jl")

@testset "allTests" begin

#Tests shrink/crossval other simpler helpers
@testset "Vector Functions Bernoulli Newsvendor" begin
	Random.seed!(8675309)
	K = 100
	ps = rand(K)
	alpha = .1
	s = .9
	p0 = .5
	Nhat = 20

	@test isapprox(JS.exp_nv_obj(ps, p0, alpha, Nhat, s),
					0.4670605847641166)
	@test isapprox(JS.exp_nv_fullInfo(ps, s),
					0.4472649361852286)
end

@testset "bestAlpha Bernoulli Newsvendor" begin
	Random.seed!(8675309)
	K = 100
	ps = rand(K)
	Nhats = rand(Poisson(20), K)
	alpha_grid = range(0, stop=80, length=80)
	ss = rand(K)/4 .+ .75
	p0 = .5
	mhats = JS.nv_sim_path(ps, Nhats)
	alphaOR, jstar, out = JS.nv_oracle_alpha(mhats, ps, p0, alpha_grid, Nhats, ss)
	@test isapprox(out[jstar], 0.4323864805448839)
	@test isapprox(alphaOR, 2.0253164556962027)

	alphaAP, jstar, out = JS.nv_apriori_alpha(ps, p0, alpha_grid, Nhats, ss)
	@test isapprox(out[jstar], 0.4338668678652777)
	@test isapprox(alphaAP, 2.025316455696202)

	alphaLOO, jstar, out = JS.nv_loo_alpha(mhats, p0, alpha_grid, Nhats, ss)
	@test isapprox(out[jstar], 7.926847617453273)
	@test isapprox(alphaLOO, 1.0126582278481013)

	p0_grid = range(0, stop=1, length=20)
	p0OR, alphaOR, jstar = JS.nv_oracle_both(mhats, ps, p0_grid, alpha_grid, Nhats, ss)
	@test isapprox(p0OR, 0.21052631578947367)
	@test isapprox(alphaOR, 14.177215189873417)
	@test jstar == 15

	p0LOO, alphaLOO, jstar = JS.nv_loo_both(mhats, p0_grid, alpha_grid, Nhats, ss)
	@test isapprox(p0LOO, 0.15789473684210525)
	@test isapprox(alphaLOO, 11.139240506329115)
	@test jstar == 12

end

@testset "simpleHelpers" begin
	K = 3
	d = 10
	N = 10
	Random.seed!(8675309)

	p0 = ones(d)/d
	ps = rand(Dirichlet(ones(d)), K)
	shrunk_p = JS.shrink(ps[:, 1], p0, .5, N)

	for i = 1:d
		@test isapprox(cShrunkP1[i], shrunk_p[i])
	end

	#generate proper data and calc mhat
	Nhats = rand(Poisson(N), K)
	mhats = JS.sim_path(ps, Nhats)

	pGM = JS.get_GM_anchor(mhats)
	for i = 1:d
		@test isapprox(cGenAnchor[i], pGM[i])
	end

	Nhats[1] = 0
	mhats = JS.sim_path(ps, Nhats)
	pGM = JS.get_GM_anchor(mhats)
	for i = 1:d
		@test isapprox(cGenAnchorZero[i], pGM[i])
	end

end

@testset "shrunkenSAA" begin
	#discrete newsvendors supported 1:d
	K = 10
	s = .95
	d = 10
	N = 20
	Random.seed!(8675309)

	#gen an "interesting" distribution of ps still centered at 1/d
	p0 = ones(d)/d
	anchor = vcat(1., zeros(d-1))
	p0 = .5 * p0 + .5 * anchor
	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	qs = rand(Dirichlet(5 * ones(d)), K-floor(Int, K/2))
	ps = [ps qs]

	alpha_grid = range(0, stop=50, length=75)
	lams = ones(K)

	Nhats = rand(Poisson(N), K)
	mhats = JS.sim_path(ps, Nhats);

	supps =  repeat(collect(1:d), outer=(1, K))
	cs = JS.getNewsVendorCosts(supps, s, K)
	xs = JS.genSSAAtrainers(supps, s)

	#First some pointwise tests
	@test isapprox(JS.z_k(xs[1], cs[:, 1], mhats[:, 1], ps[:, 1], 1, 1, (p0, 1)), 3.2530462468808903)
	@test isapprox(JS.zbar(xs, cs, mhats, ps, lams, (p0, 1)), 4.428156384118435)
	@test isapprox(JS.zstar(xs, cs, ps, lams), 4.361044771184344)

	@test isapprox(JS.zLOO_k_unsc(xs[1], cs[:, 1], mhats[:, 1], (p0, 1)), 49.0 )

	#Now check whole curve
	alphaOR, jstar, oracle_curve = JS.oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
	@test isapprox(2.027027027027027, alphaOR)
	@test isapprox(jstar, 4)

	for i = 1:length(cOracleCurve)
		@test isapprox(cOracleCurve[i], oracle_curve[i])
	end

	alphaLOO, jstar, loo_curve = JS.loo_alpha(xs, cs, mhats, p0, alpha_grid)
	@test isapprox(alphaLOO, 0.)
	@test isapprox(jstar, 1)
	for i = 1:length(cLooCurve)
		@test isapprox(cLooCurve[i], loo_curve[i])
	end

	#Check the entire cross-val curve for K-fold
	alphaCV, jstar, cv_curve = JS.cv_alpha(xs, cs, mhats, p0, alpha_grid, 5)
	@test isapprox(alphaCV, 8.783783783783784)
	@test isapprox(jstar, 14)

	for i = 1:length(cCVCurve)
		@test isapprox(cCVCurve[i], cv_curve[i])
	end

end #shrunkenSAA

@testset "KS Robust tests" begin
	#discrete newsvendors supported 1:d
	K = 10
	s = .95
	d = 10
	N = 20
	Random.seed!(8675309)

	#gen an "interesting" distribution of ps still centered at 1/d
	p0 = ones(d)/d
	anchor = vcat(1., zeros(d-1))
	p0 = .5 * p0 + .5 * anchor
	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	qs = rand(Dirichlet(5 * ones(d)), K-floor(Int, K/2))
	ps = [ps qs]

	alpha_grid = range(0, stop=50, length=75)
	lams = ones(K)
	Nhats = rand(Poisson(N), K)
	mhats = JS.sim_path(ps, Nhats);

	supps =  repeat(collect(1:d), outer=(1, K))

	#static sets
	cs = JS.getNewsVendorCosts(supps, s, K)
	xsKS = JS.genKSTrainers(supps, s, K, :aPriori)
	xsKS2 = JS.genKSTrainers(supps, s, K, :crossVal)

	#First some pointwise tests for different gamma
	@test isapprox(JS.zbar(xsKS, cs, mhats, ps, lams, 0),    4.427097278114614)
	@test isapprox(JS.zbar(xsKS, cs, mhats, ps, lams, .025), 4.38864696562911)
	@test isapprox(JS.zbar(xsKS, cs, mhats, ps, lams, .05),   4.6573883655136825)

	#one test for the cross-val version
	Gamma_grid = range(0, stop=1-s, length=51)
	@test isapprox(JS.zbar(xsKS2, cs, mhats, ps, lams, (Gamma_grid, 5)), 4.559215177976772)

end #KS Tests

@testset "optimizingAnchors" begin
	#discrete newsvendors supported 1:d
	K = 10
	s = .95
	d = 10
	N = 20
	Random.seed!(8675309)

	#gen an "interesting" distribution of ps still centered at 1/d
	p0 = ones(d)/d
	anchor = vcat(1., zeros(d-1))
	p0 = .5 * p0 + .5 * anchor
	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	qs = rand(Dirichlet(5 * ones(d)), K-floor(Int, K/2))
	ps = [ps qs]

	alpha_grid = range(0, stop=50, length=75)
	lams = ones(K)

	Nhats = rand(Poisson(N), K)
	mhats = JS.sim_path(ps, Nhats);

	supps =  repeat(collect(1:d), outer=(1, K))
	cs = JS.getNewsVendorCosts(supps, s, K)
	xs = JS.genSSAAtrainers(supps, s)

	Random.seed!(123456789)
	optp0, optAlpha, zstar = JS.loo_anchor(xs, cs, mhats, numClusters=5)
	@test isapprox(optAlpha, 17.100025815457034)
	@test isapprox(zstar, 87.0)
	for i = 1:length(cOptAnchor1)
		@test isapprox(cOptAnchor1[i], optp0[i])
	end

	optp0, optAlpha, zstar = JS.loo_anchor(xs, cs, mhats, numClusters=-1)
	@test isapprox(optAlpha, 4.174826098336609)
	@test isapprox(zstar, 86.6)
	for i = 1:length(cOptAnchor2)
		@test isapprox(cOptAnchor2[i], optp0[i])
	end

	optp0, optAlpha, zstar = JS.opt_oracle_anchor(xs, cs, ps, mhats; numClusters = 5)
	@test isapprox(optAlpha, 4.466926141845736)
	@test isapprox(zstar, 4.382226891260787)
	for i = 1:length(cOptAnchor3)
		@test isapprox(cOptAnchor3[i], optp0[i])
	end

	optp0, optAlpha, zstar = JS.opt_oracle_anchor(xs, cs, ps, mhats; numClusters = -1)
	@test isapprox(optAlpha, 3.542168883777246)
	@test isapprox(zstar, 4.361044771184344)
	for i = 1:length(cOptAnchor4)
		@test isapprox(cOptAnchor4[i], optp0[i])
	end

end #end optimizingAnchors

@testset "betaAnchors" begin
	#discrete newsvendors supported 1:d
	K = 10
	s = .95
	d = 10
	N = 20
	Random.seed!(8675309)

	#gen an "interesting" distribution of ps still centered at 1/d
	p0 = ones(d)/d
	anchor = vcat(1., zeros(d-1))
	p0 = .5 * p0 + .5 * anchor
	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	qs = rand(Dirichlet(5 * ones(d)), K-floor(Int, K/2))
	ps = [ps qs]

	alpha_grid = range(0, stop=50, length=75)
	lams = ones(K)

	Nhats = rand(Poisson(N), K)
	mhats = JS.sim_path(ps, Nhats);

	supps =  repeat(collect(1:d), outer=(1, K))
	cs = JS.getNewsVendorCosts(supps, s, K)
	xs = JS.genSSAAtrainers(supps, s)

	theta2_grid = range(1e-6, stop=3, length=5)
	mu_grid = range(1e-6, stop=1, length=5)
	alphaLOO, p0LOO, zLOO = JS.loo_betaAnchor(xs, cs, mhats, alpha_grid, theta2_grid, mu_grid)
	@test isapprox(alphaLOO, 10.81081081081081)
	@test isapprox(zLOO, 87.0)

	for i = 1:length(p0LOO)
		@test isapprox(p0LOO[i], cBetaAnchor[i])
	end

	alphaOR, p0OR, zOR = JS.oracle_betaAnchor(xs, cs, mhats, ps, lams, alpha_grid, theta2_grid, mu_grid)
	@test isapprox(alphaOR, 1.3513513513513513)
	@test isapprox(zOR, 4.361044771184344)
	for i = 1:length(p0OR)
		@test isapprox(p0OR[i], cBetaAnchorOR[i])
	end
end #end betaAnchors

end #AllTEsts


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Endpoints <: SwaggerModel
    blob::Any # spec type: Union{ Nothing, String } # spec name: blob
    queue::Any # spec type: Union{ Nothing, String } # spec name: queue
    table::Any # spec type: Union{ Nothing, String } # spec name: table
    file::Any # spec type: Union{ Nothing, String } # spec name: file
    web::Any # spec type: Union{ Nothing, String } # spec name: web
    dfs::Any # spec type: Union{ Nothing, String } # spec name: dfs
    microsoftEndpoints::Any # spec type: Union{ Nothing, StorageAccountMicrosoftEndpoints } # spec name: microsoftEndpoints
    internetEndpoints::Any # spec type: Union{ Nothing, StorageAccountInternetEndpoints } # spec name: internetEndpoints

    function Endpoints(;blob=nothing, queue=nothing, table=nothing, file=nothing, web=nothing, dfs=nothing, microsoftEndpoints=nothing, internetEndpoints=nothing)
        o = new()
        validate_property(Endpoints, Symbol("blob"), blob)
        setfield!(o, Symbol("blob"), blob)
        validate_property(Endpoints, Symbol("queue"), queue)
        setfield!(o, Symbol("queue"), queue)
        validate_property(Endpoints, Symbol("table"), table)
        setfield!(o, Symbol("table"), table)
        validate_property(Endpoints, Symbol("file"), file)
        setfield!(o, Symbol("file"), file)
        validate_property(Endpoints, Symbol("web"), web)
        setfield!(o, Symbol("web"), web)
        validate_property(Endpoints, Symbol("dfs"), dfs)
        setfield!(o, Symbol("dfs"), dfs)
        validate_property(Endpoints, Symbol("microsoftEndpoints"), microsoftEndpoints)
        setfield!(o, Symbol("microsoftEndpoints"), microsoftEndpoints)
        validate_property(Endpoints, Symbol("internetEndpoints"), internetEndpoints)
        setfield!(o, Symbol("internetEndpoints"), internetEndpoints)
        o
    end
end # type Endpoints

const _property_map_Endpoints = Dict{Symbol,Symbol}(Symbol("blob")=>Symbol("blob"), Symbol("queue")=>Symbol("queue"), Symbol("table")=>Symbol("table"), Symbol("file")=>Symbol("file"), Symbol("web")=>Symbol("web"), Symbol("dfs")=>Symbol("dfs"), Symbol("microsoftEndpoints")=>Symbol("microsoftEndpoints"), Symbol("internetEndpoints")=>Symbol("internetEndpoints"))
const _property_types_Endpoints = Dict{Symbol,String}(Symbol("blob")=>"String", Symbol("queue")=>"String", Symbol("table")=>"String", Symbol("file")=>"String", Symbol("web")=>"String", Symbol("dfs")=>"String", Symbol("microsoftEndpoints")=>"StorageAccountMicrosoftEndpoints", Symbol("internetEndpoints")=>"StorageAccountInternetEndpoints")
Base.propertynames(::Type{ Endpoints }) = collect(keys(_property_map_Endpoints))
Swagger.property_type(::Type{ Endpoints }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Endpoints[name]))}
Swagger.field_name(::Type{ Endpoints }, property_name::Symbol) =  _property_map_Endpoints[property_name]

function check_required(o::Endpoints)
    true
end

function validate_property(::Type{ Endpoints }, name::Symbol, val)
end


text:
using DifferentialEquations
using IterableTables, DataFrames
using RCall

H_comps = 4
V_comps = 3

# in sub functions, du / u are the particular relevant slices only
function F1H(du, u, p, t, βslice, I_V, N_H)
    S_H, E_H, I_H, R_H = u

    # host dynamics
    host_infection = sum(βslice .* I_V)*S_H/N_H
    host_mortality = p.μ_H .* u # include S_H, so easier to remove mortality
    host_births = sum(host_mortality)
    host_progression = p.σ_H*E_H
    recovery = p.λ*I_H

    du[1] = -host_infection + host_births
    du[2] = host_infection - host_progression
    du[3] = host_progression - recovery
    du[4] = recovery
    du[1:end] -= host_mortality
end

# in sub functions, du / u are the particular relevant slices only
function F1V(du, u, p, t, βslice, I_H, N_H)
    S_V, E_V, I_V = u
    vec_infection = sum(βslice .* I_H)*S_V/N_H
    vec_mortality = p.μ_V .* u # include S_V, so easier to remove mortality
    vec_births = sum(vec_mortality)
    vec_progression = p.σ_V*E_V

    du[1] = -vec_infection + vec_births
    du[2] = vec_infection - vec_progression
    du[3] = vec_progression
    du[1:end] -= vec_mortality
end

function F(du,u,p,t)
    dH = @view(du[1:(p.nHosts*H_comps)])
    dV = @view(du[(p.nHosts*H_comps+1):end])
    Hs = @view(u[1:(p.nHosts*H_comps)])
    Vs = @view(u[(p.nHosts*H_comps+1):end])

    I_Vs = @view(Vs[3:V_comps:V_comps*p.nVecs])
    I_Hs = @view(Hs[3:H_comps:H_comps*p.nHosts])

    for host in 0:(p.nHosts-1)
        slice = (1:H_comps).+(H_comps*host)
        F1H(@view(dH[slice]), @view(Hs[slice]), p.host[host+1], t, @view(p.β[host+1,:]), I_Vs, p.N_H)
    end
    for vec in 0:(p.nVecs-1)
        slice = (1:V_comps).+(V_comps*vec)
        F1V(@view(dV[slice]), @view(Vs[slice]), p.vec[vec+1], t, @view(p.β[:,vec+1]), I_Hs, p.N_H)
    end
end

nH = 2
nV = 2
srand(0)

S_Hs = ones(nH) .* 100.0
E_Hs = zeros(nH)
I_Hs = shuffle(vcat(zeros(nH-1),[1.0]))
R_Hs = zeros(nH)
host0 = reshape(hcat(S_Hs,E_Hs,I_Hs,R_Hs)', nH*H_comps, 1)

S_Vs = ones(nV) .* 1000.0
E_Vs = zeros(nV)
I_Vs = zeros(nV)
vec0 = reshape(hcat(S_Vs,E_Vs,I_Vs)', nV*V_comps, 1)

u0 = vcat(host0, vec0)

srand(1)

μs = 1 ./ (rand(nH) .* 360)
σs = 1 ./ (rand(nH) .* 6)
μVs = 1 ./ (rand(nV) .* 60)
σVs = 1 ./ (rand(nV) .* 14)

λs = 1 ./ (rand(nH) .* 28)
βs = rand(nH*nV) ./ 10.0

using NamedTuples
# nb: in >= Julia v0.7, can eliminate this import
#  and the @NT syntax
p = @NT(
  nHosts = nH, nVecs = nV,
  N_H = sum(host0),
  β = reshape(βs,nH,nV), # information in hosts (rows) by vectors (cols)
  vec  = [@NT(μ_V=μVs[j], σ_V=σVs[j]) for j in 1:nV],
  host = [@NT(μ_H=μs[i], σ_H=σs[i], λ=λs[i]) for i in 1:nH]
  # just building up a random collection of params for demonstration
)

tspan = (0.0, 365.0)
prob = ODEProblem(F, u0, tspan, p)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=linspace(0,365,365*10+1))

# rename!(df, Dict(:timestamp => :t,
#  :value1 => :S_H, :value2 => :E_H, :value3 => :I_H, :value4 => :R_H,
#  :value5 => :S_V, :value6 => :E_V, :value7 => :I_V
# ))
# mlt[:host] = contains.(string.(mlt[:variable]),"H"); # tag which entries are host vs vector
# df
df = DataFrame(sol)
mlt = melt(df,:timestamp) # convert results into long format for plotting
mlt[:index] = parse.(Int,replace.(string.(mlt[:variable]),r"[^\d]+"=>""))
namekey = hcat(
  reshape(["$(compartment)_H$species" for compartment in ["S","E","I","R"], species in 1:nH],1,:),
  reshape(["$(compartment)_V$species" for compartment in ["S","E","I"], species in 1:nV],1,:)
)

mlt[:name] = namekey[mlt[:index]]
mlt[:facet] = replace.(string.(mlt[:name]),r"\w+_"=>"")
mlt[:compartment] = replace.(string.(mlt[:name]),r"_\w+"=>"")
mlt

# current version RCall supports better transfers, which would simplify this mess
# but requires Julia v >= 0.7
vals = mlt[:value]
tstamps = mlt[:timestamp]
fcts = mlt[:facet]
comps = mlt[:compartment]
@rput vals tstamps fcts comps
R"
library(ggplot2)
suppressPackageStartupMessages(library(data.table))
dt <- data.table(t=tstamps, y=vals, species=fcts, compartment=comps)
ggplot(dt) + aes(x=t, y=y, color=compartment) + facet_grid(species ~ ., scale = 'free_y') +
  theme_minimal() +
  geom_line()
"

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook


text:
"""
    activate_inner!(com, constraint::ActivatorConstraint)

Activate the inner constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_inner!(com, constraint::ActivatorConstraint)
    inner_constraint = constraint.inner_constraint
    if !constraint.inner_activated
        !activate_constraint!(com, inner_constraint, inner_constraint.fct, inner_constraint.set) && return false
        constraint.inner_activated = true
        constraint.inner_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

"""
    activate_complement_inner!(com, constraint::ActivatorConstraint)

Activate the complement constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_complement_inner!(com, constraint::ActivatorConstraint)
    complement_constraint = constraint.complement_constraint
    if !constraint.complement_inner_constraint
        !activate_constraint!(com, complement_constraint, complement_constraint.fct, complement_constraint.set) && return false
        constraint.complement_inner_constraint = true
        constraint.complement_inner_constraint_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

"""
    update_best_bound_constraint!(com::CS.CoM,
        constraint::ActivatorConstraint,
        fct::Union{MOI.VectorOfVariables, VAF{T}},
        set::IS,
        vidx::Int,
        lb::Int,
        ub::Int
    ) where {T<:Real}

Update the bound constraint associated with this constraint. This means that the `bound_rhs` bounds will be changed according to
the possible values the table constraint allows. `vidx`, `lb` and `ub` don't are not considered atm.
Additionally only a rough estimated bound is used which can be computed relatively fast.
This method calls the inner_constraint method if it exists and the indicator is activated.
"""
function update_best_bound_constraint!(
    com::CS.CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    vidx::Int,
    lb::Int,
    ub::Int,
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    activator_vidx = constraint.indices[1]
    search_space = com.search_space
    activator_var = search_space[activator_vidx]
    if CS.issetto(activator_var, Int(constraint.activate_on))
        return update_best_bound_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            vidx,
            lb,
            ub,
        )
    else
        # if not activated (for example in a different subtree) we reset the bounds
        for rhs in constraint.bound_rhs
            rhs.lb = typemin(Int64)
            rhs.ub = typemax(Int64)
        end
    end
end

function single_reverse_pruning_constraint!(
    com::CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    var::Variable,
    backtrack_idx::Int,
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    # the variable must be part of the inner constraint
    if (var.idx != constraint.indices[1] || constraint.activator_in_inner)
        single_reverse_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            var,
            backtrack_idx,
        )
    end
end

function reverse_pruning_constraint!(
    com::CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    backtrack_id::Int,
) where {
    T<:Real,
}
    # check if inner constraint should be deactived again
    if constraint.inner_activated && backtrack_id == constraint.inner_activated_in_backtrack_idx
        constraint.inner_activated = false
        constraint.inner_activated_in_backtrack_idx = 0
    end
    if constraint isa ReifiedConstraint
        if constraint.complement_inner_constraint && backtrack_id == constraint.complement_inner_constraint_in_backtrack_idx
            constraint.complement_inner_constraint = false
            constraint.complement_inner_constraint_in_backtrack_idx = 0
        end
    end

    inner_constraint = constraint.inner_constraint
    reverse_pruning_constraint!(
        com,
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set,
        backtrack_id,
    )
end

function restore_pruning_constraint!(
    com::CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    prune_steps::Union{Int,Vector{Int}},
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    restore_pruning_constraint!(
        com,
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set,
        prune_steps,
    )
end

function finished_pruning_constraint!(
    com::CS.CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    finished_pruning_constraint!(
        com,
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set,
    )
end


text:
import Base: intersect

"""
    intersect(line_a::AbstractLine, line_b::AbstractLine) -> AbstractVector

Compute the intersection of two lines.

The intersection of two lines is a point.

# Examples

```jldoctest
julia> intersect(Line([0, 0], [1, 0]), Line([5, 5], [0, 1]))
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 5.0
 0.0
```

"""
function intersect(line_a::AbstractLine, line_b::AbstractLine)

    if !are_coplanar(line_a, line_b)
        throw(ArgumentError("The lines are not coplanar."))
    end

    if are_parallel(line_a.direction, line_b.direction)
        throw(ArgumentError("The lines are parallel."))
    end

    # Vector from line A to line B.
    vector_ab = Vector(line_a.point, line_b.point)

    # Vector perpendicular to both lines.
    vector_perpendicular = line_a.direction × line_b.direction

    num = (vector_ab × line_b.direction) ⋅ vector_perpendicular
    denom = norm(vector_perpendicular)^2

    # Vector along line A to the intersection point.
    vector_a_scaled = num / denom * line_a.direction

    return line_a.point + vector_a_scaled
end

"""
    intersect(plane::AbstractPlane, line::AbstractLine; kwargs...) -> AbstractVector

Compute the intersection point of a plane and a line.

An error is thrown if the plane and line are parallel
(i.e., the plane normal is perpendicular to the line direction).
The tolerances for this check can be passed as keyword arguments.

# Examples

```jldoctest
julia> intersect(Plane([0,0,0], [0, 0, 1]), Line([0, 0, 1], [1, 1, 1]))
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 -1.0
 -1.0
  0.0

julia> intersect(Plane([0,0,0], [0, 0, 1]), Line([0, 0, 1], [1, 1, 0]))
ERROR: ArgumentError: The line and plane are parallel.
```

Tolerances for the parallel check can be passed as keyword arguments.

```
julia> plane = Plane([0, 0, 0], [0, 0, 1]);

julia> line = Line([0, 0, 1], [1, 1, 1e-3]);

julia> intersect(plane, line)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 -1000.0
 -1000.0
     0.0

julia> intersect(plane, line, atol=1e-3)
ERROR: ArgumentError: The line and plane are parallel.

julia> intersect(plane, line, rtol=1e-3)
ERROR: ArgumentError: The line and plane are parallel.
```
"""
function intersect(plane::AbstractPlane, line::AbstractLine; kwargs...)

    if are_perpendicular(line.direction, plane.normal; kwargs...)
        throw(ArgumentError("The line and plane are parallel."))
    end

    vector_plane_line = Vector(plane.point, line.point)

    num = -plane.normal ⋅ vector_plane_line
    denom = plane.normal ⋅ line.direction

    # Vector along the line to the intersection point.
    vector_line_scaled = num / denom * line.direction

    return line.point + vector_line_scaled
end

"""
    intersect(plane_a::AbstractPlane, plane_b::AbstractPlane; kwargs...) -> AbstractLine

Compute the intersection of two planes. The intersection of the planes is a line.

An error is thrown if the planes are parallel.
The tolerances for this check can be passed as keyword arguments.

# Examples

```jldoctest
julia> plane_a = Plane([0, 0, 0], [0, 0, 1]);

julia> plane_b = Plane([5, 3, 2], [1, 0, 0]);

julia> intersect(plane_a, plane_b)
Line{3, Float64}([5.0, 0.0, 0.0], [0.0, 1.0, 0.0])

julia> plane_a = Plane([0, 0, 0], [0, 0, 1]);

julia> plane_b = Plane([1, 1, 1], [0, 0, -3]);

julia> intersect(plane_a, plane_b)
ERROR: ArgumentError: The planes are parallel.

julia> plane_a = Plane([0, 0, 0], [0, 0, 1]);

julia> plane_b = Plane([1, 1, 1], [0, 1e-2, 1]);

julia> intersect(plane_a, plane_b, atol=1e-3)
ERROR: ArgumentError: The planes are parallel.
```

"""
function intersect(plane_a::AbstractPlane, plane_b::AbstractPlane; kwargs...)

    if are_parallel(plane_a.normal, plane_b.normal; kwargs...)
        throw(ArgumentError("The planes are parallel."))
    end

    array_normals_stacked = vcat(plane_a.normal', plane_b.normal')

    array_11 = 2 * Matrix(1I, 3, 3)
	array_12 = array_normals_stacked'
	array_21 = array_normals_stacked
	array_22 = zeros(2, 2)

	matrix = vcat(hcat(array_11, array_12), hcat(array_21, array_22))

    dot_a = plane_a.point ⋅ plane_a.normal
    dot_b = plane_b.point ⋅ plane_b.normal

    y = [0, 0, 0, dot_a, dot_b]

    solution = matrix \ y

    point = solution[1:3]
	direction = plane_a.normal × plane_b.normal

    return Line(point, direction)
end

"""
    intersect(circle::Circle, line::AbstractLine) -> Tuple{StaticArrays.SVector, StaticArrays.SVector}

Compute the intersection of a circle and a 2D line.

The intersection is returned as two points (which can be the same).
An error is thrown if the circle and line do not intersect.

# Examples

```jldoctest
julia> circle = Circle([0, 0], 1);

julia> intersect(circle, Line([0, 0], [1, 0]))
([-1.0, 0.0], [1.0, 0.0])

julia> intersect(circle, Line([0, 1], [1, 0]))
([0.0, 1.0], [0.0, 1.0])

julia> intersect(circle, Line([0, 2], [1, 0]))
ERROR: ArgumentError: The circle and line do not intersect.
```

"""
function intersect(circle::Circle, line::AbstractLine)

    # Two points on the line.
    point_1 = line.point
    point_2 = point_1 + unit(line.direction)

    # Translate the points on the line to mimic the circle being centered on the origin.
    point_translated_1 = point_1 - circle.point
    point_translated_2 = point_2 - circle.point

    x_1, y_1 = point_translated_1
    x_2, y_2 = point_translated_2

    d_x = x_2 - x_1
    d_y = y_2 - y_1

    # Pre-compute variables common to x and y equations.
    d_r_squared = d_x^2 + d_y^2
    determinant = x_1 * y_2 - x_2 * y_1
    discriminant = circle.radius^2 * d_r_squared - determinant^2

    if discriminant < 0
        throw(ArgumentError("The circle and line do not intersect."))
    end

    root = √discriminant

    mp = [-1, 1]  # Array to compute minus/plus.
    sign = d_y < 0 ? -1 : 1

    coords_x = (determinant * d_y .+ mp * sign * d_x * root) / d_r_squared
    coords_y = (-determinant * d_x .+ mp * abs(d_y) * root) / d_r_squared

    point_translated_a = [coords_x[1], coords_y[1]]
    point_translated_b = [coords_x[2], coords_y[2]]

    # Translate the intersection points back from origin circle to real circle.
    point_a = point_translated_a + circle.point
    point_b = point_translated_b + circle.point

    return point_a, point_b
end

"""
    intersect(sphere::Sphere, line::AbstractLine) -> Tuple{StaticArrays.SVector, StaticArrays.SVector}

Compute the intersection of a sphere and a 3D line.

The intersection is returned as two points (which can be the same).
An error is thrown if the sphere and line do not intersect.

# Examples

```jldoctest
julia> sphere = Sphere([0, 0, 0], 1);

julia> intersect(sphere, Line([0, 0, 0], [1, 0, 0]))
([-1.0, 0.0, 0.0], [1.0, 0.0, 0.0])

julia> intersect(sphere, Line([0, 0, -1], [1, 0, 0]))
([0.0, 0.0, -1.0], [0.0, 0.0, -1.0])

julia> intersect(sphere, Line([0, 1, 1], [1, 0, 0]))
ERROR: ArgumentError: The sphere and line do not intersect.
```

"""
function intersect(sphere::Sphere, line::AbstractLine)

    vector_to_line = Vector(sphere.point, line.point)
    direction_unit = unit(line.direction)

    dot = direction_unit ⋅ vector_to_line

    discriminant = dot^2 - (norm(vector_to_line)^2 - sphere.radius^2)

    if discriminant < 0
        throw(ArgumentError("The sphere and line do not intersect."))
    end

    pm = [-1, 1]  # Array to compute minus/plus.
    distances = -dot .+ √discriminant * pm

    points = line.point .+ distances' .* direction_unit

    point_a = points[:, 1]
    point_b = points[:, 2]

    return point_a, point_b
end