text:
using Photometry.Aperture: circular_overlap,
                           circular_overlap_core,
                           circular_overlap_single_exact,
                           circular_overlap_single_subpixel,
                           area_arc,
                           area_triangle

@testset "circular overlap" for grid_size in [50, 500, 1000], circ_size in (0.2, 0.4, 0.8), method in [:exact, :center, (:subpixel, 2), (:subpixel, 5), (:subpixel, 10)]

    g = circular_overlap(-1, -1, 1, 1, grid_size, grid_size, circ_size, method = method)
    any(g .> 0) &&  @test maximum(g) ≈ 1.0
    any(g .< 1) &&  @test minimum(g) ≈ 0.0
end

@testset "area triangle" begin
    @test area_triangle(0, 0, 0, 0, 0, 0) ≈ 0

    @test area_triangle(0, 0, 1, 0, 0, 1) ≈ 1 / 2

    @test area_triangle(1, 1, 1, 4, 6, 1) ≈ 7.5

end

@testset "area arc (r = $r)" for r in 1:10
    @test area_arc(0, 0, 0, 0, r) ≈ 0
    @test area_arc(0, r, r, 0, r) ≈ r^2 / 2 * (π / 2 - 1)
end

@testset "overlap core (r = $r)" for r in 1:10
    @test circular_overlap_core(0, 0, 0, 0, r) ≈ 0
    @test circular_overlap_single_exact(0, 0, r, r, r) ≈ r^2 * π / 4
    @test circular_overlap_core(0, 0, r, r, r) ≈ r^2 * π / 4
end

@testset "overlap single exact" begin
    @test circular_overlap_single_exact(0, 0, 0, 0, 1) ≈ 0
    @test circular_overlap_single_exact(0, 0, 20, 20, 15) ≈ 176.714586764429
end

@testset "overlap subpixel" begin
    @test circular_overlap_single_subpixel(0, 0, 0, 0, 1, 1) ≈ 1

    @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 2) ≈ 0.25
    @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 5) ≈ 0.44
    @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 10) ≈ 0.43
    @test circular_overlap_single_subpixel(0, 0, 20, 20, 15, 100) ≈ 0.4423
end

@testset "overlap grid" begin
    @test circular_overlap(0, 0, 0, 0, 1, 1, 1) == ones(1, 1)

    @test circular_overlap(0, 20, 0, 20, 5, 5, 15) ≈ [1 1 1 0.70506901 0;
    1 1 1 0.42240835 0;
    1 1 0.74542733 0.02213982 0;
    0.70506901 0.42240835 0.02213982 0 0;
    0 0 0 0 0]

    @test circular_overlap(-10, 10, -10, 10, 5, 5, 4) ≈ [0 0         0        0        0;
    0 0.07878669 0.45661148 0.07878669 0;
    0 0.45661148 1         0.45661148 0;
    0 0.07878669 0.45661148 0.07878669 0;
    0 0         0         0        0]

end


text:
module Generator
include("constant.jl")

using Printf
using .Constant: NONO_SIZE

"""
Generate "number" children with "encodingFunc", return arrays of rows permutation according to number\n
Return a two-dimension array, which represents row as child number, col as permutated result of nonogram's rows
"""
function generate_children!(rows::Array{Array{Int,1},1}, number::Int, encodingFunc::Function, args...)::Array{UInt128,2}
    if number <= 0
        throw(ArgumentError("Number must be positive number!"))
    end
    permRows = Array{UInt128,1}(undef, 0)
    # record the encoding ignored indice of first child
    ignoreIndice = Array{Int}(undef, 0)
    for (i, row) in enumerate(rows)
        ans, isDetermine = encodingFunc(row, args...)
        if isDetermine
            push!(ignoreIndice, i)
        end
        push!(permRows, ans)
    end
    # after first child's children will follow the ignored indice
    for n in 1:(number - 1), (i, row) in enumerate(rows)
        if i in ignoreIndice
            # debug feature
            # println(@sprintf("Ignored %s index at %d in %d child.", row, i, n + 1))
            push!(permRows, permRows[i])
            continue
        end
        push!(permRows, encodingFunc(row, args...)[1])
    end
    # reshape to fit size(children numbers, nonogram size), need transpose because column-major
    return transpose(reshape(permRows, NONO_SIZE, :))
end
end

ext: jl


text:
using OrdinaryDiffEq, ParameterizedFunctions, DiffEqSensitivity

A = reshape([1,0,2,3],2,2)
function f_morris(p)
    A*p
end
m = DiffEqSensitivity.morris_sensitivity(f_morris,[[1,5],[1,5]],[10,10],len_trajectory=1500,total_num_trajectory=1000,num_trajectory=150)
@test m.means[1] ≈ A[:,1] atol=1e-12
@test m.means[2] ≈ A[:,2] atol=1e-12
@test m.variances ≈ [[0,0],[0,0]] atol=1e-12

m = DiffEqSensitivity.morris_sensitivity(f_morris,[[1,5],[1,5]],[10,10],relative_scale=true,len_trajectory=2000,total_num_trajectory=1000,num_trajectory=150)
@test m.means[1][2] ≈ 0 atol=1e-12
@test m.variances[1][2] ≈ 0 atol=1e-12
@test m.means[2][1] < m.means[2][2]


text:
function delete_cvrp(filepath)
    rm(filepath)
end

function write_cvrp(cvrp::CVRPLIB.CVRP)

    if cvrp.name == ""
        name = randstring(10)
    else
        name = cvrp.name
    end

    filepath = joinpath(pwd(), name * ".vrp")

    if cvrp.dummy > cvrp.depot # dummy is added
        org_weights = cvrp.weights[1:end-1, 1:end-1]
    else
        org_weights = cvrp.weights
    end

    open(filepath, "w") do io
        write(io, "NAME : $(name)\n")
        write(io, "COMMENT : Generated by BranchPriceCutCVRP.jl\n")
        write(io, "TYPE : CVRP\n")
        write(io, "DIMENSION : $(cvrp.dimension)\n")
        write(io, "EDGE_WEIGHT_TYPE : $(cvrp.weight_type)\n")
        write(io, "CAPACITY : $(cvrp.capacity)\n")
        if cvrp.distance < Inf
            write(io, "DISTANCE : $(cvrp.distance)\n");
        end
        if cvrp.service_time > 0.0
            write(io, "SERVICE_TIME : $(cvrp.service_time)\n");
        end

        if cvrp.weight_type == "EXPLICIT"
            @assert length(cvrp.coordinates) == 0

            write(io, "EDGE_WEIGHT_FORMAT : FULL_MATRIX\n")
            write(io, "DISPLAY_DATA_TYPE : NO_DISPLAY\n")
            write(io, "EDGE_WEIGHT_SECTION\n")
            writedlm(io, org_weights)

            # write(io, "EDGE_WEIGHT_FORMAT : LOWER_ROW\n")
            # write(io, "DISPLAY_DATA_TYPE : NO_DISPLAY\n")
            # write(io, "EDGE_WEIGHT_SECTION\n")
            # for i in 2:cvrp.dimension
            #     for j in 1:i-1
            #         write(io, "\t$(org_weights[i,j])")
            #     end
            # end
            # write(io, "\n")

        else
            write(io, "NODE_COORD_SECTION\n")
            x = cvrp.coordinates[:, 1]
            y = cvrp.coordinates[:, 2]
            for i in 1:cvrp.dimension
                write(io, "$i $(x[i]) $(y[i])\n")
            end
        end
        write(io, "DEMAND_SECTION\n")
        for i in 1:cvrp.dimension
            write(io, "$i $(cvrp.demand[i])\n")
        end
        write(io, "DEPOT_SECTION\n")
        write(io, "$(cvrp.depot)\n")
        write(io, "-1\n")
        write(io, "EOF\n")
    end

    return name, filepath
end


text:
#User OAuth data goes here
USER_AUTH_DICT = Dict(
  :consumer_key => "XXXXXXXXX",
  :consumer_secret => "XXXXXXXXX",
  :access_token => "XXXXXXXXX",
  :access_secret => "XXXXXXXXX"
)

#DBHandles
#Comment out to de-activate

include("dbhandle/TwTXTFile.jl")
#include("dbhandle/TwMongoDB.jl")

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
module ScikitLearn_

#> export the new models you're going to define (and nothing else):
export SVMClassifier, SVMRegressor
export SVMNuClassifier, SVMNuRegressor
export SVMLClassifier, SVMLRegressor
export SCElasticNet, SCElasticNetCV

#> for all Supervised models:
import MLJBase

#> for all classifiers:
using CategoricalArrays

#> import package:
import ..ScikitLearn: @sk_import
import ..ScikitLearn
@sk_import svm: SVC
@sk_import svm: NuSVC
@sk_import svm: LinearSVC

@sk_import svm: SVR
@sk_import svm: NuSVR
@sk_import svm: LinearSVR
@sk_import linear_model : ElasticNet
@sk_import linear_model : ElasticNetCV

"""
    SVMClassifier(; kwargs...)
C-Support Vector classifier from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html#sklearn.svm.SVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html#sklearn.svm.SVC). Implemented hyperparameters as per
package documentation cited above.
See also, SVMNuClassifier, SVMLClassifier, SVMRegressor
"""
mutable struct SVMClassifier <: MLJBase.Deterministic
    C::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
    decision_function_shape::String
    random_state
end

# constructor:
#> all arguments are kwargs with a default value
function SVMClassifier(
    ;C=1.0
    ,kernel="rbf"
    ,degree=3
    ,gamma="auto"
    ,coef0=0.0
    ,shrinking=true
    ,tol=1e-3
    ,cache_size=200
    ,max_iter=-1
    ,decision_function_shape="ovr"
    ,random_state=nothing)

    model = SVMClassifier(
        C
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter
        , decision_function_shape
        , random_state
        )

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SVMClassifier)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
	    model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
	    model.gamma="auto"
    end
    if(!(model.decision_function_shape in ["ovo","ovr"]))
            warning *="decision_function_shape parameter is not valid, setting to default=\"ovr\" \n"
	    model.decision_function_shape="ovr"
    end
    return warning
end

"""
    SVMNuClassifier(; kwargs...)
NU-Support Vector classifier from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVC.html#sklearn.svm.NuSVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVC.html#sklearn.svm.NuSVC). Implemented hyperparameters as per
package documentation cited above.
See also, SVMClassifier, SVMLClassifier, SVMNuRegressor
"""
mutable struct SVMNuClassifier <: MLJBase.Deterministic
    nu::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
    decision_function_shape::String
    random_state
end

# constructor:
#> all arguments are kwargs with a default value
function SVMNuClassifier(
    ;nu=0.5
    ,kernel="rbf"
    ,degree=3
    ,gamma="auto"
    ,coef0=0.0
    ,shrinking=true
    ,tol=1e-3
    ,cache_size=200
    ,max_iter=-1
    ,decision_function_shape="ovr"
    ,random_state=nothing)

    model = SVMNuClassifier(
        nu
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter
        , decision_function_shape
        , random_state
        )

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SVMNuClassifier)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
	    model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
	    model.gamma="auto"
    end
    if(!(model.decision_function_shape in ["ovo","ovr"]))
            warning *="decision_function_shape parameter is not valid, setting to default=\"ovr\" \n"
	    model.decision_function_shape="ovr"
    end
    return warning
end

"""
    SVMLClassifier(; kwargs...)
Linear-support Vector classifier from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVC.html#sklearn.svm.LinearSVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVC.html#sklearn.svm.LinearSVC). Implemented hyperparameters as per
package documentation cited above.
See also, SVMClassifier, SVMNuClassifier, SVMLRegressor
"""

mutable struct SVMLClassifier <: MLJBase.Deterministic
    C::Float64
    loss::String
    dual::Bool
    penalty::String
    tol::Float64
    max_iter::Int
    intercept_scaling::Float64
    random_state
end

# constructor:
#> all arguments are kwargs with a default value

function SVMLClassifier(
    ;C=1.0
    ,loss="squared_hinge"
    ,dual=true
    ,penalty="l2"
    ,tol=1e-3
    ,max_iter=-1
    ,intercept_scaling=1.
    ,random_state=nothing)

    model = SVMLClassifier(
        C
        , loss
	, dual
	, penalty
        , tol
        , max_iter
	, intercept_scaling
        , random_state
        )

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SVMLClassifier)
    warning = ""
    if(!(model.loss in ["hinge","squared_hinge"]))
            warning *="loss parameter is not valid, setting to default=\"squared_hinge\" \n"
	    model.loss="squared_hinge"
    end
    if(!(model.penalty in ["l1","l2"]))
            warning *="penalty parameter is not valid, setting to default=\"l2\" \n"
	    model.penalty="l2"
    end
    return warning
end

"""
    SVMRegressor(; kwargs...)
Epsilon-Support Vector Regression from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html#sklearn.svm.SVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html#sklearn.svm.SVR). Implemented hyperparameters as per
package documentation cited above.
See also, SVMClassifier, SVMNuRegressor, SVMLRegressor
"""
mutable struct SVMRegressor <: MLJBase.Deterministic
    C::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
    epsilon::Float64
end

# constructor:
#> all arguments are kwargs with a default value
function SVMRegressor(
    ;C=1.0
    ,kernel="rbf"
    ,degree=3
    ,gamma="auto"
    ,coef0=0.0
    ,shrinking=true
    ,tol=1e-3
    ,cache_size=200
    ,max_iter=-1
    ,epsilon=0.1)

    model = SVMRegressor(
        C
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter
        , epsilon)

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SVMRegressor)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
	    model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
	    model.gamma="auto"
    end
    return warning
end

"""
    SVMNuRegressor(; kwargs...)
Nu Support Vector Regression from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVR.html#sklearn.svm.NuSVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVR.html#sklearn.svm.NuSVR). Implemented hyperparameters as per
package documentation cited above.
See also, SVMNuClassifier, SVMRegressor, SVMLRegressor
"""

mutable struct SVMNuRegressor <: MLJBase.Deterministic
    nu::Float64
    C::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
end

# constructor:
#> all arguments are kwargs with a default value
function SVMNuRegressor(
    ;nu=0.5
    ,C=1.0
    ,kernel="rbf"
    ,degree=3
    ,gamma="auto"
    ,coef0=0.0
    ,shrinking=true
    ,tol=1e-3
    ,cache_size=200
    ,max_iter=-1)

    model = SVMNuRegressor(
        nu
	, C
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter)

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SVMNuRegressor)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
	    model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
	    model.gamma="auto"
    end
    return warning
end

"""
    SVMLRegressor(; kwargs...)
Linear Support Vector Regression from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR). Implemented hyperparameters as per
package documentation cited above.
See also, SVMRegressor, SVMNuRegressor, SVMLClassifier
"""

mutable struct SVMLRegressor <: MLJBase.Deterministic
    C::Float64
    loss::String
    fit_intercept::Bool
    dual::Bool
    tol::Float64
    max_iter::Int
    epsilon::Float64
end

# constructor:
#> all arguments are kwargs with a default value
function SVMLRegressor(
    ;C=1.0
    ,loss="epsilon_insensitive"
    ,fit_intercept=true
    ,dual=true
    ,tol=1e-3
    ,max_iter=-1
    ,epsilon=0.1)

    model = SVMLRegressor(
        C
	, loss
	, fit_intercept
	, dual
        , tol
        , max_iter
        , epsilon)

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SVMLRegressor)
    warning = ""
    if(!(model.loss in ["epsilon_insensitive","squared_epsilon_insensitive"]))
            warning *="loss parameter is not valid, setting to default=\"epsilon_insensitive\" \n"
	    model.loss="epsilon_insensitive"
    end
    return warning
end

function MLJBase.fit(model::SVMClassifier
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    y_plain = MLJBase.int(y)
    decode  = MLJBase.decoder(y[1]) # for predict method

    cache = SVC(C=model.C,
            kernel=model.kernel,
            degree=model.degree,
            coef0=model.coef0,
            shrinking=model.shrinking,
            gamma=model.gamma,
            tol=model.tol,
            cache_size=model.cache_size,
            max_iter=model.max_iter,
            decision_function_shape=model.decision_function_shape,
            random_state=model.random_state
    )

    result = ScikitLearn.fit!(cache, Xmatrix, y_plain)
    fitresult = (result, decode)
    report = NamedTuple()

    return fitresult, nothing, report

end

function MLJBase.fit(model::SVMNuClassifier
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    y_plain = MLJBase.int(y)
    decode  = MLJBase.decoder(y[1]) # for predict method

    cache = NuSVC(nu=model.nu,
            kernel=model.kernel,
            degree=model.degree,
            coef0=model.coef0,
            shrinking=model.shrinking,
	    gamma=model.gamma,
            tol=model.tol,
            cache_size=model.cache_size,
            max_iter=model.max_iter,
            decision_function_shape=model.decision_function_shape,
            random_state=model.random_state
    )

    result = ScikitLearn.fit!(cache,Xmatrix,y_plain)
    fitresult = (result, decode)
    report = NamedTuple()

    return fitresult, nothing, report

end

function MLJBase.fit(model::SVMLClassifier
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    y_plain = MLJBase.int(y)
    decode  = MLJBase.decoder(y[1]) # for predict method

    cache = LinearSVC(C=model.C,
	    loss = model.loss,
            dual=model.dual,
            penalty=model.penalty,
            intercept_scaling=model.intercept_scaling,
            tol=model.tol,
            max_iter=model.max_iter,
            random_state=model.random_state
    )

    result = ScikitLearn.fit!(cache,Xmatrix,y_plain)
    fitresult = (result, decode)
    report = NamedTuple()

    return fitresult, nothing, report

end

function MLJBase.fit(model::SVMRegressor
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    cache = SVR(C=model.C,
            kernel=model.kernel,
            degree=model.degree,
            coef0=model.coef0,
            shrinking=model.shrinking,
	    gamma=model.gamma,
            tol=model.tol,
            cache_size=model.cache_size,
            max_iter=model.max_iter,
            epsilon=model.epsilon)

    fitresult = ScikitLearn.fit!(cache,Xmatrix,y)
    report = NamedTuple()

    return fitresult, nothing, report
end

function MLJBase.fit(model::SVMNuRegressor
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    cache = NuSVR(nu=model.nu,
            C=model.C,
            kernel=model.kernel,
            degree=model.degree,
            coef0=model.coef0,
            shrinking=model.shrinking,
	    gamma=model.gamma,
            tol=model.tol,
            cache_size=model.cache_size,
            max_iter=model.max_iter)

    fitresult = ScikitLearn.fit!(cache,Xmatrix,y)
    report = NamedTuple()

    return fitresult, nothing, report
end

function MLJBase.fit(model::SVMLRegressor
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    cache = LinearSVR(C=model.C,
            loss=model.loss,
	    fit_intercept=model.fit_intercept,
	    dual=model.dual,
	    tol=model.tol,
	    max_iter=model.max_iter,
	    epsilon=model.epsilon)

    fitresult = ScikitLearn.fit!(cache,Xmatrix,y)
    report = NamedTuple()

    return fitresult, nothing, report
end

# placeholder types for predict dispatching
SVMC = Union{SVMClassifier, SVMNuClassifier, SVMLClassifier}
SVMR = Union{SVMRegressor, SVMNuRegressor, SVMLRegressor}
SVM = Union{SVMC, SVMR}

function MLJBase.predict(model::SVMC
                         , fitresult
                         , Xnew)

    xnew = MLJBase.matrix(Xnew)
    result, decode = fitresult
    prediction = ScikitLearn.predict(result, xnew)
    return decode(prediction)
end

function MLJBase.predict(model::SVMR
                         , fitresult
                         , Xnew)
    xnew = MLJBase.matrix(Xnew)
    prediction = ScikitLearn.predict(fitresult,xnew)
    return prediction
end

"""
   SCElasticNet(; kwargs...)
ElasticNet from
[https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html). Implemented hyperparameters as per
package documentation cited above.

"""
mutable struct SCElasticNet <: MLJBase.Deterministic
    alpha::Float64
    l1_ratio::Float64
    fit_intercept::Bool
    normalize::Bool
    precompute::Bool
    max_iter::Int
    copy_X::Bool
    tol::Float64
    warm_start::Bool
    positive::Bool
    selection::String
end

# constructor:
#> all arguments are kwargs with a default value
function SCElasticNet(
    ;alpha=1.0
    ,l1_ratio = 0.5
    ,fit_intercept = true
    ,normalize=false
    ,precompute=false
    ,max_iter=1000
    ,copy_X=true
    ,tol=0.0001
    ,warm_start=false
    ,positive=false
    ,selection="cyclic")

    model = SCElasticNet(
        alpha
        , l1_ratio
        , fit_intercept
        , normalize
        , precompute
        , max_iter
        , copy_X
        , tol
        , warm_start
        , positive
        , selection
        )

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SCElasticNet)
    warning = ""
	if(model.alpha<0.0)
		warning *="alpha must be stricly positive, set to 1"
		model.alpha=1
	end
	if(model.l1_ratio!=nothing)
		for (iter,val) in enumerate(model.l1_ratio)
			if(!(0<val<=1))
				warning *="l1 must be in (0,1], set to 1"
				model.l1_ratio[iter]=1
			end
		end
	end
	return warning
end

function MLJBase.fit(model::SCElasticNet
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    cache = ElasticNet(alpha=model.alpha,
            l1_ratio=model.l1_ratio,
            fit_intercept=model.fit_intercept,
            normalize=model.normalize,
            precompute=model.precompute,
            max_iter=model.max_iter,
	    copy_X=model.copy_X,
            tol=model.tol,
            warm_start=model.warm_start,
            positive=model.positive )

    result = ScikitLearn.fit!(cache,Xmatrix,y)
    fitresult = result
	report = NamedTuple{(:n_iter,:dual_gap)}((fitresult.n_iter_, fitresult.dual_gap_))
    return fitresult, nothing, report

end

function MLJBase.fitted_params(model::SCElasticNet, fitresult)
	 return NamedTuple{(:intercept,:coef)}((fitresult.intercept_,fitresult.coef_))
end

function MLJBase.predict(model::SCElasticNet
                         , fitresult
                         , Xnew)
    xnew = MLJBase.matrix(Xnew)
    prediction = ScikitLearn.predict(fitresult,xnew)
    return prediction
end

"""
   SCElasticNetCV(; kwargs...)
ElasticNetCV from
[https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNetCV.html](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNetCV.html). Implemented hyperparameters as per
package documentation cited above.

"""
mutable struct SCElasticNetCV <: MLJBase.Deterministic
    l1_ratio
	eps::Float64
	n_alphas::Int
	alphas
	fit_intercept::Bool
	normalize::Bool
	precompute
	max_iter::Int
	tol::Float64
	cv
	copy_X::Bool
	positive::Bool
	selection::String
end

# constructor:
#> all arguments are kwargs with a default value
function SCElasticNetCV(
    ;l1_ratio = 0.5
	,eps=0.001
	,n_alphas=100
	,alphas=nothing
    ,fit_intercept = true
    ,normalize=false
    ,precompute=false
    ,max_iter=1000
    ,tol=0.0001
	,cv=nothing
	,copy_X=true
    ,positive=false
    ,selection="cyclic")

    model = SCElasticNetCV(
		l1_ratio
		, eps
		, n_alphas
		, alphas
		, fit_intercept
		, normalize
		, precompute
		, max_iter
		, tol
		, cv
		, copy_X
		, positive
		, selection
        )

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJBase.clean!(model::SCElasticNetCV)
    warning = ""
	if(model.alphas!=nothing)
		if(model.alphas<0.0)
			warning *="alpha must be stricly positive, set to 1"
			model.alphas=1
		end
	end
	if(model.l1_ratio!=nothing)
		for (iter,val) in enumerate(model.l1_ratio)
			if(!(0<val<=1))
				warning *="l1 must be in (0,1], set to 1"
				model.l1_ratio[iter]=1
			end
		end
	end
	return warning
end

function MLJBase.fit(model::SCElasticNetCV
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MLJBase.matrix(X)

    cache = ElasticNetCV(l1_ratio=model.l1_ratio
						, eps=model.eps
						, n_alphas = model.n_alphas
						, alphas = model.alphas
						, fit_intercept = model.fit_intercept
						, normalize = model.normalize
						, precompute = model.precompute
						, max_iter = model.max_iter
						, tol = model.tol
						, cv = model.cv
						, copy_X = model.copy_X
						, positive = model.positive
						, selection = model.selection
						)

    result = ScikitLearn.fit!(cache,Xmatrix,y)
    fitresult = result
    report = NamedTuple{(:l1_ratio,:alpha,:n_iter,:dual_gap)}((fitresult.l1_ratio_,fitresult.alpha_,fitresult.n_iter_,fitresult.dual_gap_))

    return fitresult, nothing, report

end

function MLJBase.fitted_params(model::SCElasticNetCV, fitresult)
	 return NamedTuple{(:intercept,:coef)}((fitresult.intercept_,fitresult.coef_))
end

function MLJBase.predict(model::SCElasticNetCV
                         , fitresult
                         , Xnew)
    xnew = MLJBase.matrix(Xnew)
    prediction = ScikitLearn.predict(fitresult,xnew)
    return prediction
end

# metadata:
MLJBase.load_path(::Type{<:SVMClassifier}) = "MLJModels.ScikitLearn_.SVMClassifier"
MLJBase.load_path(::Type{<:SVMNuClassifier}) = "MLJModels.ScikitLearn_.SVMNuClassifier"
MLJBase.load_path(::Type{<:SVMLClassifier}) = "MLJModels.ScikitLearn_.SVMLClassifier"
MLJBase.load_path(::Type{<:SVMRegressor}) = "MLJModels.ScikitLearn_.SVMRegressor"
MLJBase.load_path(::Type{<:SVMNuRegressor}) = "MLJModels.ScikitLearn_.SVMNuRegressor"
MLJBase.load_path(::Type{<:SVMRegressor}) = "MLJModels.ScikitLearn_.SVMRegressor"
MLJBase.load_path(::Type{<:SVMLRegressor}) = "MLJModels.ScikitLearn_.SVMLRegressor"

MLJBase.package_name(::Type{<:SVM}) = "ScikitLearn"
MLJBase.package_uuid(::Type{<:SVM}) = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
MLJBase.is_pure_julia(::Type{<:SVM}) = false
MLJBase.package_url(::Type{<:SVM}) = "https://github.com/cstjean/ScikitLearn.jl"
MLJBase.input_scitype_union(::Type{<:SVM}) = MLJBase.Continuous
MLJBase.input_is_multivariate(::Type{<:SVM}) = true
MLJBase.target_scitype_union(::Type{<:SVMC}) = MLJBase.Finite
MLJBase.target_scitype_union(::Type{<:SVMR}) = MLJBase.Continuous
MLJBase.input_is_multivariate(::Type{<:SVM}) = true

MLJBase.load_path(::Type{<:SCElasticNet}) = "MLJModels.ScikitLearn_.SCElasticNet"
MLJBase.package_name(::Type{<:SCElasticNet}) = "ScikitLearn"
MLJBase.package_uuid(::Type{<:SCElasticNet}) = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
MLJBase.is_pure_julia(::Type{<:SCElasticNet}) = false
MLJBase.package_url(::Type{<:SCElasticNet}) = "https://github.com/cstjean/ScikitLearn.jl"
MLJBase.input_scitype_union(::Type{<:SCElasticNet}) = MLJBase.Continuous
MLJBase.target_scitype_union(::Type{<:SCElasticNet}) = MLJBase.Continuous
MLJBase.input_is_multivariate(::Type{<:SCElasticNet}) = true

MLJBase.load_path(::Type{<:SCElasticNetCV}) = "MLJModels.ScikitLearn_.SCElasticNet"
MLJBase.package_name(::Type{<:SCElasticNetCV}) = "ScikitLearn"
MLJBase.package_uuid(::Type{<:SCElasticNetCV}) = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
MLJBase.is_pure_julia(::Type{<:SCElasticNetCV}) = false
MLJBase.package_url(::Type{<:SCElasticNetCV}) = "https://github.com/cstjean/ScikitLearn.jl"
MLJBase.input_scitype_union(::Type{<:SCElasticNetCV}) = MLJBase.Continuous
MLJBase.target_scitype_union(::Type{<:SCElasticNetCV}) = MLJBase.Continuous
MLJBase.input_is_multivariate(::Type{<:SCElasticNetCV}) = true
end # module


text:
using Flux: glorot_uniform

"""
Build a model of the architecture introduced in the Xie and Grossman 2018 paper: https://arxiv.org/abs/1710.10324

Input to the resulting model is a FeaturedGraph with feature matrix with `input_feature_length` rows and one column for each node in the input graph.

Network has convolution layers, then pooling to some fixed length, followed by Dense layers leading to output.

# Arguments
- `input_feature_length::Integer`: length of feature vector at each node
- `num_conv::Integer`: number of convolutional layers
- `conv_activation::F`: activation function on convolutional layers
- `atom_conv_feature_length::Integer`: length of output of conv layers
- `pool_type::String`: type of pooling after convolution (mean or max)
- `pool_width::Float`: fraction of atom_conv_feature_length that pooling window should span
- `pooled_feature_length::Integer`: feature length to pool down to
- `num_hidden_layers::Integer`: how many Dense layers before output? Note that if this is set to 1 there will be no nonlinearity imposed on these layers
- `hidden_layer_activation::F`: activation function on hidden layers
- `output_layer_activation::F`: activation function on output layer; should generally be identity for regression and something that normalizes appropriately (e.g. softmax) for classification
- `output_length::Integer`: length of output vector
- `initW::F`: function to use to initialize weights in trainable layers
"""
function build_CGCNN(
    input_feature_length;
    num_conv = 2,
    conv_activation = softplus,
    atom_conv_feature_length = 80,
    pool_type = "mean",
    pool_width = 0.1,
    pooled_feature_length = 40,
    num_hidden_layers = 1,
    hidden_layer_activation = softplus,
    output_layer_activation = identity,
    output_length = 1,
    initW = glorot_uniform,
)
    @assert atom_conv_feature_length >= pooled_feature_length "Feature length after pooling must be <= feature length before pooling!"
    model = Chain(
        AGNConv(
            input_feature_length => atom_conv_feature_length,
            conv_activation,
            initW = initW,
        ),
        [
            AGNConv(
                atom_conv_feature_length => atom_conv_feature_length,
                conv_activation,
                initW = initW,
            ) for i = 1:num_conv-1
        ]...,
        AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width),
        [
            Dense(
                pooled_feature_length,
                pooled_feature_length,
                hidden_layer_activation,
                init = initW,
            ) for i = 1:num_hidden_layers-1
        ]...,
        Dense(pooled_feature_length, output_length, output_layer_activation, init = initW),
    )
end

# these dispatches can be removed once https://github.com/FluxML/Flux.jl/issues/1673 is fixed...
# (m::Parallel)(x) = mapreduce(f -> f(x), m.connection, m.layers)
# (m::Parallel)(xs...) = mapreduce((f, x) -> f(x), m.connection, m.layers, xs)
"""
This is a helper function to the main model builder below. Takes in the inputs and models for the two "parallel" CGCNN-like models at the start of the SGCNN architecture and outputs the concatenated final result.
"""
function slab_graph_layer(
    bulk_graph::FeaturizedAtoms{AtomGraph,GraphNodeFeaturization},
    bulk_model,
    surface_graph::FeaturizedAtoms{AtomGraph,GraphNodeFeaturization},
    surface_model,
)
    bulk_output = bulk_model(bulk_graph)
    surface_output = surface_model(surface_graph)
    vcat(bulk_output, surface_output)
end

"""
Build a slab graph network, based off of Kim et al. 2020: https://pubs.acs.org/doi/10.1021/acs.chemmater.9b03686

For now, fixing both "parallel" convolutional paths to have same hyperparams for simplicity and to basically be copies of CGCNN. Might relax this later.

# Arguments:
Same as [`build_CGCNN`](@ref) except for additional parameter of `hidden_layer_width`
"""
function build_SGCNN(
    input_feature_length::Integer;
    num_conv = 2,
    conv_activation = softplus,
    atom_conv_feature_length = 80,
    pool_type = "mean",
    pool_width = 0.1,
    pooled_feature_length = 40,
    hidden_layer_width = 40,
    num_hidden_layers = 3,
    hidden_layer_activation = softplus,
    output_layer_activation = identity,
    output_length = 1,
    initW = glorot_uniform,
)
    bulk_model = Chain(
        AGNConv(
            input_feature_length => atom_conv_feature_length,
            conv_activation,
            initW = initW,
        ),
        [
            AGNConv(
                atom_conv_feature_length => atom_conv_feature_length,
                conv_activation,
                initW = initW,
            ) for i = 1:num_conv-1
        ]...,
        AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width),
    )
    surface_model = Chain(
        AGNConv(
            input_feature_length => atom_conv_feature_length,
            conv_activation,
            initW = initW,
        ),
        [
            AGNConv(
                atom_conv_feature_length => atom_conv_feature_length,
                conv_activation,
                initW = initW,
            ) for i = 1:num_conv-1
        ]...,
        AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width),
    )
    model = Chain(
        Parallel(vcat, surface_model, bulk_model),
        Dense(
            2 * pooled_feature_length,
            hidden_layer_width,
            hidden_layer_activation,
            init = initW,
        ),
        [
            Dense(
                hidden_layer_width,
                hidden_layer_width,
                hidden_layer_activation,
                init = initW,
            ) for i = 1:num_hidden_layers-1
        ]...,
        Dense(hidden_layer_width, output_length, output_layer_activation, init = initW),
    )
end

"""
Similar to [`build_CGCNN`](@ref) but in the style of Deep Equilibrium Models: https://arxiv.org/abs/1909.01377

Can be thought of as equivalent to a CGCNN with infinite convolutional layers.

Input to the resulting model is a FeaturedGraph with feature matrix with `input_feature_length` rows and one column for each node in the input graph.

# Arguments
- `input_feature_length::Integer`: length of feature vector at each node
- `conv_activation::F`: activation function on convolutional layers
- `pool_type::String`: type of pooling after convolution (mean or max)
- `pool_width::Float`: fraction of atom_conv_feature_length that pooling window should span
- `pooled_feature_length::Integer`: feature length to pool down to
- `num_hidden_layers::Integer`: how many Dense layers before output? Note that if this is set to 1 there will be no nonlinearity imposed on these layers
- `hidden_layer_activation::F`: activation function on hidden layers
- `output_layer_activation::F`: activation function on output layer; should generally be identity for regression and something that normalizes appropriately (e.g. softmax) for classification
- `output_length::Integer`: length of output vector
- `initW::F`: function to use to initialize weights in trainable layers
"""
function build_CGCNN_DEQ(
    input_feature_length;
    conv_activation = softplus,
    pool_type = "mean",
    pool_width = 0.1,
    pooled_feature_length = 40,
    num_hidden_layers = 1,
    hidden_layer_activation = softplus,
    output_layer_activation = identity,
    output_length = 1,
    initW = glorot_uniform,
)
    @assert input_feature_length >= pooled_feature_length "Feature length after pooling must be <= feature length before pooling!"
    model = Chain(
        AGNConvDEQ(input_feature_length, conv_activation, initW = initW),
        AGNPool(pool_type, input_feature_length, pooled_feature_length, pool_width),
        [
            Dense(
                pooled_feature_length,
                pooled_feature_length,
                hidden_layer_activation,
                init = initW,
            ) for i = 1:num_hidden_layers-1
        ]...,
        Dense(pooled_feature_length, output_length, output_layer_activation, init = initW),
    )
end


text:
using Pkg

org = ENV["ORG"]
project = ENV["PROJECT"]
branch = ENV["BRANCH"]
entrypoint = ENV["ENTRYPOINT"]

println("ORG: $org")
println("PROJECT: $project")
println("BRANCH: $branch")
println("ENTRYPOINT: $entrypoint")

package = replace(project, ".jl" => "")
println("PACKAGE: $package")

repo="git@github.com:$org/$project.git"
println("REPO: $repo")

println("Setup TimeZones")
Pkg.add("TimeZones")
Pkg.build("TimeZones")

println("Git checkout")
mkpath("/root/.julia/dev/")
run(`git clone --branch=$branch $repo /root/.julia/dev/$package`)

println("Creating symlink: ", "/root/.julia/dev/$package/$entrypoint", "/run.sh")
symlink("/root/.julia/dev/$package/$entrypoint", "/run.sh")

println("Activating package: $package")
Pkg.activate("/root/.julia/dev/$package")
println("Instantiating package: $package")
Pkg.instantiate()
println("Using package: $package")
@eval using $(Symbol(package))

max_stars_repo_name: chipkent/action-julia-docker

max_issues_repo_name: chipkent/action-julia-docker

max_forks_repo_name: chipkent/action-julia-docker


text:
@testset "Operators under mixed barrier conditions" begin
    @testset "1-element xs" begin
        # right at the elements at xs
        @test findnearestindex([1], 1) == 1

        # values out of the range
        @test findnearestindex([1], 0.8) == 1
        @test findnearestindex([1], 9) == 1
    end

    @testset "2-element xs" begin
        # right at the elements at xs
        @test findnearestindex([1;3], 1) == 1
        @test findnearestindex([1;3], 3) == 2

        # values in between
        @test findnearestindex([1;3], 1.2) == 1
        @test findnearestindex([1;3], 1.5) == 1
        @test findnearestindex([1;3], 1.8) == 1
        @test findnearestindex([1;3], 2) == 2
        @test findnearestindex([1;3], 2.0) == 2
        @test findnearestindex([1;3], 2.3) == 2
        @test findnearestindex([1;3], 2.8) == 2

        # values out of the range
        @test findnearestindex([1;3], 0.8) == 1
        @test findnearestindex([1;3], 9) == 2
    end

    @testset "3-element xs" begin
        # right at the elements at xs
        @test findnearestindex([1;2;3], 1) == 1
        @test findnearestindex([1;2;3], 2) == 2
        @test findnearestindex([1;2;3], 3) == 3

        # values in between
        @test findnearestindex([1;2;3], 1.2) == 1
        @test findnearestindex([1;2;3], 1.5) == 2
        @test findnearestindex([1;2;3], 1.8) == 2
        @test findnearestindex([1;2;3], 2.3) == 2
        @test findnearestindex([1;2;3], 2.8) == 3
        @test findnearestindex([1;2;3], 3.3) == 3

        # values out of the range
        @test findnearestindex([1;2;3], 0.8) == 1
        @test findnearestindex([1;2;3], 9) == 3
    end
end


text:
@testset "attention.jl" begin
    @testset "Attention" begin
        dim_x = 2
        dim_y = 3
        dim_embedding = 20
        num_heads = 5
        batch_size = 6
        n = 7
        m = 8

        layer = attention(
            dim_x=dim_x,
            dim_y=dim_y,
            dim_embedding=dim_embedding,
            num_heads=num_heads
        )

        xc = randn(Float32, n, dim_x, batch_size)
        yc = randn(Float32, n, dim_y, batch_size)
        xt = randn(Float32, m, dim_x, batch_size)

        # Perform encodings.
        keys = layer.encoder_x(xc)
        queries = layer.encoder_x(xt)
        values = layer.encoder_xy(cat(xc, yc, dims=2))

        dim_head = div(dim_embedding, num_heads)

        # Compute with layer.
        _, out = code(layer, xc, yc, xt)

        # Brute-force the attention computation.
        embeddings = zeros(Float32, m, dim_head, num_heads, batch_size)
        for c = 1:num_heads
            for b = 1:batch_size
                # Calculate weights.
                weights = Array{Float32}(undef, n, m)
                for i = 1:n, j = 1:m
                    key = keys[i, :, c, b]
                    query = queries[j, :, c, b]
                    # Normalise to control variance of product.
                    product = dot(key, query) / Float32(sqrt(length(key)))
                    weights[i, j] = exp(product)
                end
                for j = 1:m
                    weights[:, j] ./= sum(weights[:, j])
                end

                # Calculate embeddings.
                for i = 1:n, j = 1:m
                    embeddings[j, :, c, b] .+= values[i, :, c, b] .* weights[i, j]
                end
            end
        end
        ref = layer.transformer(layer.mixer(embeddings), queries)

        # Check that the layer lines up with the brute-force reference.
        @test out ≈ ref
    end

    @testset "BatchedMLP" begin
        layer = NeuralProcesses.batched_mlp(
            dim_in=2,
            dim_out=3,
            dim_hidden=10,
            num_layers=3
        )
        x = randn(10, 2, 4, 5)
        @test size(layer(x)) == (10, 3, 4, 5)
    end
end


text:
function LJForceMagnitude(r::T,d::T) where T<:Real
    if r > d || r < 0.8*d
        return zero(T)
    else
        return 2*d/r^2*(1 - d/r)
    end
end

function ForceCalculation(d::T,Δx::T,Δy::T) where T<:Real
    r = sqrt(Δx^2 + Δy^2)
    tx = Δx/r
    ty = Δy/r

    LJmag = LJForceMagnitude(r,d)

    Fx = tx*LJmag
    Fy = ty*LJmag

    return Fx,Fy
end

function compute_body_repulsion(ind::Int,olist)
    FX = 0.0
    FY = 0.0

    for ti=1:(ind-1)
        Δx = olist[ind].center.x - olist[ti].center.x
        Δy = olist[ind].center.y - olist[ti].center.y
        d = olist[ti].radius + olist[ind].radius +
            (1 + olist[ti].buffer_percent + olist[ind].buffer_percent)/100

        fx,fy = ForceCalculation(d,Δx,Δy)
        FX += fx; FY += fy
    end

    return FX, FY
end

function compute_wall_repulsion(ind,olist,wlist)
    FX = 0.0
    FY = 0.0

    for tw=1:length(wlist)
        pointOnWall = get_nearest_point(olist[ind].center,wlist[tw])
        if is_point_in_wall(pointOnWall,wlist[tw])
            Δx = pointOnWall.x - olist[ind].center.x
            Δy = pointOnWall.y - olist[ind].center.y

            d = olist[ind].radius + wlist[tw].thickness/2

            fx,fy = ForceCalculation(d,Δx,Δy)
            FX += fx
            FY += fy
        end
    end

    return FX, FY
end

function shuffle_object!(ind::Int,olist,wlist)
    fXr,fYr = compute_body_repulsion(ind,olist)
    fXw,fYw = compute_wall_repulsion(ind,olist,wlist)
    olist[ind].center.x += (fXr + fXw)*1e-6
    olist[ind].center.y += (fYr + fYw)*1e-6
    nothing
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
using LightDarkPOMDPs
using POMDPToolbox
using POMCP
using Plots
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using ParticleFilters
using Reel
using JLD

pomdp = LightDark2DTarget(term_radius=0.1,
                    init_dist=SymmetricNormal2([2.0, 2.0], 5.0),
                    discount=0.95)
titles = []
hists = []

rng_seed = 10
rng = MersenneTwister(rng_seed)

## VANILLA ##

rng3 = copy(rng)
action_gen = AdaptiveRadiusRandom(max_radius=10.0, to_zero_first=true, rng=rng3)
rollout_policy = RadiusRandom(radius=10.0, rng=rng3)
tree_queries = 50_000
updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(1000), 0.05, rng)
solver = POMCPDPWSolver(next_action=action_gen,
                        tree_queries=tree_queries,
                        c=20.0,
                        max_depth=40,
                        k_action=10.0,
                        alpha_action=1/8,
                        k_observation=4.0,
                        alpha_observation=1/8,
                        estimate_value=RolloutEstimator(rollout_policy),
                        rng=rng3
                       )
policy = solve(solver, pomdp)

#=
ib = initial_state_distribution(pomdp)
node = RootNode(ib)
a = action(policy, node)
POMCP.blink(node)
=#

hr = HistoryRecorder(max_steps=40, rng=MersenneTwister(rng_seed), initial_state=Vec2(2.0, 3.0), show_progress=true)
hist = simulate(hr, pomdp, policy, updater)
@show n_steps(hist)
println(discounted_reward(hist))
push!(titles, "Vanilla_POMCP")
push!(hists, hist)

## WITH FEEDBACK ##

#=
rng3 = copy(rng)
action_gen = AdaptiveRadiusRandom(max_radius=10.0, to_zero_first=true, rng=rng3)
rollout_policy = SimpleFeedback(max_radius=10.0)
rollout_updater = FastPreviousObservationUpdater{Vec2}()
tree_queries = 50_000
updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(1000), 0.05, rng)
solver = POMCPDPWSolver(next_action=action_gen,
                        tree_queries=tree_queries,
                        c=20.0,
                        max_depth=40,
                        k_action=10.0,
                        alpha_action=1/8,
                        k_observation=4.0,
                        alpha_observation=1/8,
                        estimate_value=PORolloutEstimator(rollout_policy, rollout_updater),
                        rng=rng3
                       )
policy = solve(solver, pomdp)

#=
ib = initial_state_distribution(pomdp)
node = RootNode(ib)
a = action(policy, node)
POMCP.blink(node)
=#

hr = HistoryRecorder(max_steps=40, rng=MersenneTwister(rng_seed), initial_state=Vec2(2.0, 3.0), show_progress=true)
hist = simulate(hr, pomdp, policy, updater)
println(discounted_reward(hist))
# write_gif(pomdp, hist, "")
=#

## WITH BELIEF ##

rng3 = copy(rng)
action_gen = AdaptiveRadiusRandom(max_radius=10.0, to_zero_first=true, rng=rng3)
rollout_policy = SimpleFeedback(max_radius=10.0)
tree_queries = 50_000
updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(1000), 0.05, rng)
node_updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100), 0.05, rng)
solver = POMCPDPWSolver(next_action=action_gen,
                        tree_queries=tree_queries,
                        c=20,
                        max_depth=40,
                        k_action=10.0,
                        alpha_action=1/8,
                        k_observation=4.0,
                        alpha_observation=1/8,
                        estimate_value=PORolloutEstimator(rollout_policy, node_updater),
                        node_belief_updater=node_updater,
                        rng=rng3
                       )
policy = solve(solver, pomdp)

#=
ib = initial_state_distribution(pomdp)
node = RootNode(ib)
a = action(policy, node)
POMCP.blink(node)
=#

hr = HistoryRecorder(max_steps=40, rng=MersenneTwister(rng_seed), initial_state=Vec2(2.0, 3.0), show_progress=true)
hist = simulate(hr, pomdp, policy, updater)
@show n_steps(hist)
println(discounted_reward(hist))
push!(titles, "Modified_POMCP")
push!(hists, hist)

@save(joinpath(tempdir(), "gifplots.jld"), pomdp, titles, hists)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@everywhere module TestDecisionConstraint

using StochasticPrograms
using Test

function test_SingleDecision_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @decision(sp, z[1:2])
        @variable(sp, t)
        @constraint(sp, con11, x in MOI.LessThan(10.0))
        @constraint(sp, con12[i in 1:2], z[i] in MOI.LessThan(float(i)))
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, r)
        @recourse(sp, y)
        @recourse(sp, w[1:2])
        @constraint(sp, con21, y in MOI.LessThan(a))
        @constraint(sp, con22[i in 1:2], w[i] in MOI.LessThan(float(i)))
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    z = DecisionRef.(sp[1,:z])
    con11 = sp[1,:con11]
    @test "con11" == @inferred JuMP.name(con11)
    @test index(con11) == JuMP.constraint_by_name(owner_model(con11), 1, "con11").index
    c = JuMP.constraint_object(con11)
    @test c.func == x
    @test c.set == MOI.LessThan(10.0)
    con12 = sp[1,:con12]
    @test "con12[1]" == @inferred JuMP.name(con12[1])
    @test index(con12[1]) == JuMP.constraint_by_name(owner_model(con12[1]), 1, "con12[1]").index
    c = JuMP.constraint_object(con12[1])
    @test c.func == z[1]
    @test c.set == MOI.LessThan(1.0)
    # Scenario-dependent
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    w1 = DecisionRef.(sp[2,:w], 1)
    w2 = DecisionRef.(sp[2,:w], 2)
    con21 = sp[2,:con21]
    @test "con21" == @inferred JuMP.name(con21)
    @test index(con21) == JuMP.constraint_by_name(owner_model(con21), 2, "con21").index
    c1 = JuMP.constraint_object(con21, 1)
    c2 = JuMP.constraint_object(con21, 2)
    @test c1.func == y1
    @test c2.func == y2
    @test c1.set == MOI.LessThan(1.0)
    @test c2.set == MOI.LessThan(2.0)
    con22 = sp[2,:con22]
    @test "con22[1]" == @inferred JuMP.name(con22[1])
    @test index(con22[1]) == JuMP.constraint_by_name(owner_model(con22[1]), 2, "con22[1]").index
    c1 = JuMP.constraint_object(con22[1], 1)
    c2 = JuMP.constraint_object(con22[1], 2)
    @test c1.func == w1[1]
    @test c2.func == w2[1]
    @test c.set == MOI.LessThan(1.0)
end

function test_VectorOfDecisions_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x[1:2])
        @variable(sp, w)
        @constraint(sp, con11, x in MOI.Zeros(2))
        @constraint(sp, con12, [x[2],x[1]] in MOI.Zeros(2))
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y[1:2])
        @constraint(sp, con21, y in MOI.Zeros(2))
        @constraint(sp, con22, [y[2],y[1]] in MOI.Zeros(2))
    end
    # First-stage
    x = DecisionRef.(sp[1,:x])
    sp_cref = sp[1,:con11]
    c = JuMP.constraint_object(sp_cref)
    @test c.func == x
    @test c.set == MOI.Zeros(2)
    sp_cref = sp[1,:con12]
    c = JuMP.constraint_object(sp_cref)
    @test c.func == [x[2],x[1]]
    @test c.set == MOI.Zeros(2)
    # Scenario-dependent
    y = DecisionRef.(sp[2,:y], 2)
    sp_cref = sp[2,:con21]
    c = JuMP.constraint_object(sp_cref, 2)
    @test c.func == y
    @test c.set == MOI.Zeros(2)
    sp_cref = sp[2,:con22]
    c = JuMP.constraint_object(sp_cref, 2)
    @test c.func == [y[2],y[1]]
    @test c.set == MOI.Zeros(2)
end

function test_DecisionAffExpr_scalar_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @constraint(sp, con11, x >= 1.)
        @constraint(sp, con12, 2x <= 10)
        @constraint(sp, con13, 3x + 1 >= 10)
        @constraint(sp, con14, 1 == -x)
        @constraint(sp, con15, 2 == 1)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @constraint(sp, con21, y >= a)
        @constraint(sp, con22, 2x + a*y <= 10)
        @constraint(sp, con23, a*x + 3y + 1 >= 10)
        @constraint(sp, con24, 1 == -y - x)
        @constraint(sp, con25, 2 == 1)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    c = JuMP.constraint_object(sp[1,:con11])
    @test JuMP.isequal_canonical(c.func, 1.0x)
    @test c.set == MOI.GreaterThan(1.0)
    c = JuMP.constraint_object(sp[1,:con12])
    @test JuMP.isequal_canonical(c.func, 2x)
    @test c.set == MOI.LessThan(10.0)
    c = JuMP.constraint_object(sp[1,:con13])
    @test JuMP.isequal_canonical(c.func, 3x)
    @test c.set == MOI.GreaterThan(9.0)
    c = JuMP.constraint_object(sp[1,:con14])
    @test JuMP.isequal_canonical(c.func, 1.0x)
    @test c.set == MOI.EqualTo(-1.0)
    @test_throws ErrorException JuMP.constraint_object(sp[1,:con15])
    # Scenario-dependent
    x = DecisionRef(sp[1,:x], 2, 1)
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    c1 = JuMP.constraint_object(sp[2,:con21], 1)
    c2 = JuMP.constraint_object(sp[2,:con21], 2)
    @test JuMP.isequal_canonical(c1.func, 1.0y1)
    @test JuMP.isequal_canonical(c2.func, 1.0y2)
    @test c1.set == MOI.GreaterThan(1.0)
    @test c2.set == MOI.GreaterThan(2.0)
    c1 = JuMP.constraint_object(sp[2,:con22], 1)
    c2 = JuMP.constraint_object(sp[2,:con22], 2)
    @test JuMP.isequal_canonical(c1.func, 2x + y1)
    @test JuMP.isequal_canonical(c2.func, 2x + 2y2)
    @test c1.set == MOI.LessThan(10.0)
    @test c2.set == MOI.LessThan(10.0)
    c1 = JuMP.constraint_object(sp[2,:con23], 1)
    c2 = JuMP.constraint_object(sp[2,:con23], 2)
    @test JuMP.isequal_canonical(c1.func, x + 3y1)
    @test JuMP.isequal_canonical(c2.func, 2x + 3y2)
    @test c1.set == MOI.GreaterThan(9.0)
    @test c2.set == MOI.GreaterThan(9.0)
    c1 = JuMP.constraint_object(sp[2,:con24], 1)
    c2 = JuMP.constraint_object(sp[2,:con24], 2)
    @test JuMP.isequal_canonical(c1.func, 1.0x + 1.0y1)
    @test JuMP.isequal_canonical(c2.func, 1.0x + 1.0y2)
    @test c1.set == MOI.EqualTo(-1.0)
    @test c2.set == MOI.EqualTo(-1.0)
    @test_throws ErrorException JuMP.constraint_object(sp[2,:con25], 1)
end

function test_DecisionAffExpr_vectorized_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @constraint(sp, con1, [x, 2x] .== [1-x, 3])
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @constraint(sp, con2, [a*x + y, 2x + a*y] .== [1-x-y, 3])
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    c = JuMP.constraint_object.(sp[1,:con1])
    @test JuMP.isequal_canonical(c[1].func, 2.0x)
    @test c[1].set == MOI.EqualTo(1.0)
    @test JuMP.isequal_canonical(c[2].func, 2.0x)
    @test c[2].set == MOI.EqualTo(3.0)
    # Scenario-dependent
    x = DecisionRef(sp[1,:x], 2, 1)
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    c1 = JuMP.constraint_object.(sp[2,:con2], 1)
    c2 = JuMP.constraint_object.(sp[2,:con2], 2)
    @test JuMP.isequal_canonical(c1[1].func, 2.0x + 2.0y1)
    @test JuMP.isequal_canonical(c2[1].func, 3.0x + 2.0y2)
    @test c1[1].set == MOI.EqualTo(1.0)
    @test c2[1].set == MOI.EqualTo(1.0)
    @test JuMP.isequal_canonical(c1[2].func, 2.0x + y1)
    @test JuMP.isequal_canonical(c2[2].func, 2.0x + 2.0y2)
    @test c1[2].set == MOI.EqualTo(3.0)
    @test c2[2].set == MOI.EqualTo(3.0)
end

function test_delete_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @constraint(sp, con1, 2x <= 1)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @constraint(sp, con2, a*x + 2y <= 1)
    end
    # First-stage
    sp_cref = sp[1,:con1]
    @test JuMP.is_valid(sp, sp_cref)
    JuMP.delete(sp, sp_cref)
    @test !JuMP.is_valid(sp, sp_cref)
    second_sp = StochasticProgram([Scenario()], Structure...)
    @test_throws Exception JuMP.delete(second_sp, sp_cref)
    # Scenario-dependent
    sp_cref = sp[2,:con2]
    @test_throws ErrorException JuMP.is_valid(sp, sp_cref)
    @test JuMP.is_valid(sp, sp_cref, 1)
    @test JuMP.is_valid(sp, sp_cref, 2)
    @test_throws ErrorException JuMP.delete(sp, sp_cref)
    JuMP.delete(sp, sp_cref, 1)
    @test !JuMP.is_valid(sp, sp_cref, 1)
    @test JuMP.is_valid(sp, sp_cref, 2)
    JuMP.delete(sp, sp_cref, 2)
    @test !JuMP.is_valid(sp, sp_cref, 2)
    @first_stage sp = begin
        @decision(sp, x[1:9])
        @variable(sp, w)
        @constraint(sp, con11, sum(x[1:2:9]) <= 3)
        @constraint(sp, con12, sum(x[2:2:8]) <= 2)
        @constraint(sp, con13, sum(x[1:3:9]) <= 1)
    end
    @second_stage sp = begin
        @known(sp, x)
        @variable(sp, z)
        @recourse(sp, y[1:9])
        @constraint(sp, con21, sum(x[1:2:9]) + sum(y[1:2:9]) <= 3)
        @constraint(sp, con22, sum(x[2:2:8]) + sum(y[2:2:8]) <= 2)
        @constraint(sp, con23, sum(x[1:3:9]) + sum(y[1:3:9]) <= 1)
    end
    # First-stage
    cons = all_constraints(sp, 1, DecisionAffExpr{Float64}, MOI.LessThan{Float64})
    @test all(JuMP.is_valid.(sp, cons))
    JuMP.delete(sp, cons[[1, 3]])
    @test all((!JuMP.is_valid).(sp, cons[[1, 3]]))
    @test JuMP.is_valid(sp, cons[2])
    # Scenario-dependent
    cons = all_constraints(sp, 2, DecisionAffExpr{Float64}, MOI.LessThan{Float64})
    @test_throws ErrorException all(JuMP.is_valid.(sp, cons))
    @test all(JuMP.is_valid.(sp, cons, 1))
    @test all(JuMP.is_valid.(sp, cons, 2))
    @test_throws ErrorException JuMP.delete(sp, cons[[1, 3]])
    JuMP.delete(sp, cons[[1, 3]], 1)
    @test all((!JuMP.is_valid).(sp, cons[[1, 3]], 1))
    @test JuMP.is_valid(sp, cons[2], 1)
    @test all(JuMP.is_valid.(sp, cons, 2))
    JuMP.delete(sp, cons[[1, 3]], 2)
    @test all((!JuMP.is_valid).(sp, cons[[1, 3]], 2))
    @test JuMP.is_valid(sp, cons[2], 2)
end

function test_DecisionQuadrExpr_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @constraint(sp, con1, x^2 + x <= 1)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @constraint(sp, con2, y^2 + a*y*x - 1.0 == 0.0)
    end
    # First-stage
    x = DecisionRef(sp[1,:x])
    sp_cref = sp[1,:con1]
    c = JuMP.constraint_object(sp_cref)
    @test JuMP.isequal_canonical(c.func, x^2 + x)
    @test c.set == MOI.LessThan(1.0)
    # Scenario-dependent
    sp_cref = sp[2,:con2]
    x = DecisionRef(sp[1,:x], 2, 1)
    y1 = DecisionRef(sp[2,:y], 1)
    y2 = DecisionRef(sp[2,:y], 2)
    c1 = JuMP.constraint_object(sp_cref, 1)
    c2 = JuMP.constraint_object(sp_cref, 2)
    @test JuMP.isequal_canonical(c1.func, y1^2 + x*y1)
    @test JuMP.isequal_canonical(c2.func, y2^2 + 2*x*y2)
    @test c1.set == MOI.EqualTo(1.0)
    @test c2.set == MOI.EqualTo(1.0)
    # sp_cref = sp[2,:con3]
    # c = JuMP.constraint_object(cref)
    # @test JuMP.isequal_canonical(c.func[1], -1 + 3x^2 - 4x*y + 2x)
    # @test JuMP.isequal_canonical(c.func[2],  1 - 2x^2 + 2x*y - 3y)
    # @test c.set == MOI.SecondOrderCone(2)
end

function test_all_decision_constraints(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x >= 0)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y >= 0)
        @constraint(sp, x + a*y >= 0)
    end
    # First-stage
    x = sp[1,:x]
    @test 1 == @inferred num_constraints(sp, 1, DecisionRef, MOI.GreaterThan{Float64})
    ref = @inferred all_constraints(sp, 1, DecisionRef, MOI.GreaterThan{Float64})
    @test ref == [LowerBoundRef(x)]
    @test 0 == @inferred num_constraints(sp, 1, DecisionAffExpr{Float64},
                                            MOI.GreaterThan{Float64})
    aff_constraints = all_constraints(sp, 1, DecisionAffExpr{Float64},
                                        MOI.GreaterThan{Float64})
    @test isempty(aff_constraints)
    # Scenario-dependent
    y = sp[2,:y]
    ref = @inferred all_constraints(sp, 2, DecisionRef, MOI.GreaterThan{Float64})
    @test ref == [LowerBoundRef(y)]
    @test 1 == @inferred num_constraints(sp, 2, DecisionAffExpr{Float64},
                                            MOI.GreaterThan{Float64})
    aff_constraints = all_constraints(sp, 2, DecisionAffExpr{Float64},
                                      MOI.GreaterThan{Float64})
    x = DecisionRef(x, 2, 1)
    y1 = DecisionRef(y, 1)
    y2 = DecisionRef(y, 2)
    c1 = constraint_object(aff_constraints[1], 1)
    c2 = constraint_object(aff_constraints[1], 2)
    @test JuMP.isequal_canonical(c1.func, x + y1)
    @test JuMP.isequal_canonical(c2.func, x + 2y2)
end

function test_list_of_constraint_types(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x >= 0, Bin)
        @variable(sp, w)
        @constraint(sp, 2x <= 1)
        @constraint(sp, [x, x] in SecondOrderCone())
        @constraint(sp, x^2- x <= 2)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y >= 0, Bin)
        @constraint(sp, 2x + 2y <= a)
        @constraint(sp, [x + y, x] in SecondOrderCone())
        @constraint(sp, x^2+y^2 - x*y <= 1)
    end
    # First-stage
    constraint_types = @inferred list_of_constraint_types(sp, 1)
    @test Set(constraint_types) == Set(
        [(DecisionRef, MOI.ZeroOne),
         (DecisionRef, MOI.GreaterThan{Float64}),
         (DecisionAffExpr{Float64}, MOI.LessThan{Float64}),
         (Vector{DecisionRef}, MOI.SecondOrderCone),
         (DecisionQuadExpr{Float64}, MOI.LessThan{Float64})])
    # First-stage
    constraint_types = @inferred list_of_constraint_types(sp, 1)
    @test Set(constraint_types) == Set(
        [(DecisionRef, MOI.ZeroOne),
         (DecisionRef, MOI.GreaterThan{Float64}),
         (DecisionAffExpr{Float64}, MOI.LessThan{Float64}),
         (Vector{DecisionRef}, MOI.SecondOrderCone),
         (DecisionQuadExpr{Float64}, MOI.LessThan{Float64})])
    # Second-stage
    constraint_types = @inferred list_of_constraint_types(sp, 2)
    @test Set(constraint_types) == Set(
        [(DecisionRef, MOI.ZeroOne),
        (DecisionRef, MOI.GreaterThan{Float64}),
        (DecisionAffExpr{Float64}, MOI.LessThan{Float64}),
        (Vector{DecisionAffExpr{Float64}}, MOI.SecondOrderCone),
        (DecisionQuadExpr{Float64}, MOI.LessThan{Float64})])
end

function test_change_decision_coefficient(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @constraint(sp, con1, 2 * x == -1)
        @constraint(sp, quadcon1, x^2 == 0)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @constraint(sp, con2, 2x + a*y == 0)
        @constraint(sp, quadcon2, y^2 + a*x == 0)
    end
    # First-stage
    x = sp[1,:x]
    sp_cref = sp[1,:con1]
    @test JuMP.normalized_coefficient(sp_cref, x) == 2.0
    JuMP.set_normalized_coefficient(sp_cref, x, 1.0)
    @test JuMP.normalized_coefficient(sp_cref, x) == 1.0
    JuMP.set_normalized_coefficient(sp_cref, x, 3)
    @test JuMP.normalized_coefficient(sp_cref, x) == 3.0
    quad_con = sp[1,:quadcon1]
    @test JuMP.normalized_coefficient(quad_con, x) == 0.0
    JuMP.set_normalized_coefficient(quad_con, x, 2)
    @test JuMP.normalized_coefficient(quad_con, x) == 2.0
    x = DecisionRef(x)
    @test JuMP.isequal_canonical(
        JuMP.constraint_object(quad_con).func, x^2 + 2x)
    # Scenario-dependent
    x = sp[1,:x]
    y = sp[2,:y]
    sp_cref = sp[2,:con2]
    @test_throws ErrorException JuMP.normalized_coefficient(sp_cref, y) == 2.0
    @test JuMP.normalized_coefficient(sp_cref, y, 1) == 1.0
    @test JuMP.normalized_coefficient(sp_cref, y, 2) == 2.0
    @test_throws ErrorException JuMP.set_normalized_coefficient(sp_cref, y, 1.0)
    JuMP.set_normalized_coefficient(sp_cref, y, 1, 2.0)
    JuMP.set_normalized_coefficient(sp_cref, y, 2, 1.0)
    @test JuMP.normalized_coefficient(sp_cref, y, 1) == 2.0
    @test JuMP.normalized_coefficient(sp_cref, y, 2) == 1.0
    JuMP.set_normalized_coefficient(sp_cref, y, 1, 3)
    JuMP.set_normalized_coefficient(sp_cref, y, 2, 3)
    @test JuMP.normalized_coefficient(sp_cref, y, 1) == 3.0
    @test JuMP.normalized_coefficient(sp_cref, y, 2) == 3.0
    quad_con = sp[2,:quadcon2]
    @test_throws ErrorException JuMP.normalized_coefficient(quad_con, y) == 0.0
    @test JuMP.normalized_coefficient(quad_con, y, 1) == 0.0
    @test JuMP.normalized_coefficient(quad_con, y, 2) == 0.0
    @test_throws ErrorException JuMP.set_normalized_coefficient(quad_con, y, 2)
    JuMP.set_normalized_coefficient(quad_con, y, 1, 2)
    JuMP.set_normalized_coefficient(quad_con, y, 2, 2)
    @test JuMP.normalized_coefficient(quad_con, y, 1) == 2.0
    @test JuMP.normalized_coefficient(quad_con, y, 2) == 2.0
    x = DecisionRef(x, 2, 1)
    y1 = DecisionRef(y, 1)
    y2 = DecisionRef(y, 2)
    @test JuMP.isequal_canonical(
        JuMP.constraint_object(quad_con, 1).func, y1^2 + 2y1 + x)
    @test JuMP.isequal_canonical(
        JuMP.constraint_object(quad_con, 2).func, y2^2 + 2y2 + 2x)
end

function test_change_decision_rhs(Structure)
    ξ₁ = @scenario a = 1. probability = 0.5
    ξ₂ = @scenario a = 2. probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
        @constraint(sp, con1, 2 * x <= 1)
    end
    @second_stage sp = begin
        @known(sp, x)
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y)
        @constraint(sp, con2, 2x + y <= a)
    end
    # First-stage
    sp_cref = sp[1,:con1]
    @test JuMP.normalized_rhs(sp_cref) == 1.0
    JuMP.set_normalized_rhs(sp_cref, 2.0)
    @test JuMP.normalized_rhs(sp_cref) == 2.0
    JuMP.set_normalized_rhs(sp_cref, 3)
    @test JuMP.normalized_rhs(sp_cref) == 3.0
    # Scenario-dependent
    sp_cref = sp[2,:con2]
    @test_throws ErrorException JuMP.normalized_rhs(sp_cref)
    @test JuMP.normalized_rhs(sp_cref, 1) == 1.0
    @test JuMP.normalized_rhs(sp_cref, 2) == 2.0
    @test_throws ErrorException JuMP.set_normalized_rhs(sp_cref, 2.0)
    JuMP.set_normalized_rhs(sp_cref, 1, 2.0)
    JuMP.set_normalized_rhs(sp_cref, 2, 1.0)
    @test JuMP.normalized_rhs(sp_cref, 1) == 2.0
    @test JuMP.normalized_rhs(sp_cref, 2) == 1.0
    JuMP.set_normalized_rhs(sp_cref, 1, 3)
    JuMP.set_normalized_rhs(sp_cref, 2, 3)
    @test JuMP.normalized_rhs(sp_cref, 1) == 3.0
    @test JuMP.normalized_rhs(sp_cref, 2) == 3.0
end

function runtests()
    @testset "DecisionConstraint" begin
        for structure in [(Deterministic(),),
                          (StageDecomposition(),),
                          (ScenarioDecomposition(),),
                          (Deterministic(), () -> MOIU.MockOptimizer(MOIU.Model{Float64}()))]
            @testset "$(structure)" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

function run_dtests()
    @testset "DecisionConstraint" begin
        for structure in [(DistributedStageDecomposition(),), (DistributedScenarioDecomposition(),)]
            @testset "$(structure)" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

end


text:
struct StringStrategy <: SearchStrategy
    alphabet::Vector{Char}
    minsize::UInt32
    maxsize::UInt32

    StringStrategy(alphabet, minsize, maxsize) = new(alphabet, minsize, maxsize)
end

function draw(strategy::StringStrategy)
    length = rand(strategy.minsize:strategy.maxsize)
    return string(rand(strategy.alphabet, length)...)
end

Base.eltype(::StringStrategy) = String

strings(; minsize=0, maxsize=5000) = strings(:alpha; minsize=minsize, maxsize=maxsize)
strings(sym::Symbol; minsize=0, maxsize=5000) = strings(Val(sym); minsize=minsize, maxsize=maxsize)
strings(::Val{:alpha}; minsize=0, maxsize=5000) = StringStrategy(collect(Iterators.flatten(('A':'Z', 'a':'z'))), minsize, maxsize)
strings(::Val{:alphanum}; minsize=0, maxsize=5000) = StringStrategy(collect(Iterators.flatten(('0':'9', 'A':'Z', 'a':'z'))), minsize, maxsize)
strings(alphabet; minsize=0, maxsize=5000) = StringStrategy(alphabet, minsize, maxsize)

struct SymbolStrategy <: SearchStrategy
    alphabet::Vector{Char}
    minsize::UInt32
    maxsize::UInt32

    SymbolStrategy(alphabet, minsize, maxsize) = new(alphabet, minsize, maxsize)
end

function draw(strategy::SymbolStrategy)
    length = rand(strategy.minsize:strategy.maxsize)
    return Symbol(rand(strategy.alphabet, length)...)
end

Base.eltype(::SymbolStrategy) = Symbol

symbols(; minsize=0, maxsize=5000) = symbols(:alpha; minsize=minsize, maxsize=maxsize)
symbols(sym::Symbol; minsize=0, maxsize=5000) = symbols(Val(sym); minsize=minsize, maxsize=maxsize)
symbols(::Val{:alpha}; minsize=0, maxsize=5000) = SymbolStrategy(collect(Iterators.flatten(('A':'Z', 'a':'z'))), minsize, maxsize)
symbols(::Val{:alphanum}; minsize=0, maxsize=5000) = SymbolStrategy(collect(Iterators.flatten(('0':'9', 'A':'Z', 'a':'z'))), minsize, maxsize)
symbols(alphabet; minsize=0, maxsize=5000) = SymbolStrategy(alphabet, minsize, maxsize)


text:
function get_valid_log(instance_name, outpath=joinpath("data", "mosek_logs", instance_name))
    index =  0
    file =  "$instance_name.log"
    while ispath(joinpath(outpath, file))
        index += 1
        file = "$(instance_name)_$(index).log"
    end
    return file
end

"""
    solve_sdp(instance_name, cliques, cliquetree, path_opf_ctr, path_opf_mat; time_limit=100000.0)

Solve an SDP cliques decompositions.
"""
function solve_sdp(instance_name, cliques, cliquetree, path_opf_ctr, path_opf_mat; time_limit=100000.0)
    REPO = "OPF"
    objscale = 4 #0 if no scaling
    originalSTDOUT = stdout
    outpath = joinpath("data/mosek_logs", instance_name)
    isdir(outpath) || mkpath(outpath)
    log_file = get_valid_log(instance_name, outpath)
    outlog = open(joinpath(outpath,log_file), "w")
    redirect_stdout(outlog)
    println("objscale : 10^(-$objscale) \n")
    Xoptimal, obj, status, time, mem = sdp_model_solver_to_specify(cliques, cliquetree, path_opf_ctr, path_opf_mat, 10.0^(-objscale); time_limit=time_limit)
    println("time : $time \n memory : $mem")
    close(outlog)
    redirect_stdout(originalSTDOUT)
    nb_lc, m, time, objective, nb_iter = read_mosek_log(joinpath(outpath, log_file))
    #time, nb_iter, objective, m, nlc = read_mosek_log(joinpath(outpath, log_file))
    return time, nb_iter, joinpath(outpath, log_file), objective, m, nb_lc
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using DataFrames, Random, WeakRefStrings
using JDF

@testset "JDF.jl nrow" begin
    df = DataFrame([collect(1:100) for i = 1:3000])
    df[!, :int_missing] =
        rand([rand(rand([UInt, Int, Float64, Float32, Bool])), missing], nrow(df))
    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i = 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i = 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i = 1:nrow(df)]
    df[!, :stringarray_missing] =
        StringVector([rand([missing, randstring(8)]) for i = 1:nrow(df)])

    JDF.save("a.jdf", df)

    dfjdf = jdf"a.jdf"

    @test ncol(jdf"a.jdf") == 3006
    @test nrow(jdf"a.jdf") == 100
    @test size(jdf"a.jdf") == (100, 3006)
    @test size(jdf"a.jdf", 1) == 100
    @test size(jdf"a.jdf", 2) == 3006

    # clean up
    rm("a.jdf", force = true, recursive = true)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
module cartoonLike

using Plots
using LinearAlgebra
using Interpolations
using NearestNeighbors

export create_circle, given_points, distinguishing_curvature, cartoonLikeImage

# still to make: piecewise smooth curves
#                fill in the interior, exterior with something smooth

struct cartoonLikeImage{T}
  xSize::Int
  ySize::Int
  ninterp::Int
  ε::T
  function cartoonLikeImage{T}(; xSize=100, ySize=100, ninterp = xSize*ySize,
                               ε = 1/(xSize*ySize)^(3/4)) where T<: Real
    return new(xSize, ySize, ninterp, T(ε))
  end
end

"""
    create_circle(r, rate = 10)
create a circle of radius r, sampled at rate `rate`
"""
function create_circle(r, rate = 10)
    t = range(0,1,length = 10)
    x = r .* sin.(2π*t)
    y = r .* cos.(2π*t)
    A = hcat(x,y)
end

"""
    image, xlocs, ylocs = given_points(points; xSize=100, ySize=100, ninterp =
        xSize*ySize, ε = .01)

given a collection of nodes, i.e. row vectors of the form (x,y), make a
(xSize)×(ySize) image of path going through them, which has been blurred with a
gaussian of standard deviation of ε (\epsilon)
"""
function given_points(points::AbstractArray{T}; cart::cartoonLikeImage{<:Real} =
                      cartoonLikeImage{T}()) where T<:Real
    xSize = cart.xSize
    ySize = cart.ySize
    ninterp = cart.ninterp
    ε = cart.ε
    t= range(0, 1, length = size(points, 1))
    itp = scale(interpolate(points, (BSpline(Cubic(Natural(OnGrid()))), NoInterp())), t, 1:2)
    tfine = range(0, 1, length = ceil(Int, ninterp))
    finePoints = itp(tfine, 1:2)'
    tree = KDTree(Array{Float64,2}(finePoints))

    xRange = range(-1, 1, length=xSize)
    yRange = range(-1, 1, length=ySize)

    imagePoints = zeros(2, xSize, ySize)
    for (ix,x) =enumerate(xRange), (iy,y) = enumerate(yRange)
        imagePoints[1, ix, iy] = x
        imagePoints[2, ix, iy] = y
    end
    image = zeros(xSize, ySize)
    idxs, dists = knn(tree, reshape(imagePoints,(2,:)),1)
    dists = [d[1] for d in dists]       # the distances are given as lists, and we want
    # it as just an array
    dists = reshape(dists, size(image))
    image = exp.(-(dists).^2 ./ε)
    return (image, xRange, yRange)
end

"""
make a collection of `n` images that have a range of different curvatures at a
  central point, equally spaced
"""
function distinguishing_curvature(nPoints; cart::cartoonLikeImage{<:Real} =
                                  cartoonLikeImage{Float32}())
    yRange = range(-1,1,length=2*nPoints+1)
    points = make_points.(yRange)
    image_groups = given_points.(points, cart = cart)
    xLocs = image_groups[1][2]

    yLocs = image_groups[1][3]

    images = zeros(cart.xSize, cart.ySize, length(image_groups))
    for (i,tup) in enumerate(image_groups)
        images[:, :, i] = tup[1]
    end
    return images
end

function make_points(y)
  return [-.5 0 .5;
           y  0  y]'
end

end


text:
# # Absorption: Spectral Line Shapes

# Part of ESE 156 Remote Sensing Class; 2020; Christian Frankenberg

# The purpose of this tutorial is to demonstrate how to compute line-shapes (in our case mostly of rotational-vibrational transition lines), what processes determine line-shapes and how the dependencies in the line-shape on pressure and temperature could actually be used to our advantage in a retrieval system.

# ### Basic Tools

# We are making use of the [HITRAN database](https://www.cfa.harvard.edu/hitran/), see a list of tutorials [here](https://www.youtube.com/playlist?list=PLqOG3cBizTUFQw5lwpkCpQFCVuxHTZgrT)

# There are other spectroscopic linelists but we are focusing on HITRAN only here. In addition, there are more complex line-shapes that we are not treating here, including [collisional narrowing](https://journals.aps.org/pr/abstract/10.1103/PhysRev.89.472), [line-mixing](https://aip.scitation.org/doi/abs/10.1063/1.478724), [collision-induced absorption (CIA)](https://www.sciencedirect.com/science/article/abs/pii/S0022407311003773).

# ### Line-shape (background covered in class):
# If we have a line-strength S (in cm$^{-1}$ cm$^2$/molecule) for a specific transition at $\nu_0$, we can compute the cross section as:

# $$\sigma(\nu) = S\Phi(\nu-\nu_0)\,,$$

# where $\Phi(\nu-\nu_0)$ denotes the line-shape function (in 1/cm$^{-1}$), which is normalized to 1:

# $\int_{-\infty}^\infty \Phi(\nu-\nu_0) d\nu=1$

# There are several processes that affect the shape and width of $\Phi$ and we will walk through the most important ones here now.

# ##### Doppler Broadening
# Doppler broadening is caused by a simple doppler shift of emitted and absorbed frequencies, caused by the relative velocities of the molecules along the line of sight. A doppler shifted apparent frequency from the centroid frequency $\nu_0$ can be described as:

# $$\nu = \nu_0\left(1+\frac{v_r}{c}\right)\,$$

# where $v_r$ is the relative velocity of the absorbing photon along the line of sight. The doppler shift is then simply

# $$\Delta\nu = \nu_0\frac{v_r}{c}$$

# Let us take a simple example with the satellite flying at 7km/s and either staring into the flight direction (technically, it wouldn't see the atmosphere then but let's ignore this) or looking into the back:

# ____

## Speed of light (in m/s)
c = 299792458.0

## Relative velocity (in m/s)
vᵣ = 7000.0

## Center wavenumber (say 1600nm, which is 1e7/1600=6250 cm$^{-1}$)
v₀ = 6250.0

## Doppler shift:
Δ_ν = ( v₀ * vᵣ) / c

## Just writing out doppler shift in wavenumbers and wavelengths
println("Doppler shift = $Δ_ν cm-1")
println("Doppler shift = $(1e7/(v₀-Δ_ν)-1e7/v₀) nm")

# #### Random movements of molecules in gases lead to doppler broadening effects

# In the one-dimensional case along say the x-axis (we only observe the distribution along one line of sight), the speed of molecules is distributed according to the Maxwell-Boltzmann distribution:

# $$f(v_x) = \sqrt{\frac{m}{2\pi kT}}\exp{\left(-\frac{mv_x^2}{2kT}\right)}.$$

# We can then define a Doppler width $\Delta \nu_D$ as

# $$\Delta \nu_D = \frac{\nu_0}{c}\sqrt{\frac{2kT}{m}}= \frac{\nu_0}{c}\sqrt{\frac{2RT}{M}}\,$$

# which yields the following line-shape:

# $$\phi_D(\nu) = \frac{1}{\Delta \nu_D \sqrt{\pi}}\exp{\left(-\frac{(\nu-\nu_0)^2}{\Delta \nu^2_D}\right)}=\,,$$

# which is a Gaussian distribution.

# Let us put in some numbers with R=8.3144598 J/K/mol at 6000cm$^{-1}$:

#
T = 220K, 290K

#
M = 16g/mol (CH$4$) or 44g/mol (CO$_2$)

# ---

# $\Delta \nu_D(290K,CO_2)=0.0066cm^{-1}$
# $\Delta \nu_D(220K,CO_2)=0.0058cm^{-1}$
# $\Delta \nu_D(290K,CH_4)=0.0110cm^{-1}$
# $\Delta \nu_D(290K,CH_4)=0.0096cm^{-1}$

# Multitply with about 1.6585 (2$\sqrt{ln(2)}$) to get FWHM

# ----

# ### Natural and collisional broadening

# Due to the Heisenberg principle

# $$\Delta E \Delta t \sim h/2\pi\,,$$

# As $\Delta E$ is $h\Delta\nu$, we can write:

# $$\Delta\nu = \frac{h/2\pi}{\tau}$$

# The natural line-width is defined by using the radiative lifetime but this is mostly negligible as the natural lifetime of the upper state is usually much much smaller than the "perturbed" lifetime in the presence of quencher (e.g. through collisions). Again, there are exceptions.

# Collisions between molecules reduce the lifetime of the upper state, thus widening the line width. This behavor gives rise to the so-called Lorentz lineshape.

# $$\phi_L(\nu) = \frac{\alpha_L}{\pi \left[(\nu-\nu_0)^2+\alpha_L^2\right]}$$

#  $\alpha_L$ depends linearly on the number density of the perturbing molecules and the relative speed of the collision partners (thus scales linearly with pressure and with $\sqrt{T}$, basically both density and velocity affect the number of collisions per time).

# ----

# ### Voigt lineshape
# The Voigt line-shape is the combination of Doppler and Pressure broadening (convolution of the two) but cannot be evaluated analytically. However, there are numerical routines to compute it efficiently.

# ----

# ### Other more complex lineshapes

# Once you dig deeper, there are various other more complex line-shapes (and line-mixing effects), which we ignore for now as the Voigt line-shape can provide very reasonable results. See, for instance, [here](https://www.degruyter.com/view/j/pac.2014.86.issue-12/pac-2014-0208/pac-2014-0208.xml)

using Plots
using Pkg.Artifacts
using vSmartMOM
using vSmartMOM.Absorption

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

co2_par = Absorption.read_hitran(artifact("CO2"), mol=2, iso=1, ν_min=6214.4, ν_max=6214.8);
line_voigt   = make_hitran_model(co2_par, Voigt(), architecture=CPU())
line_doppler = make_hitran_model(co2_par, Doppler(), architecture=CPU())
line_lorentz = make_hitran_model(co2_par, Lorentz(), architecture=CPU())
## Specify our wavenumber grid
ν = 6210:0.001:6220;

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

cs_co2_1atm   = absorption_cross_section(line_voigt, ν, 1013.0     , 296.0);
cs_co2_075atm = absorption_cross_section(line_voigt, ν, 0.75*1013.0, 296.0);
cs_co2_05atm  = absorption_cross_section(line_voigt, ν, 0.5*1013.0 , 296.0);
cs_co2_025atm = absorption_cross_section(line_voigt, ν, 0.25*1013.0, 296.0);
cs_co2_01atm  = absorption_cross_section(line_voigt, ν, 0.1*1013.0 , 296.0);

#Get some more line-shapes just for Doppler and Voigt
cs_co2_01atm    = absorption_cross_section(line_voigt,   ν, 0.1*1013.0 , 296.0);
cs_co2_doppler  = absorption_cross_section(line_doppler, ν, 0.1*1013.0 , 296.0);
cs_co2_lorentz  = absorption_cross_section(line_lorentz, ν, 0.1*1013.0 , 296.0);

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

#plotly()
ff = 1e20;
plot(ν,  ff*cs_co2_1atm,   label="Voigt, 1atm", yformatter = :scientific)
plot!(ν, ff*cs_co2_075atm, label="Voigt, 0.75atm")
plot!(ν, ff*cs_co2_05atm,  label="Voigt, 0.5atm")
plot!(ν, ff*cs_co2_025atm, label="Voigt, 0.25atm")

xlims!((6214,6215.2))
xlabel!("Wavenumber (cm⁻¹))")

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

plot( ν, cs_co2_01atm   /maximum(cs_co2_01atm) ,label="Voigt, 296K, 0.1atm")
plot!(ν, cs_co2_doppler /maximum(cs_co2_01atm) ,label="Doppler, 296K")
plot!(ν, cs_co2_lorentz /maximum(cs_co2_01atm) ,label="Lorentz, 296K, 0.1atm")
xlims!((6214.4,6214.8))

ylabel!("σ/max(σ)")
xlabel!("Wavenumber (cm⁻¹)")

# ### From an individual line to a band

# Here, we will just compute an entire band of CO$_2$ (a few to be precise) and look at some simple behavior, e.g. the re-distribution of individual lines in the P and R branch with changing temperature.

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

co2_par_band = Absorption.read_hitran(artifact("CO2"), mol=2, iso=1, ν_min=6000.0, ν_max=6400.0);
band_voigt   = make_hitran_model(co2_par_band , Voigt(), architecture=CPU())

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

ν_band = 6300:0.01:6400;
σ_co2_Voigt220 = absorption_cross_section(band_voigt, ν_band, 1013.0 , 220.0);
σ_co2_Voigt290 = absorption_cross_section(band_voigt, ν_band, 1013.0 , 290.0);

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

plot( ν_band, ff*σ_co2_Voigt220, alpha=0.5, label="220K")
plot!(ν_band, ff*σ_co2_Voigt290, alpha=0.5, label="290K")
xlims!((6300,6380))

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

plot(ν_band, ff*(σ_co2_Voigt220 - σ_co2_Voigt290), label="220K-290K")
xlims!((6300,6380))

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

# Think about this temperature dependence!!
#
 Why is it there?
#
 Can you see the P and R branch?
#
 Can you identify different J levels?

# ---
# ### Some fun stuff

T = 290.0
@gif for p = 10:10:1100
    σ = absorption_cross_section(band_voigt, ν_band, p , T);
    plot(ν_band, ff*σ, yaxis=:log,label="p=$p")
    ylims!((1e-7, 1e-1))
end

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

p = 900.0
@gif for T = 10:10:320
    σ = absorption_cross_section(band_voigt, ν_band, p , T);
    plot(ν_band, ff*σ, yaxis=:log, label="T=$T")
    ylims!((1e-7, 1e-1))
end

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

## More extreme case, let's take 10 atmospheres (10,000hPa)
σ = absorption_cross_section(band_voigt, ν_band, 10000.0 , 300.0);
plot(ν_band, ff*σ, label="p=10000.0hPa")

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


text:
using SamplingHS3
using Plots
using GR
# 3d gaussian mixture test..

function plot_3d_view( x_3d::Array{Float64,2} ; xf = [] )
    if(isempty(xf)); xf = ones(size(x_3d)); end
    gr()
    p1 = scatter( x_3d[1,:] , x_3d[2 ,:] , x_3d[3,:] , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none)
    p2 = scatter( x_3d[1,:] , x_3d[2 ,:]  , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none)
    p3 = scatter( x_3d[1,:] , x_3d[3 ,:]  , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none )
    p4 = scatter( x_3d[2,:] , x_3d[3 ,:]  , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none )

    display(plot(p1,p2,p3,p4,layout=(2,2)))
end

## Create Gausisan mixture
gm3_mix = []
for zi=1:10;
    found_es = false
    vd_pre   = NaN * ones(3,3)
    while(!found_es)
        vd = qr(rand(3,3))[1]
        vd_pre = vd*diagm(rand(3))*vd'
        found_es = issymmetric(vd_pre)
    end
    gpi = SamplingHS3.GaussParam( 8*rand(3) , vd_pre )
    gm3_mix = [ gm3_mix ;  gpi ]
end
gm3 = SamplingHS3.GaussMixture(gm3_mix,0.5+rand(10))
bb3_G = [eye(3);-eye(3)]
bb3_h = [10.;10;10;10;10;10]

## Create cost function (i.e. -log(f(x)) for a pdf given by f(x) )
f_cost_3d = (x) -> -log(SamplingHS3.pdf_gaussmixture_nonnormalized(gm3,x))

## Test Plot
xx3  = SamplingHS3.draw_gaussmixture( gm3 , 100_000 )
plot_3d_view(xx3[:,1:10_000],xf=1:10000)
#scatter(xx3[1,1:2000] , xx3[2 ,1:2000] , xx3[3,1:2000] , ma = 0.1 , m=2, msw=0 , c=:lightrainbow , lab="grad")

## Sample with SamplingHS3.jl
X_init = SamplingHS3.sample_unif_hr( bb3_G , bb3_h , zeros(3) , 400 , 2000 )
ns_result = SamplingHS3.nested_sampling( X_init , f_cost_3d , bb3_G , bb3_h , SamplingHS3.NSConfig(8,0.95,SamplingHS3.StopCriterion(max_time=10.0),SamplingHS3.MCMCConfig()) )

## To generate an actual sample, we have to resample the ns run:
rs_result = SamplingHS3.resample_ns_data( ns_result.ns_iter_data )

## The function sim_shell_compression simulates the shrinking process, and eval_integral
ssc_data = SamplingHS3.sim_shell_compression(ns_result.ns_iter_data)
volume_bounding_box = 20. * 20. * 20.
ei_result = SamplingHS3.eval_integral(volume_bounding_box,ns_result.ns_iter_data)

# plot the shell compression process:
p_comp = plot( log.( ssc_data[2]' ) )
#plot the integration process:
p_integ = plot( ei_result[3]'  , la=0.1 , leg=:none )
p_cum_integ = plot( cumsum( ei_result[3] , 2 )' , lc=:grey , la=0.1 , leg=:none )

## Plot the animcated sampling process (just because it's cool)
x_live_1 = SamplingHS3.reconstruct_live_samples(ns_result.ns_iter_data)
x_live  = x_live_1[1]
xf_live = x_live_1[2]
for zi in 1:length(xf_live)
    #xf_zi = 1:length(xf_live) )
    xf_zi = xf_live[zi]
    display( plot_3d_view( x_live[zi] , xf = xf_zi ) )
end

#binned_a = SamplingHS3.bin2d( ns_x_sequence[[1,2],:] , ones(size(ns_x_sequence,2)) , 40,40 , x -> sum(x) )
binned_a = SamplingHS3.bin2d( ns_x_sequence[[1,2],1:3000] , ones(3000) , 40,40 , x -> sum(x) )
heatmap(binned_a)
binned_b = SamplingHS3.bin2d( ns_x_sequence[[1,3],1:3000] , ones(3000) , 40,40 , x -> sum(x) )
heatmap(binned_b)
scatter!(rs_result[1][1,:],rs_result[1][2,:],rs_result[1][3,:])


text:
@testset "Int" begin
    @testset "affirm: 1-10" begin
    @test all(1:10) do a
        @match a begin
            1 => a == 1
            2 => a == 2
            3 => a == 3
            4 => a == 4
            5 => a == 5
            6 => a == 6
            7 => a == 7
            8 => a == 8
            9 => a == 9
            10 => a == 10
        end
    end
    end
    @testset "negate: 1-10" begin
    @test !any(1:10) do a
        @match a begin
            1 => a != 1
            2 => a != 2
            3 => a != 3
            4 => a != 4
            5 => a != 5
            6 => a != 6
            6 => a != 6
            7 => a != 7
            8 => a != 8
            9 => a != 9
            10 => a != 10
        end
    end
    end
end

@testset "Generic Number" begin
    _1_10 = map(Float64, 1:10)
    @testset "affirm float32: 1-10" begin
    @test all(_1_10) do a
        @match a begin
            1.0 => a == 1.0
            2.0 => a == 2.0
            3.0 => a == 3.0
            4.0 => a == 4.0
            5.0 => a == 5.0
            6.0 => a == 6.0
            7.0 => a == 7.0
            8.0 => a == 8.0
            9.0 => a == 9.0
            10.0 => a == 10.0
        end
    end
    end
    @testset "negate float32: 1-10" begin
    foreach(_1_10) do a
        @test !(@match a begin
            1.0 => a != 1.0
            2.0 => a != 2.0
            3.0 => a != 3.0
            4.0 => a != 4.0
            5.0 => a != 5.0
            6.0 => a != 6.0
            6.0 => a != 6.0
            7.0 => a != 7.0
            8.0 => a != 8.0
            9.0 => a != 9.0
            10.0 => a != 10.0
        end)
    end
    end
end

@testset "String" begin
    @test @match "aaaa" begin
        "1234" => false
        "aaaa" => true
        1234   => false
    end
end

@testset "OrPattern" begin
    @test @match 1 begin
        2 => false
        3 => false
        (4 || 5 || 1) => true
    end
end

@testset "AndPattern" begin
    @test @match 2 begin
        "1" || 1 => false
        2   && 1 => false
        10  || 2 && ::Int => true
    end
end

@testset "Predicate" begin
    @test @match 10 begin
        x && if x^2 > 200 end => false
        x && if (x-3)^3 < 125 end => false
        x => x === 10
    end
end

@testset "Wildcard" begin
    @test @match 100 begin
        1 && 100 => false
        _        => true
    end
end

@testset "RefPattern" begin
    @test @match 100 begin
        (&100)   =>  true
        _        =>  false
    end
    a = "123"
    @test @match a begin
        (&a)     =>  true
        _        =>  false
    end
end

@testset "Recognizer(AppPattern)" begin
    @data internal TestRecog begin
        TestRecog_A(Int, Int)
        TestRecog_B(a :: Float64, b :: String)
    end

    a = TestRecog_A(1, 2)
    b = TestRecog_B(1.0, "2")

    @testset "Wildcard" begin
        @test @match a begin
            TestRecog_A(_) => true
            _ => false
        end
    end
    @testset "Ordered Anonymous Fields" begin
        @test @match a begin
            TestRecog_A(1, 2) => true
            _ => false
        end

        @test @match b begin
            TestRecog_B(1.0, "2") => true
            _ => false
        end
    end

    @testset "Named Fields" begin
        @test @match a begin
            TestRecog_A(_2 = 2) => true
            _ => false
        end

        @test @match b begin
            TestRecog_B(a = 1.0) => true
            _ => false
        end
    end
end

@testset "null-value dict pattern" begin
    x = Dict(:a => nothing)
    @test @match x begin
        Dict(:a => a) => a === nothing
        _ => false
    end

end

struct TestGH end

@testset "Generalized Recognizer(GAppPattern)" begin
    @use GADT
    @data internal TestGRecog{T} begin
        TestGRecog_A{T, A} :: (A, T) => TestGRecog{T}
        TestGRecog_B{T, B} :: (a :: T, b :: B) => TestGRecog{T}
    end

    a = TestGRecog_A(1, TestGH())
    b = TestGRecog_B([1], "2")

    @testset "Matching for specialization" begin
        @test @match a begin
            TestGRecog_A(_) => true
            _ => false
        end

        @test @match a begin
            TestGRecog_A{TestGH, Int}(_) => true
            _ => false
        end

        @test @match b begin
            TestGRecog_B{Vector{Int}, String}(_) => true
            _ => false
        end

        @test @match b begin
            ::TestGRecog{Vector{Int}} => true
            _ => false
        end

        @test @match a begin
            ::TestGRecog{TestGH} => true
            _ => false
        end

        @test @match a begin
            TestGRecog_A{TestGH}(_) => true
            _ => false
        end

    end

    @testset "Matching for generalization" begin

        @test @match a begin
            TestGRecog_A{T, A}(::A, ::T) where {A, T} => true
            _ => false
        end

        @test @match a begin
            TestGRecog_A{TestGH, A}(_) where A <: Number => true
            _ => false
        end

        @test @match a begin
            TestGRecog_A{TestGH}(_) => true
            _ => false
        end

        @test @match b begin
            TestGRecog_B{A, B}(::B, ::A) where {A, B} => false
            TestGRecog_B{A, B}(::A, ::B) where {A, B} => true
            _ => false
        end

    end

end