text:
module PeriodicDarcyTests

using Gridap
using Test
using LinearAlgebra

u(x) = VectorValue(x[1]*(x[1]-1)*(2x[2]-1.0),-x[2]*(x[2]-1.0)*(2x[1]-1.0))
p(x) = x[2]-0.5
f(x) = VectorValue(x[1]*(x[1]-1)*(2x[2]-1.0),-x[2]*(x[2]-1.0)*(2x[1]-1.0)+1.0)
g(x) = 0.0

domain = (0,1,0,1,0)
order = 1
partition = (4,4)
model = CartesianDiscreteModel(domain,partition; isperiodic=(true,false))

V = FESpace(
reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64},
conformity=:Hdiv, model=model,dirichlet_tags=collect(1:6))

Q = FESpace(
reffe=:QLagrangian, order=order, valuetype=Float64,
conformity=:L2, model=model, constraint=:zeromean)

U = TrialFESpace(V,u)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

trian = Triangulation(model)
degree = 2*(order)
quad = CellQuadrature(trian,degree)
x = get_physical_coordinate(trian)

function a(x,y)
  u, p = x
  v, q = y
  v⋅u - p*(∇⋅v) + q*(∇⋅u)
end

function l(y)
  v, q = y
  v⋅f + q*g
end

t_Ω = AffineFETerm(a,l,trian,quad)
op = AffineFEOperator(X,Y,t_Ω)
xh = solve(op)
uh, ph = xh

eu = u - uh
ep = p - ph

l2(v) = v⋅v
eu_l2 = sum(integrate(l2(eu),trian,quad))
ep_l2 = sum(integrate(l2(ep),trian,quad))

tol = 1.0e-9
@test eu_l2 < tol
@test ep_l2 < tol

end # module


text:
export BatchNorm
import Knet
using Knet.Ops21
using Knet.KnetArrays: DevArray
using AutoGrad: Param, value
using CUDA.CUDNN:
    cudnnNormalizationForward,
    cudnnNormalizationForward!,
    cudnnNormalizationForwardInference,
    cudnnNormalizationForwardTraining,
    cudnnNormalizationBackward,
    cudnnGetNormalizationForwardTrainingWorkspaceSize,
    cudnnGetNormalizationTrainingReserveSpaceSize,
    cudnnActivationDescriptor,
    cudnnTensorDescriptor,
    cudnnNormMode_t,
        CUDNN_NORM_PER_ACTIVATION, # 0, bnScale, bnBias tensor dims are 1xCxHxWx.. (one value per CHW...-slice, normalized over N slice)
        CUDNN_NORM_PER_CHANNEL,    # 1, bnScale, bnBias tensor dims are 1xCx1x1 (one value per C-dim normalized over Nx1xHxW subtensors)
    cudnnNormOps_t,
        CUDNN_NORM_OPS_NORM,                # 0, /* do normalization only */
        CUDNN_NORM_OPS_NORM_ACTIVATION,     # 1, /* do Norm, then activation */
        CUDNN_NORM_OPS_NORM_ADD_ACTIVATION, # 2, /* do Norm, then elemWiseAdd, then activation */
    cudnnNormAlgo_t,
        CUDNN_NORM_ALGO_STANDARD, # 0
        CUDNN_NORM_ALGO_PERSIST,  # 1
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,      # 0
        CUDNN_ACTIVATION_RELU,         # 1
        CUDNN_ACTIVATION_TANH,         # 2
        CUDNN_ACTIVATION_CLIPPED_RELU, # 3
        CUDNN_ACTIVATION_ELU,          # 4
        CUDNN_ACTIVATION_IDENTITY,     # 5
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C, # 2, /* each image point is vector of element of C, vector length in data type */
    handle

"""
    b = BatchNorm(; kwargs...)
    b(x)
"""
mutable struct BatchNorm
    dims
    mean
    var
    bias
    scale
    use_estimates
    update::Float64
    epsilon::Float64
    format::cudnnTensorFormat_t
    mode::cudnnNormMode_t
    savedMean
    savedVar
    workspace
    reserveSpace
    dx::Ref{Any}
    dscale::Ref{Any}
    dbias::Ref{Any}
end

function BatchNorm(
    ;
    dims = nothing,
    mean = nothing,
    var = nothing,
    bias = nothing,
    scale = nothing,
    use_estimates = nothing,
    update::Real = Cdouble(0.1),
    epsilon::Real = Cdouble(1e-5),
)
    format = CUDNN_TENSOR_NCHW
    mode = CUDNN_NORM_PER_CHANNEL
    savedMean = nothing
    savedVar = nothing
    workspace = nothing
    reserveSpace = nothing
    dx = Ref{Any}(nothing)
    dscale = Ref{Any}(nothing)
    dbias = Ref{Any}(nothing)
    BatchNorm(dims, mean, var, bias, scale, use_estimates, update, epsilon, format, mode, savedMean, savedVar, workspace, reserveSpace, dx, dscale, dbias)
end

# Some fields can only be initialized after seeing the first input x

function initBatchNorm(b::BatchNorm, x; use_estimates)
    n = ndims(x)
    if b.dims === nothing || b.dims === (1:(n-2)...,n)
        b.mode, b.format, bsize = CUDNN_NORM_PER_CHANNEL, CUDNN_TENSOR_NCHW, ntuple(i->(i===n-1 ? size(x,i) : 1), n)
    elseif b.dims === (2:n...,)
        b.mode, b.format, bsize = CUDNN_NORM_PER_CHANNEL, CUDNN_TENSOR_NHWC, ntuple(i->(i===1 ? size(x,i) : 1), n)
    elseif length(b.dims) === 1 && b.dims[1] === n
        b.mode, b.format, bsize = CUDNN_NORM_PER_ACTIVATION, CUDNN_TENSOR_NCHW, ntuple(i->(i===n ? 1 : size(x,i)), n)
    else
        error("x=$(size(x)) dims=$dims not supported")
    end
    issimilar(u,v,s=size(v))=(typeof(value(u)) === typeof(value(v)) && size(u) === s)
    b.mean === nothing ? b.mean = fill!(similar(x, bsize), 0) : @assert issimilar(b.mean, x, bsize) "x=$(summary(x)) mean=$(summary(b.mean)) bsize=$bsize"
    b.var === nothing ? b.var = fill!(similar(x, bsize), 1) : @assert issimilar(b.var, x, bsize)
    b.bias === nothing ? b.bias = Param(fill!(similar(x, bsize), 0)) : @assert issimilar(b.bias, x, bsize)
    b.scale === nothing ? b.scale = Param(fill!(similar(x, bsize), 1)) : @assert issimilar(b.scale, x, bsize)
    if !use_estimates && x isa DevArray # cudnnNormalizationForwardTraining will be called
        b.savedMean === nothing ? b.savedMean = similar(x, bsize) : @assert issimilar(b.savedMean, x, bsize)
        b.savedVar === nothing ? b.savedVar = similar(x, bsize) : @assert issimilar(b.savedVar, x, bsize)
        workspaceSize, reserveSpaceSize = cudnnNormalizationTempSpaceSizes(b, x)
        if sizeof(b.reserveSpace) < reserveSpaceSize; b.reserveSpace = cudnnTempSpace(reserveSpaceSize); end
        if sizeof(b.workspace) < workspaceSize; b.workspace = cudnnTempSpace(workspaceSize); end
        b.dx[] === nothing ? b.dx[] = similar(x) : @assert issimilar(x, b.dx[])
        b.dscale[] === nothing ? b.dscale[] = similar(b.scale) : @assert issimilar(b.scale, b.dscale[])
        b.dbias[] === nothing ? b.dbias[] = similar(b.bias) : @assert issimilar(b.bias, b.dbias[])
    end
end

function cudnnNormalizationTempSpaceSizes(b::BatchNorm, x)
    workspaceSize, reserveSpaceSize = Ref{Csize_t}(0), Ref{Csize_t}(0)
    xDesc, bDesc = (u->cudnnTensorDescriptor(value(u); format=b.format)).((x, b.mean))
    mode, normOps, algo, groupCnt = b.mode, CUDNN_NORM_OPS_NORM, CUDNN_NORM_ALGO_STANDARD, 1
    cudnnGetNormalizationForwardTrainingWorkspaceSize(handle(), mode, normOps, algo, xDesc, C_NULL, xDesc, bDesc, C_NULL, bDesc, workspaceSize, groupCnt)
    cudnnGetNormalizationTrainingReserveSpaceSize(handle(), mode, normOps, algo, C_NULL, xDesc, reserveSpaceSize, groupCnt)
    workspaceSize[], reserveSpaceSize[]
end

function (b::BatchNorm)(x)
    use_estimates = (b.use_estimates === nothing ? !Knet.training() : b.use_estimates)
    update = Knet.training() ? b.update : 0
    initBatchNorm(b, x; use_estimates)
    batchnorm(x, b.mean, b.var, b.bias, b.scale; use_estimates, update,
              b.epsilon, b.mode, b.format, b.savedMean, b.savedVar, b.workspace, b.reserveSpace, b.dx, b.dscale, b.dbias)
end

ext: jl


text:
using OREnvironment
using ORInterface
using Test

#######################
# PRELIMINARIES FOR THE TESTS
#######################
mutable struct MyProblem <: Problem
    costs::Array{Float64,1};
    constraints::Array{<:Constraint,1};
    variablesConstraints::Array{Array{Int,1},1};
end

function constructorProblem()
    cost = collect(1.0:6.0);

    variables1 = [1, 3, 4, 6];
    variables2 = [1, 3, 5, 6];
    coefficients1 = [2.3, 3.2, 3.1, 12.34];
    coefficients2 = coefficients1 .+ 1.0;
    constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefficients1);
    constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefficients2);
    constraints = [constraint1, constraint2];

    variablesConstraints = [[1,2], Int[], [1,2], [1], [2], [1,2]];

    return MyProblem(cost, constraints, variablesConstraints);
end

function constructorSolution(numConstraints=2)
    status = OREnvironment.constructStatus(numConstraints);
    Tvariables = Float64; sizeArray = 6;
    sol  = OREnvironment.constructSolution(:FixedLengthArray, (Tvariables, sizeArray, status));
    return sol;
end

function constructorProblemNegative()
    cost = collect(1.0:6.0);

    variables1 = [1, 3, 4, 6];
    variables2 = [1, 3, 5, 6];
    coefficients1 = -[2.3, 3.2, 3.1, 12.34];
    coefficients2 = coefficients1 .- 1.0;
    constraint1 = OREnvironment.constructConstraint(15.0, :greaterOrEq, variables1, coefficients1);
    constraint2 = OREnvironment.constructConstraint(9.0, :greaterOrEq, variables2, coefficients2);
    constraints = [constraint1, constraint2];

    variablesConstraints = [[1,2], Int[], [1,2], [1], [2], [1,2]];

    return MyProblem(cost, constraints, variablesConstraints);
end

#######################
# TESTS
#######################
@testset "Building Solution" begin
    @testset "FixedLengthArray" begin
        numConstraints = 5;
        sol = constructorSolution(numConstraints);

        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.is_optimal(sol) == false;
        @test OREnvironment.get_objfunction(sol) == 0.0;
        for i in 1:numConstraints
            @test OREnvironment.get_constraint_consumption(sol, i) == 0.0;
        end

        OREnvironment.set_feasible!(sol, true);
        OREnvironment.set_optimal!(sol, true);
        OREnvironment.set_objfunction!(sol, 16.0);
        for i in 1:numConstraints
            OREnvironment.set_constraint_consumption!(sol, i, 12.0);
        end
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.is_optimal(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 16.0;
        for i in 1:numConstraints
            @test OREnvironment.get_constraint_consumption(sol, i) == 12.0;
        end

        st = OREnvironment.constructStatus(numConstraints);
        OREnvironment.set_feasible!(st, false);
        OREnvironment.set_objfunction!(st, 125.63);
        OREnvironment.set_constraint_consumption!(st, 1, 3.0);
        OREnvironment.set_constraint_consumption!(st, 5, 4.0);
        OREnvironment.update_status!(sol.status, st);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.is_optimal(sol) == false;
        @test OREnvironment.get_objfunction(sol) == 125.63;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 3.0;
        @test OREnvironment.get_constraint_consumption(sol, 5) == 4.0;
    end
end

@testset "create_empty_solution" begin
    p = constructorProblem()
    sol = OREnvironment.create_empty_solution(p, Int)
    @test OREnvironment.is_feasible(sol) == false;
    @test OREnvironment.is_optimal(sol) == false;
    @test OREnvironment.get_objfunction(sol) == 0.0;
    @test length(sol.sol) == 6;
end

@testset "Adding and removing solutions" begin
    value = 12.3;
    valueDummy = 0.0;
    @testset "FixedLengthArray" begin
        sol  = constructorSolution();
        OREnvironment.add_solution!(sol, 3, value);
        OREnvironment.add_solution!(sol, 1, 2value);
        OREnvironment.add_solution!(sol, 2, 3value);
        @test OREnvironment.get_solution(sol, 3) == value;
        @test OREnvironment.get_solution(sol, 1) == 2value;
        @test OREnvironment.get_solution(sol, 2) == 3value;
        OREnvironment.remove_solution!(sol, 3);
        @test OREnvironment.get_solution(sol, 3) == 0.0;
        @test OREnvironment.get_solution(sol, 1) == 2value;
        @test OREnvironment.get_solution(sol, 2) == 3value;
        OREnvironment.remove_all_solutions!(sol);
        @test length(sol.sol) == 6;
        for v in sol.sol
            @test v == 0.0;
        end
    end
end

@testset "is first solution better" begin
    numConstraints = 5;
    status1 = OREnvironment.constructStatus(numConstraints);
    status2 = OREnvironment.constructStatus(numConstraints);
    Tvariables = Int; sizeArray = 6;
    s1 = OREnvironment.constructSolution(:FixedLengthArray, (Tvariables, sizeArray, status1));
    s2 = OREnvironment.constructSolution(:FixedLengthArray, (Tvariables, sizeArray, status2));
    OREnvironment.set_objfunction!(s1, 34.0);
    OREnvironment.set_objfunction!(s2, 3.0);

    @testset "feasibility not required" begin
        feasibility = false;
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, false);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, false);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;
    end

    @testset "feasibility required" begin
        feasibility = true;
        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0)
        OREnvironment.set_feasible!(s1, false);
        OREnvironment.set_feasible!(s2, false);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, false);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == true;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, false);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == true;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_solution_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_solution_better(s2,s1,:min, feasibility) == false;
    end
end

@testset "copy_first_solution_to_second" begin
    s1 = constructorSolution();
    s2 = constructorSolution();
    p = constructorProblem();

    OREnvironment.add_solution!(s1, 1, 7.0);
    OREnvironment.add_solution!(s1, 2, 7.4);
    OREnvironment.add_solution!(s1, 3, 3.3);
    OREnvironment.set_constraint_consumption!(s1, 1, 26.66);
    OREnvironment.set_constraint_consumption!(s1, 2, 36.96);
    OREnvironment.set_objfunction!(s1, 31.7);

    OREnvironment.add_solution!(s2, 4, 73.4);
    OREnvironment.add_solution!(s2, 5, 8.0);
    OREnvironment.set_constraint_consumption!(s2, 1, 227.54);
    OREnvironment.set_constraint_consumption!(s2, 2, 32.8);
    OREnvironment.set_objfunction!(s2, 333.6);

    @test OREnvironment.get_objfunction(s1) == 31.7;
    @test OREnvironment.get_objfunction(s2) == 333.6;
    @test OREnvironment.is_feasible(s1) == false;
    @test OREnvironment.is_feasible(s2) == false;
    @test OREnvironment.is_optimal(s1) == false;
    @test OREnvironment.is_optimal(s2) == false;
    @test OREnvironment.get_constraint_consumption(s1, 1) ≈ 26.66;
    @test OREnvironment.get_constraint_consumption(s1, 2) ≈ 36.96;
    @test OREnvironment.get_constraint_consumption(s2, 1) ≈ 227.54;
    @test OREnvironment.get_constraint_consumption(s2, 2) == 32.8;

    @test OREnvironment.get_solution(s1, 1) == 7.0;
    @test OREnvironment.get_solution(s1, 2) == 7.4;
    @test OREnvironment.get_solution(s1, 3) == 3.3;
    @test OREnvironment.get_solution(s1, 4) == 0.0;
    @test OREnvironment.get_solution(s1, 5) == 0.0;
    @test OREnvironment.get_solution(s1, 6) == 0.0;

    @test OREnvironment.get_solution(s2, 1) == 0.0;
    @test OREnvironment.get_solution(s2, 2) == 0.0;
    @test OREnvironment.get_solution(s2, 3) == 0.0;
    @test OREnvironment.get_solution(s2, 4) == 73.4;
    @test OREnvironment.get_solution(s2, 5) == 8.0;
    @test OREnvironment.get_solution(s2, 6) == 0.0;

    OREnvironment.set_feasible!(s1, true);
    OREnvironment.set_optimal!(s1, true);
    @test OREnvironment.is_feasible(s1) == true;
    @test OREnvironment.is_feasible(s2) == false;
    @test OREnvironment.is_optimal(s1) == true;
    @test OREnvironment.is_optimal(s2) == false;
    OREnvironment.copy_first_solution_to_second!(s1,s2);
    @test OREnvironment.get_solution(s1, 1) == 7.0;
    @test OREnvironment.get_solution(s1, 2) == 7.4;
    @test OREnvironment.get_solution(s1, 3) == 3.3;
    @test OREnvironment.get_solution(s1, 4) == 0.0;
    @test OREnvironment.get_solution(s1, 5) == 0.0;
    @test OREnvironment.get_solution(s1, 6) == 0.0;
    @test OREnvironment.get_solution(s2, 1) == 7.0;
    @test OREnvironment.get_solution(s2, 2) == 7.4;
    @test OREnvironment.get_solution(s2, 3) == 3.3;
    @test OREnvironment.get_solution(s2, 4) == 0.0;
    @test OREnvironment.get_solution(s2, 5) == 0.0;
    @test OREnvironment.get_solution(s2, 6) == 0.0;

    @test OREnvironment.get_objfunction(s1) == 31.7;
    @test OREnvironment.get_objfunction(s2) == 31.7;
    @test OREnvironment.is_feasible(s1) == true;
    @test OREnvironment.is_feasible(s2) == true;
    @test OREnvironment.is_optimal(s1) == true;
    @test OREnvironment.is_optimal(s2) == true;
    @test OREnvironment.get_constraint_consumption(s1, 1) ≈ 26.66;
    @test OREnvironment.get_constraint_consumption(s1, 2) ≈ 36.96;
    @test OREnvironment.get_constraint_consumption(s2, 1) ≈ 26.66;
    @test OREnvironment.get_constraint_consumption(s2, 2) ≈ 36.96;
end

#######################
# TESTS WHEN DEALING WITH CONSTRAINTS
#######################

@testset "update_constraint_consumption" begin
    sol  = constructorSolution();
    p = constructorProblem();

    @testset "feasible solution" begin
        sol.sol[1] = 1.0;
        sol.sol[2] = 10.0; # check that this doesn't appear in the constraints
        sol.sol[3] = 1.0;
        OREnvironment.update_constraint_consumption!(sol, p.constraints);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 5.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
    end

    @testset "infeasible solution" begin
        sol.sol[1] = 2.0;
        OREnvironment.update_constraint_consumption!(sol, p.constraints);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 7.8;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 10.8;
    end

    @testset "empty array of constraints" begin
        # the return is true because it would be the case of problem with no constraints.
        voidArray = Array{OREnvironment.DefaultConstraint, 1}();
        OREnvironment.update_constraint_consumption!(sol, voidArray);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 7.8;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 10.8;
    end
end

@testset "update_constraint_consumption_and_feasibility" begin
    sol  = constructorSolution();
    p = constructorProblem();

    @testset "feasible solution" begin
        sol.sol[1] = 1.0;
        sol.sol[2] = 10.0; # check that this doesn't appear in the constraints
        sol.sol[3] = 1.0;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 5.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
    end

    @testset "infeasible solution" begin
        sol.sol[1] = 2.0;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 7.8;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 10.8;
    end

    @testset "empty array of constraints" begin
        # the return is true because it would be the case of problem with no constraints.
        voidArray = Array{OREnvironment.DefaultConstraint, 1}();
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, voidArray);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 7.8;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 10.8;
    end
end

@testset "should_global_feasibility_be_checked" begin
  isSolutionFeasible  = true;
  isIncrementFeasible = true;
  @test OREnvironment.should_global_feasibility_be_checked(isSolutionFeasible, isIncrementFeasible) == false;

  isSolutionFeasible  = true;
  isIncrementFeasible = false;
  @test OREnvironment.should_global_feasibility_be_checked(isSolutionFeasible, isIncrementFeasible) == false;

  isSolutionFeasible  = false;
  isIncrementFeasible = false;
  @test OREnvironment.should_global_feasibility_be_checked(isSolutionFeasible, isIncrementFeasible) == false;

  isSolutionFeasible  = false;
  isIncrementFeasible = true;
  @test OREnvironment.should_global_feasibility_be_checked(isSolutionFeasible, isIncrementFeasible) == true;
end

@testset "update_constraint_consumption for increments" begin
    sol = constructorSolution();
    p = constructorProblem();

    sol.sol[1] = 1.0;
    sol.sol[3] = 1.0;
    OREnvironment.update_constraint_consumption!(sol, p.constraints);
    @test OREnvironment.get_constraint_consumption(sol, 1) == 5.5;
    @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;

    @testset "feasible solution" begin
        Δ = 0.2;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.32;
        sol.sol[1] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 1, Δ, p.variablesConstraints[1]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.58;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.98;
    end

    @testset "feasible solution negative sign" begin
        Δ = -0.2;
        sol.sol[1] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 1, Δ, p.variablesConstraints[1]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.32;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 5.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
    end

    @testset "infeasible solution" begin
        Δ = 1.5;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.15;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.15;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 13.65;
        sol.sol[1] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 1, Δ, p.variablesConstraints[1]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 13.6;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;
    end

    @testset "empty Array" begin
        Δ = 1.5;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 2, Δ, p.variablesConstraints[2]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 13.6;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;
    end

    @testset "empty constraints" begin
        Δ = 1.5;
        voidArray = Array{OREnvironment.DefaultConstraint, 1}();
        OREnvironment.update_constraint_consumption!(sol, voidArray, 2, Δ, p.variablesConstraints[2]);
    end

    @test OREnvironment.get_constraint_consumption(sol, 1) == 13.6;
    @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;

    @testset "from infeasible to infeasible" begin
        Δ = -1.0;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;
    end

    @testset "from infeasible to feasible" begin
        Δ = -3.0;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 6.3;
    end

    @testset "from feasible to feasible" begin
        Δ = -1.0;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.get_constraint_consumption(sol, 1) == 7.4;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 6.3;
    end

    @testset "from feasible to infeasible" begin
        Δ = 3.0;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 16.7;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 6.3;
    end
end

@testset "update_constraint_consumption_and_feasibility for increments" begin
    sol = constructorSolution();
    p = constructorProblem();

    sol.sol[1] = 1.0;
    sol.sol[3] = 1.0;
    OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints);
    @test OREnvironment.get_constraint_consumption(sol, 1) == 5.5;
    @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;

    @testset "feasible solution" begin
        Δ = 0.2;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.32;
        sol.sol[1] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 1, Δ, p.variablesConstraints[1]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.58;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.98;
    end

    @testset "feasible solution negative sign" begin
        Δ = -0.2;
        sol.sol[1] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 1, Δ, p.variablesConstraints[1]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.32;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.12;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 5.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
    end

    @testset "infeasible solution" begin
        Δ = 1.5;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.15;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 7.5;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.15;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 13.65;
        sol.sol[1] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 1, Δ, p.variablesConstraints[1]);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 13.6;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;
    end

    @testset "empty Array" begin
        Δ = 1.5;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 2, Δ, p.variablesConstraints[2]);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 13.6;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;
    end

    @testset "empty constraints" begin
        Δ = 1.5;
        voidArray = Array{OREnvironment.DefaultConstraint, 1}();
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, voidArray, 2, Δ, p.variablesConstraints[2]);
        @test OREnvironment.is_feasible(sol) == true;
    end

    OREnvironment.set_feasible!(sol, false);
    @test OREnvironment.get_constraint_consumption(sol, 1) == 13.6;
    @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;

    @testset "from infeasible to infeasible" begin
        Δ = -1.0;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 18.6;
    end

    @testset "from infeasible to feasible" begin
        Δ = -3.0;
        sol.sol[5] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 5, Δ, p.variablesConstraints[5]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 10.5;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 6.3;
    end

    @testset "from feasible to feasible" begin
        Δ = -1.0;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 7.4;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 6.3;
    end

    @testset "from feasible to infeasible" begin
        Δ = 3.0;
        sol.sol[4] += Δ;
        OREnvironment.update_constraint_consumption_and_feasibility!(sol, p.constraints, 4, Δ, p.variablesConstraints[4]);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 16.7;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 6.3;
    end
end

#######################
# TESTS WHEN DEALING WITH PROBLEMS
#######################
@testset "Adding and removing solution with constraint update" begin
    sol = constructorSolution();
    p = constructorProblem();

    @test OREnvironment.is_feasible(sol) == false;
    @test OREnvironment.get_objfunction(sol) == 0.0;
    @test OREnvironment.get_constraint_consumption(sol, 1) == 0.0;
    @test OREnvironment.get_constraint_consumption(sol, 2) == 0.0;

    @testset "Adding" begin
        # from 0 to 2.5 => increase 2.5
        OREnvironment.add_solution_and_update_status!(sol, 1, 2.5, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 2.5;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 5.75;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 8.25;
        @test OREnvironment.get_solution(sol, 1) == 2.5;
        for i in 2:6
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end
        # from 2.5 to 3 => increase 0.5
        OREnvironment.add_solution_and_update_status!(sol, 1, 3.0, p);
        @test OREnvironment.is_feasible(sol) == false;
        @test OREnvironment.get_objfunction(sol) == 3.0;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 6.9;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 9.9;
        @test OREnvironment.get_solution(sol, 1) == 3.0;
        for i in 2:6
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end

        # from 3 to 1.5 => reduce of 1.5
        OREnvironment.add_solution_and_update_status!(sol, 1, 1.5, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 1.5;
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 3.45;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 4.95;
        @test OREnvironment.get_solution(sol, 1) == 1.5;
        for i in 2:6
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end

        # from 0 to 1.75 in a variable that does not affect constraints
        OREnvironment.add_solution_and_update_status!(sol, 2, 1.75, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 5.0;
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 3.45;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 4.95;
        @test OREnvironment.get_solution(sol, 1) == 1.5;
        @test OREnvironment.get_solution(sol, 2) == 1.75;
        for i in 3:6
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end

        # from 0 to 2.25 in variable 4 (affects only to the 1st constraint)
        OREnvironment.add_solution_and_update_status!(sol, 4, 2.25, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 14.0;
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 10.425;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 4.95;
        @test OREnvironment.get_solution(sol, 1) == 1.5;
        @test OREnvironment.get_solution(sol, 2) == 1.75;
        @test OREnvironment.get_solution(sol, 4) == 2.25;
        for i in [3,5,6]
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end

        # from 2.25 to 1 in variable 4
        OREnvironment.add_solution_and_update_status!(sol, 4, 1.0, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 9.0;
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 6.55;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 4.95;
        @test OREnvironment.get_solution(sol, 1) == 1.5;
        @test OREnvironment.get_solution(sol, 2) == 1.75;
        @test OREnvironment.get_solution(sol, 4) == 1.0;
        for i in [3,5,6]
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end
    end
    @testset "Removing" begin
        OREnvironment.remove_solution_and_update_status!(sol, 4, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 5.0;
        @test OREnvironment.get_constraint_consumption(sol, 1) ≈ 3.45;
        @test OREnvironment.get_constraint_consumption(sol, 2) ≈ 4.95;
        @test OREnvironment.get_solution(sol, 1) == 1.5;
        @test OREnvironment.get_solution(sol, 2) == 1.75;
        for i in 3:6
            @test OREnvironment.get_solution(sol, i) == 0.0;
        end
    end
    @testset "Removing all" begin
        OREnvironment.set_feasible!(sol, false);
        OREnvironment.remove_all_solutions_and_update_status!(sol, p);
        @test OREnvironment.is_feasible(sol) == true;
        @test OREnvironment.get_objfunction(sol) == 0.0;
        @test OREnvironment.get_constraint_consumption(sol, 1) == 0.0;
        @test OREnvironment.get_constraint_consumption(sol, 2) == 0.0;
    end
end

@testset "Constraints with all negative consumption" begin
    s = constructorSolution();
    p = constructorProblemNegative();

    OREnvironment.add_solution_and_update_status!(s, 1, -2.0, p);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 4.6;
    @test lhs == OREnvironment.get_constraint_consumption(s, 1);
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 6.6;
    @test lhs == OREnvironment.get_constraint_consumption(s, 2);
    @test OREnvironment.is_feasible(s) == false;
    @test OREnvironment.is_feasible(s, p.constraints) == false;

    OREnvironment.add_solution_and_update_status!(s, 2, -4.5, p);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 4.6;
    @test lhs == OREnvironment.get_constraint_consumption(s, 1);
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 6.6;
    @test lhs == OREnvironment.get_constraint_consumption(s, 2);
    @test OREnvironment.is_feasible(s) == false;
    @test OREnvironment.is_feasible(s, p.constraints) == false;

    OREnvironment.add_solution_and_update_status!(s, 3, -1.5, p);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 9.4;
    @test lhs == OREnvironment.get_constraint_consumption(s, 1);
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 12.9;
    @test lhs == OREnvironment.get_constraint_consumption(s, 2);
    @test OREnvironment.is_feasible(s) == false;
    @test OREnvironment.is_feasible(s, p.constraints) == false;

    OREnvironment.add_solution_and_update_status!(s, 3, -3.5, p);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 15.8;
    @test lhs == OREnvironment.get_constraint_consumption(s, 1);
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 21.3;
    @test lhs == OREnvironment.get_constraint_consumption(s, 2);
    @test OREnvironment.is_feasible(s) == true;
    @test OREnvironment.is_feasible(s, p.constraints) == true;

    OREnvironment.remove_all_solutions_and_update_status!(s, p);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 0.0;
    @test lhs == OREnvironment.get_constraint_consumption(s, 1);
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 0.0;
    @test lhs == OREnvironment.get_constraint_consumption(s, 2);
    @test OREnvironment.is_feasible(s) == false;
    @test OREnvironment.is_feasible(s, p.constraints) == false;
end


text:
using FactCheck, Base.Test, Images, Color, FixedPointNumbers

facts("Algorithms") do
    # Comparison of each element in arrays with a scalar
    approx_equal(ar, v) = all(abs(ar.-v) .< sqrt(eps(v)))
    approx_equal(ar::Images.AbstractImage, v) = approx_equal(Images.data(ar), v)

    context("Arithmetic") do
        img = convert(Images.Image, zeros(3,3))
        img2 = (img .+ 3)/2
        @fact all(img2 .== 1.5) --> true
        img3 = 2img2
        @fact all(img3 .== 3) --> true
        img3 = copy(img2)
        img3[img2 .< 4] = -1
        @fact all(img3 .== -1) --> true
        img = convert(Images.Image, rand(3,4))
        A = rand(3,4)
        img2 = img .* A
        @fact all(Images.data(img2) == Images.data(img).*A) --> true
        img2 = convert(Images.Image, A)
        img2 = img2 .- 0.5
        img3 = 2img .* img2
        img2 = img ./ A
        img2 = (2img).^2
        # Same operations with ColorValue images
        img = Images.colorim(zeros(Float32,3,4,5))
        img2 = (img .+ RGB{Float32}(1,1,1))/2
        @fact all(img2 .== RGB{Float32}(1,1,1)/2) --> true
        img3 = 2img2
        @fact all(img3 .== RGB{Float32}(1,1,1)) --> true
        A = fill(2, 4, 5)
        @fact all(A.*img2 .== fill(RGB{Float32}(1,1,1), 4, 5)) --> true
        img2 = img2 .- RGB{Float32}(1,1,1)/2
        A = rand(Uint8,3,4)
        img = reinterpret(Images.Gray{Ufixed8}, Images.grayim(A))
        imgm = mean(img)
        imgn = img/imgm
        @fact reinterpret(Float32, Images.data(imgn)) --> roughly(convert(Array{Float32}, A/mean(A)))
    end

    context("Reductions") do
        A = rand(5,5,3)
        img = Images.colorim(A, "RGB")
        s12 = sum(img, (1,2))
        @fact Images.colorspace(s12) --> "RGB"
        s3 = sum(img, (3,))
        @fact Images.colorspace(s3) --> "Unknown"
        A = [NaN, 1, 2, 3]
        @fact Images.meanfinite(A, 1) --> roughly([2])
        A = [NaN 1 2 3;
             NaN 6 5 4]
        @test_approx_eq Images.meanfinite(A, 1) [NaN 3.5 3.5 3.5]
        @test_approx_eq Images.meanfinite(A, 2) [2, 5]'
        @test_approx_eq Images.meanfinite(A, (1,2)) [3.5]
        @fact Images.minfinite(A) --> 1
        @fact Images.maxfinite(A) --> 6
        @fact Images.maxabsfinite(A) --> 6
        A = rand(Float32,3,5,5)
        img = Images.colorim(A, "RGB")
        dc = Images.data(Images.meanfinite(img, 1))-reinterpret(RGB{Float32}, mean(A, 2), (1,5))
        @fact maximum(map(abs, dc)) --> less_than(1e-6)

        a = convert(Array{Uint8}, [1, 1, 1])
        b = convert(Array{Uint8}, [134, 252, 4])
        @fact Images.sad(a, b) --> 387
        @fact Images.ssd(a, b) --> 80699
        af = reinterpret(Ufixed8, a)
        bf = reinterpret(Ufixed8, b)
        @fact Images.sad(af, bf) --> roughly(387f0/255)
        @fact Images.ssd(af, bf) --> roughly(80699f0/255^2)
        ac = reinterpret(RGB{Ufixed8}, a)
        bc = reinterpret(RGB{Ufixed8}, b)
        @fact Images.sad(ac, bc) --> roughly(387f0/255)
        @fact Images.ssd(ac, bc) --> roughly(80699f0/255^2)
        ag = reinterpret(RGB{Ufixed8}, a)
        bg = reinterpret(RGB{Ufixed8}, b)
        @fact Images.sad(ag, bg) --> roughly(387f0/255)
        @fact Images.ssd(ag, bg) --> roughly(80699f0/255^2)
    end

    context("fft and ifft") do
        A = rand(Float32, 3, 5, 6)
        img = Images.colorim(A)
        imgfft = fft(img)
        @fact Images.data(imgfft) --> roughly(fft(A, 2:3))
        @fact Images.colordim(imgfft) --> 1
        img2 = ifft(imgfft)
        @fact img2 --> roughly(reinterpret(Float32, img))
    end

    context("Array padding") do
        A = [1 2; 3 4]
        @fact Images.padindexes(A, 1, 0, 0, "replicate") --> [1,2]
        @fact Images.padindexes(A, 1, 1, 0, "replicate") --> [1,1,2]
        @fact Images.padindexes(A, 1, 0, 1, "replicate") --> [1,2,2]
        @fact Images.padarray(A, (0,0), (0,0), "replicate") --> A
        @fact Images.padarray(A, (1,2), (2,0), "replicate") --> [1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4]
        @fact Images.padindexes(A, 1, 1, 0, "circular") --> [2,1,2]
        @fact Images.padindexes(A, 1, 0, 1, "circular") --> [1,2,1]
        @fact Images.padarray(A, [2,1], [0,2], "circular") --> [2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4]
        @fact Images.padindexes(A, 1, 1, 0, "symmetric") --> [1,1,2]
        @fact Images.padindexes(A, 1, 0, 1, "symmetric") --> [1,2,2]
        @fact Images.padarray(A, (1,2), (2,0), "symmetric") --> [2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2]
        @fact Images.padarray(A, (1,2), (2,0), "value", -1) --> [-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1]
        A = [1 2 3; 4 5 6]
        @fact Images.padindexes(A, 2, 1, 0, "reflect") --> [2,1,2,3]
        @fact Images.padindexes(A, 2, 0, 1, "reflect") --> [1,2,3,2]
        @fact Images.padarray(A, (1,2), (2,0), "reflect") --> [6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6]
        A = [1 2; 3 4]
        @fact Images.padarray(A, (1,1)) --> [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4]
        @fact Images.padarray(A, (1,1), "replicate", "both") --> [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4]
        @fact Images.padarray(A, (1,1), "circular", "pre") --> [4 3 4; 2 1 2; 4 3 4]
        @fact Images.padarray(A, (1,1), "symmetric", "post") --> [1 2 2; 3 4 4; 3 4 4]
        A = ["a" "b"; "c" "d"]
        @fact Images.padarray(A, (1,1)) --> ["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"]
        @fact_throws ErrorException Images.padindexes(A, 1, 1, 1, "unknown")
        @fact_throws ErrorException Images.padarray(A, (1,1), "unknown")
        # issue #292
        A = trues(3,3)
        @fact typeof(Images.padarray(A, (1,2), (2,1), "replicate")) --> BitArray{2}
        @fact typeof(Images.padarray(Images.grayim(A), (1,2), (2,1), "replicate")) --> BitArray{2}
    end

    context("Filtering") do
        EPS = 1e-14
        imgcol = Images.colorim(rand(3,5,6))
        imgcolf = convert(Images.Image{RGB{Ufixed8}}, imgcol)
        for T in (Float64, Int)
            A = zeros(T,3,3); A[2,2] = 1
            kern = rand(3,3)
            @fact maximum(abs(Images.imfilter(A, kern) - rot180(kern))) --> less_than(EPS)
            kern = rand(2,3)
            @fact maximum(abs(Images.imfilter(A, kern)[1:2,:] - rot180(kern))) --> less_than(EPS)
            kern = rand(3,2)
            @fact maximum(abs(Images.imfilter(A, kern)[:,1:2] - rot180(kern))) --> less_than(EPS)
        end
        kern = zeros(3,3); kern[2,2] = 1
        @fact maximum(map(abs, imgcol - Images.imfilter(imgcol, kern))) --> less_than(EPS)
        @fact maximum(map(abs, imgcolf - Images.imfilter(imgcolf, kern))) --> less_than(EPS)
        for T in (Float64, Int)
            # Separable kernels
            A = zeros(T,3,3); A[2,2] = 1
            kern = rand(3).*rand(3)'
            @fact maximum(abs(Images.imfilter(A, kern) - rot180(kern))) --> less_than(EPS)
            kern = rand(2).*rand(3)'
            @fact maximum(abs(Images.imfilter(A, kern)[1:2,:] - rot180(kern))) --> less_than(EPS)
            kern = rand(3).*rand(2)'
            @fact maximum(abs(Images.imfilter(A, kern)[:,1:2] - rot180(kern))) --> less_than(EPS)
        end
        A = zeros(3,3); A[2,2] = 1
        kern = rand(3,3)
        @fact maximum(abs(Images.imfilter_fft(A, kern) - rot180(kern))) --> less_than(EPS)
        kern = rand(2,3)
        @fact maximum(abs(Images.imfilter_fft(A, kern)[1:2,:] - rot180(kern))) --> less_than(EPS)
        kern = rand(3,2)
        @fact maximum(abs(Images.imfilter_fft(A, kern)[:,1:2] - rot180(kern))) --> less_than(EPS)
        kern = zeros(3,3); kern[2,2] = 1
        @fact maximum(map(abs, imgcol - Images.imfilter_fft(imgcol, kern))) --> less_than(EPS)
        @fact maximum(map(abs, imgcolf - Images.imfilter_fft(imgcolf, kern))) --> less_than(EPS)

        @fact approx_equal(Images.imfilter(ones(4,4), ones(3,3)), 9.0) --> true
        @fact approx_equal(Images.imfilter(ones(3,3), ones(3,3)), 9.0) --> true
        @fact approx_equal(Images.imfilter(ones(3,3), [1 1 1;1 0.0 1;1 1 1]), 8.0) --> true
        img = convert(Images.Image, ones(4,4))
        @fact approx_equal(Images.imfilter(img, ones(3,3)), 9.0) --> true
        A = zeros(5,5,3); A[3,3,[1,3]] = 1
        @fact Images.colordim(A) --> 3
        kern = rand(3,3)
        kernpad = zeros(5,5); kernpad[2:4,2:4] = kern
        Af = Images.imfilter(A, kern)

        @fact cat(3, rot180(kernpad), zeros(5,5), rot180(kernpad)) --> roughly(Af)
        Aimg = permutedims(convert(Images.Image, A), [3,1,2])
        @fact Images.imfilter(Aimg, kern) --> roughly(permutedims(Af, [3,1,2]))
        @fact approx_equal(Images.imfilter(ones(4,4),ones(1,3),"replicate"), 3.0) --> true

        A = zeros(5,5); A[3,3] = 1
        kern = rand(3,3)
        Af = Images.imfilter(A, kern, "inner")
        @fact Af --> rot180(kern)
        Afft = Images.imfilter_fft(A, kern, "inner")
        @fact Af --> roughly(Afft)
        h = [0.24,0.87]
        hfft = Images.imfilter_fft(eye(3), h, "inner")
        hfft[abs(hfft) .< 3eps()] = 0
        @fact Images.imfilter(eye(3), h, "inner") --> roughly(hfft)  # issue #204

        @fact approx_equal(Images.imfilter_gaussian(ones(4,4), [5,5]), 1.0) --> true
        A = fill(convert(Float32, NaN), 3, 3)
        A[1,1] = 1
        A[2,1] = 2
        A[3,1] = 3
        @fact Images.imfilter_gaussian(A, [0,0]) --> exactly(A)
        @test_approx_eq Images.imfilter_gaussian(A, [0,3]) A
        B = copy(A)
        B[isfinite(B)] = 2
        @test_approx_eq Images.imfilter_gaussian(A, [10^3,0]) B
        @fact maximum(map(abs, Images.imfilter_gaussian(imgcol, [10^3,10^3]) - mean(imgcol))) --> less_than(1e-4)
        @fact maximum(map(abs, Images.imfilter_gaussian(imgcolf, [10^3,10^3]) - mean(imgcolf))) --> less_than(1e-4)
        A = rand(4,5)
        img = reinterpret(Images.Gray{Float64}, Images.grayim(A))
        imgf = Images.imfilter_gaussian(img, [2,2])
        @fact reinterpret(Float64, Images.data(imgf)) --> roughly(Images.imfilter_gaussian(A, [2,2]))
        A = rand(3,4,5)
        img = Images.colorim(A)
        imgf = Images.imfilter_gaussian(img, [2,2])
        @fact reinterpret(Float64, Images.data(imgf)) --> roughly(Images.imfilter_gaussian(A, [0,2,2]))

        A = zeros(Int, 9, 9); A[5, 5] = 1
        @fact maximum(abs(Images.imfilter_LoG(A, [1,1]) - Images.imlog(1.0))) --> less_than(EPS)
        @fact maximum(Images.imfilter_LoG([0 0 0 0 1 0 0 0 0], [1,1]) - sum(Images.imlog(1.0),1)) --> less_than(EPS)
        @fact maximum(Images.imfilter_LoG([0 0 0 0 1 0 0 0 0]', [1,1]) - sum(Images.imlog(1.0),2)) --> less_than(EPS)

        @fact Images.imaverage() --> fill(1/9, 3, 3)
        @fact Images.imaverage([3,3]) --> fill(1/9, 3, 3)
        @fact_throws ErrorException Images.imaverage([5])
    end

    context("Restriction") do
        imgcol = Images.colorim(rand(3,5,6))
        A = reshape([convert(UInt16, i) for i = 1:60], 4, 5, 3)
        B = Images.restrict(A, (1,2))
        Btarget = cat(3, [ 0.96875  4.625   5.96875;
                           2.875   10.5    12.875;
                           1.90625  5.875   6.90625],
                      [ 8.46875  14.625 13.46875;
                        17.875    30.5   27.875;
                        9.40625  15.875 14.40625],
                      [15.96875  24.625 20.96875;
                       32.875    50.5   42.875;
                       16.90625  25.875 21.90625])
        @fact B --> roughly(Btarget)
        Argb = reinterpret(RGB, reinterpret(Ufixed16, permutedims(A, (3,1,2))))
        B = Images.restrict(Argb)
        Bf = permutedims(reinterpret(Float64, B), (2,3,1))
        @fact Bf --> roughly(Btarget/reinterpret(one(Ufixed16)), 1e-12)
        Argba = reinterpret(RGBA{Ufixed16}, reinterpret(Ufixed16, A))
        B = Images.restrict(Argba)
        @fact reinterpret(Float64, B) --> roughly(Images.restrict(A, (2,3))/reinterpret(one(Ufixed16)), 1e-12)
        A = reshape(1:60, 5, 4, 3)
        B = Images.restrict(A, (1,2,3))
        @fact cat(3, [ 2.6015625  8.71875 6.1171875;
                       4.09375   12.875   8.78125;
                       3.5390625 10.59375 7.0546875],
                     [10.1015625 23.71875 13.6171875;
                      14.09375   32.875   18.78125;
                      11.0390625 25.59375 14.5546875]) --> roughly(B)
        Images.restrict(imgcol, (1,2))
    end

    context("Erode/ dilate") do
        A = zeros(4,4,3)
        A[2,2,1] = 0.8
        A[4,4,2] = 0.6
        Ae = Images.erode(A)
        @fact Ae --> zeros(size(A))
        Ad = Images.dilate(A)
        Ar = [0.8 0.8 0.8 0;
              0.8 0.8 0.8 0;
              0.8 0.8 0.8 0;
              0 0 0 0]
        Ag = [0 0 0 0;
              0 0 0 0;
              0 0 0.6 0.6;
              0 0 0.6 0.6]
        @fact Ad --> cat(3, Ar, Ag, zeros(4,4))
        Ae = Images.erode(Ad)
        Ar = [0.8 0.8 0 0;
              0.8 0.8 0 0;
              0 0 0 0;
              0 0 0 0]
        Ag = [0 0 0 0;
              0 0 0 0;
              0 0 0 0;
              0 0 0 0.6]
        @fact Ae --> cat(3, Ar, Ag, zeros(4,4))
        # issue #311
        @fact Images.dilate(trues(3)) --> trues(3)
    end

    context("Opening / closing") do
        A = zeros(4,4,3)
        A[2,2,1] = 0.8
        A[4,4,2] = 0.6
        Ao = Images.opening(A)
        @fact Ao --> zeros(size(A))
        A = zeros(10,10)
        A[4:7,4:7] = 1
        B = copy(A)
        A[5,5] = 0
        Ac = Images.closing(A)
        @fact Ac --> B
    end

    context("Label components") do
        A = [true  true  false true;
             true  false true  true]
        lbltarget = [1 1 0 2;
                     1 0 2 2]
        lbltarget1 = [1 2 0 4;
                      1 0 3 4]
        @fact Images.label_components(A) --> lbltarget
        @fact Images.label_components(A, [1]) --> lbltarget1
        connectivity = [false true  false;
                        true  false true;
                        false true  false]
        @fact Images.label_components(A, connectivity) --> lbltarget
        connectivity = trues(3,3)
        lbltarget2 = [1 1 0 1;
                      1 0 1 1]
        @fact Images.label_components(A, connectivity) --> lbltarget2
    end

    context("Phantoms") do
        P = [ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
              0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
              0.0  0.0  0.2  0.3  0.3  0.2  0.0  0.0;
              0.0  0.0  0.2  0.0  0.2  0.2  0.0  0.0;
              0.0  0.0  0.2  0.0  0.0  0.2  0.0  0.0;
              0.0  0.0  0.2  0.2  0.2  0.2  0.0  0.0;
              0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
              0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 ]
        Q = Images.shepp_logan(8)
        @fact norm((P-Q)[:]) --> less_than(1e-10)
        P = [ 0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0;
              0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
              0.0  0.0  1.02  1.03  1.03  1.02  0.0  0.0;
              0.0  0.0  1.02  1.0   1.02  1.02  0.0  0.0;
              0.0  0.0  1.02  1.0   1.0   1.02  0.0  0.0;
              0.0  0.0  1.02  1.02  1.02  1.02  0.0  0.0;
              0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
              0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0 ]
        Q = Images.shepp_logan(8, highContrast=false)
        @fact norm((P-Q)[:]) --> less_than(1e-10)
    end

    context("Image resize") do
        img = convert(Images.Image, zeros(10,10))
        img2 = Images.imresize(img, (5,5))
        @fact length(img2) --> 25
    end

context("Contrast") do
    # Issue #282
    img = convert(Images.Image{Gray{Ufixed8}}, eye(2,2))
    imgs = Images.imstretch(img, 0.3, 0.4)
    @fact data(imgs) --> roughly(1./(1 + (0.3./(eye(2,2) + eps())).^0.4))
end

end


text:
using Distances
using LightGraphs
using NearestNeighbors
using SimpleWeightedGraphs
using StaticArrays

using Random
using SparseArrays

struct GeodesicComplex{T, P<:AbstractVector, M<:Metric,
                       I<:AbstractVector{Int}, K<:NNTree{P, M}}
    points    ::Vector{P}
    radius    ::T
    metric    ::M
    landmarks ::I
    graph     ::SimpleWeightedGraph{Int, T} # <- rab se samo dists in paths
    tree      ::K
end

function GeodesicComplex(pts::AbstractVector{<:AbstractVector}, r; landmarks = nothing,
                         metric = Euclidean(), tree = KDTree)
    kdt = tree(pts, metric)
    landmark_idxs, cover = getcover(landmarks, pts, r, kdt)
    I = Int[]
    J = Int[]
    T = result_type(metric, pts[1], pts[1])
    V = T[]
    n = length(landmark_idxs)

    for i in 1:n, j in 1:(i - 1)
        if !isempty(intersect(cover[i], cover[j]))
            d = evaluate(metric, pts[landmark_idxs[i]], pts[landmark_idxs[j]])
            append!(I, (i, j))
            append!(J, (j, i))
            append!(V, (d, d))
        end
    end
    graph = SimpleWeightedGraph(sparse(I, J, V, n, n))

    GeodesicComplex(pts, T(r), metric, landmark_idxs, graph, kdt)
end

function getcover(::Nothing, pts, r, tree)
    covered = falses(length(pts))
    idxs = shuffle(eachindex(pts))
    landmarks = Int[]
    cover = Set{Int}[]

    for i in idxs
        covered[i] && continue
        covered[i] = true
        inball = inrange(tree, pts[i], r)
        covered[inball] .= true
        push!(landmarks, i)
        push!(cover, Set(inball))
    end
    ord = sortperm(landmarks)
    landmarks[ord], cover[ord]
end

function getcover(landmarks, pts, r, tree)
    cover = Set{Int}[]
    landmarks = sort(landmarks)

    for l in landmarks
        inball = inrange(tree, pts[l], r)
        push!(cover, Set(inball))
    end
    landmarks, cover
end

function Base.show(io::IO, gc::GeodesicComplex{T, P}) where {T, P}
    print(io, "GeodesicComplex{$T, $P} with $(n_points(gc)) points, " *
          "$(n_landmarks(gc)) landmarks and $(ne(landmark_graph(gc))) edges")
end

# new interface:
landmarks(gc::GeodesicComplex) = gc.landmarks
landmarks(gc::GeodesicComplex, idxs) = gc.landmarks[idxs]

n_landmarks(gc::GeodesicComplex) = length(gc.landmarks)

points(gc::GeodesicComplex) = gc.points
points(gc::GeodesicComplex, idxs) = gc.points[idxs]

n_points(gc::GeodesicComplex) = length(gc.points)

radius(gc::GeodesicComplex) = gc.radius

landmark_graph(gc::GeodesicComplex) = gc.graph

function landmark_shortest_paths(gc)
    fw = floyd_warshall_shortest_paths(landmark_graph(gc))
    fw.dists, fw.parents
end

pairwise_ambient_distance(gc::GeodesicComplex{T, D}, is, js) where {T, D} =
    pairwise(gc.metric,
             reshape(reinterpret(T, points(gc, is)), (length(D), length(is))),
             reshape(reinterpret(T, points(gc, js)), (length(D), length(js))))

nearby_points(gc::GeodesicComplex, i, r = radius(gc)) = inrange(gc.tree, points(gc, i), r)

distance_result_type(gc::GeodesicComplex{T}) where T = T

max_stars_repo_name: mtsch/IntrinsicPersistence

max_issues_repo_name: mtsch/IntrinsicPersistence

max_forks_repo_name: mtsch/IntrinsicPersistence


text:
#=%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Given surface boundary conditions and interior sources of phosphate, produce 3D fields of various tracers.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% =#

using Revise, TMI, GoogleDrive
using PyPlot, PyCall, Test
#, Distributions, LinearAlgebra,  Zygote, ForwardDiff, Optim

TMIversion = "modern_90x45x33_G14_v2"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
qPO₄ = readtracer(TMIfile,"qPO₄")

# PO₄ inversion
PO₄obs = readtracer(TMIfile,"PO₄")
PO₄sfc = PO₄obs[:,:,1] # nutrient surface boundary condition
r = 1.0
PO₄model = steady_inversion(Alu,PO₄sfc,qPO₄,r,γ.wet)
@test isapprox(PO₄model[50,30,10],PO₄obs[50,30,10])

# NO₃ inversion
NO₃obs = readtracer(TMIfile,"NO₃")
NO₃sfc = NO₃obs[:,:,1] # nutrient surface boundary condition
r = 15.5
NO₃model = steady_inversion(Alu,NO₃sfc,qPO₄,r,γ.wet)
@test isapprox(NO₃model[50,30,10],NO₃obs[50,30,10])

# O₂ inversion
O₂obs = readtracer(TMIfile,"O₂")
O₂sfc = O₂obs[:,:,1] # nutrient surface boundary condition
r = -170.0
O₂model = steady_inversion(Alu,O₂sfc,qPO₄,r,γ.wet)
@test isapprox(O₂model[50,30,10],O₂obs[50,30,10])

# δ¹³C inversion
δ¹³Cobs = readtracer(TMIfile,"δ¹³C")
δ¹³Csfc = δ¹³Cobs[:,:,1] # nutrient surface boundary condition
r = -1.1
δ¹³Cmodel = steady_inversion(Alu,δ¹³Csfc,qPO₄,r,γ.wet)
@test isapprox(δ¹³Cmodel[50,30,10],δ¹³Cobs[50,30,10])


text:
import GeoStatsBase: solvesingle

# define a dummy solver for testing
@simsolver Dummy begin end
function solvesingle(problem::SimulationProblem, covars::NamedTuple, ::Dummy, preproc)
  mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))
  reals = map(covars.names) do var
    V = mactypeof[var]
    var => fill(one(V), nelements(domain(problem)))
  end
  Dict(reals)
end


text:
# Computes the function cel(kc,p,a,b) from Bulirsch (1969).
function cel_bulirsch(k2,p,a,b)
@assert (k2 <= 1)
ca = sqrt(eps(k2))
# Avoid undefined k2=1 case:
if k2 != 1.0
  kc = sqrt(1.0-k2)
else
  kc = eps(k2)
end
return cel_bulirsch(k2,kc,p,a,b)
end

# Version called with kc (this is to improve precision of computation):
function cel_bulirsch(k2,kc,p,a,b)
@assert (k2 <= 1)
ca = sqrt(eps(k2))
# Avoid undefined k2=1 case:
if k2 == 1.0 || kc == 0.0
  kc = eps(k2)
end
# Initialize values:
ee = kc; m=1.0
if p > 0.0
  p = sqrt(p); pinv = inv(p); b *= pinv
else
  q=k2; g=1.0-p; f = g-k2
  q *= (b-a*p); ginv = inv(g); p=sqrt(f*ginv); a=(a-b)*ginv
  pinv = inv(p)
  b = -q*ginv^2*pinv+a*p
end
# Compute recursion:
f=a; a += b*pinv; g=ee*pinv; b += f*g; b +=b; p +=g; g=m; m += kc
iter = 0; itmax = 50
while abs(g-kc) > g*ca && iter < itmax
  kc=sqrt(ee)
  kc += kc
  ee = kc*m
  f=a
  pinv = inv(p)
  a += b*pinv
  g=ee*pinv
  b += f*g
  b +=b
  p +=g
  g=m
  m += kc
  iter +=1
end
if iter == itmax
  println("k2 ",k2," kc ",kc," abs(g-kc) ",abs(g-kc)," g*ca ",g*ca," cel ",0.5*pi*(a*m+b)/(m*(m+p)))
end
return 0.5*pi*(a*m+b)/(m*(m+p))
end

# Vector version to attempt to improve speed for multiple elliptic integrals with same k_c:
function cel_bulirsch(k2,kc,p,a1,a2,a3,b1,b2,b3)
# This assumes first value of a and b uses p; the rest have p=1.
@assert (k2 <= 1)
ca = sqrt(eps(k2))
# Avoid undefined k2=1 case:
if k2 == 1.0 || kc == 0.0
  kc = eps(k2)
end
# Initialize values:
ee = kc; m=1.0
if p > 0.0
  p = sqrt(p); pinv = inv(p); b1 *= pinv
else
  q=k2; g=1.0-p; f = g-k2
  q *= (b1-a1*p); ginv = inv(g); p=sqrt(f*ginv); a1=(a1-b1)*ginv
  pinv = inv(p)
  b1 = -q*ginv^2*pinv+a1*p
end
# Compute recursion:
f1=a1
# First compute the first integral with p:
a1 += b1*pinv; g=ee*pinv; b1 += f1*g; b1 += b1; p +=g; g=m;
# Next, compute the remainder with p = 1:
p1 = 1; g1=ee
f2 = a2; f3 = a3
a2 += b2; b2 += f2*g1; b2 += b2
a3 += b3; b3 += f3*g1; b3 += b3
p1 +=g1
g1 = m
m += kc
iter = 0; itmax = 50
while (abs(g-kc) > g*ca) || (abs(g1-kc) > g1*ca)  && iter < itmax
  kc = sqrt(ee)
  kc += kc
  ee = kc*m
  f1 = a1; f2=a2; f3=a3
  pinv = inv(p)
  pinv1 = inv(p1)
  a1 += b1*pinv
  a2 += b2*pinv1
  a3 += b3*pinv1
  g = ee*pinv
  g1= ee*pinv1
  b1 += f1*g
  b2 += f2*g1
  b3 += f3*g1
  b1 += b1
  b2 += b2
  b3 += b3
  p  += g
  p1 += g1
  g  = m
  m += kc
  iter +=1
end
if iter == itmax
  println("k2 ",k2," kc ",kc," abs(g-kc) ",abs(g-kc)," g*ca ",g*ca," cel ",0.5*pi*(a1*m+b1)/(m*(m+p)))
end
f1 = 0.5*pi*(a1*m+b1)/(m*(m+p))
f2 = 0.5*pi*(a2*m+b2)/(m*(m+p1))
f3 = 0.5*pi*(a3*m+b3)/(m*(m+p1))
return f1,f2,f3
end

# More efficient "Vector" version of cel:
function cel_bulirsch_nodiv(k2,kc,p,a1,a2,a3,b1,b2,b3)
# This assumes first value of a and b uses p; the rest have p=1.
@assert (k2 <= 1)
ca = sqrt(eps(k2))
# Avoid undefined k2=1 case:
if k2 == 1.0 || kc == 0.0
  kc = eps(k2)
end
# Initialize values:
ee = kc; m=1.0
if p > 0.0
  d = sqrt(p); a1 *= d; d1 = 1
else
  println("p negative")
  q=k2; g=1.0-p; f = g-k2
  q *= (b1-a1*p); ginv = inv(g); p=sqrt(f*ginv); a1=(a1-b1)*ginv
  pinv = inv(p)
  b1 = -q*ginv^2*pinv+a1*p
end
# Compute recursion:
f1=a1
# First compute the first integral with p:
a1*= p; a1 += b1*d; b1 *=p; b1 += f1*d*ee; b1 *=2;
t = p; p = p^2 + d^2*ee; d *= t
# Next, compute the remainder with p = 1:
p1 = 1
f2 = a2; a2 *= p1; a2 += b2*d1; b2 *= p1; b2 += f2*d1*ee; b2 *=2
f3 = a3; a3 *= p1; a3 += b3*d1; b3 *= p1; b3 += f3*d1*ee; b3 *=2
t1 = p1; p1 = p1^2 + d1^2*ee; d1 *= t1
g = m; m += kc
iter = 0; itmax = 50
while abs(g-kc) > g*ca  && iter < itmax
  kc = sqrt(ee)
  kc += kc
  ee = kc*m
  f1 = a1; f2 = a2; f3 = a3
  a1 *= p;  a1 += b1*d
  a2 *= p1; a2 += b2*d1
  a3 *= p1; a3 += b3*d1
  b1 *= p;  b1 += f1*d*ee
  b2 *= p1; b2 += f2*d1*ee
  b3 *= p1; b3 += f3*d1*ee
  b1 *= 2; b2 *= 2; b3 *= 2
  t = p; t1 = p1
  p = p^2 + d^2*ee
  p1 = p1^2 + d1^2*ee
  d *= t
  d1 *= t1
  g = m
  m += kc
  iter +=1
end
if iter == itmax
  println("k2 ",k2," kc ",kc," abs(g-kc) ",abs(g-kc)," g*ca ",g*ca," cel ",0.5*pi*(a1*m+b1)/(m*(m+p)))
end
println("a1: ",a1," m: ",m," b1: ",b1," d: ",d," p: ",p)
f1 = 0.5*pi*(a1*m+b1)/(m*(d*m+p))
f2 = 0.5*pi*(a2*m+b2)/(m*(d1*m+p1))
f3 = 0.5*pi*(a3*m+b3)/(m*(d1*m+p1))
return f1,f2,f3
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
__precompile__()
"""
Module with reachability algorithms.
"""
module ReachSets

using ..Utils, ..Properties
using LazySets, MathematicalSystems, HybridSystems, Expokit, Optim,
      ProgressMeter

# fix namespace conflicts with MathematicalSystems
using LazySets: LinearMap
using Reachability: info, warn

include("../compat.jl")

using LazySets.Approximations: symmetric_interval_hull,
                               decompose,
                               overapproximate,
                               box_approximation,
                               AbstractDirections

import LazySets.Approximations: project

using Reachability: @timing,
                    Options, OptionSpec, TwoLayerOptions, AbstractOptions,
                    validate_and_wrap_options, print_option_spec,
                    haskey_specified,
                    validate_solver_options_and_add_default_values!

# ========================================
# Discretize and compute bloating factors
# ========================================
include("discretize.jl")

export discretize

# ==========================
# Property checking results
# ==========================
include("CheckSolution.jl")

export CheckSolution

# =============================
# Property checking algorithms
# =============================

# dictionary of registered algorithms
available_algorithms_check = Dict{String, Dict{String, Any}}()

include("ContinuousPost/BFFPSV18/check_blocks.jl")
include("ContinuousPost/BFFPSV18/check_property.jl")
include("ContinuousPost/BFFPSV18/partitions.jl")
include("ContinuousPost/BFFPSV18/inout_map_reach.jl")

# "explicit" backends
push!(available_algorithms_check, "explicit_blocks"=>Dict("func"=>check_blocks,
                                                          "is_explicit"=>true))

export available_algorithms_check,
       check_property

# ====================================================
# Algorithms to find a threshold for property checking
# ====================================================
include("tune.jl")
export tune_δ

# =====================
# Reachability results
# =====================
include("ReachSet.jl")
include("ReachSolution.jl")

export ReachSet,
       ReachSolution

# ===============
# Post operators
# ===============
include("PostOperator.jl")
include("ContinuousPost/ContinuousPost.jl")
include("DiscretePost/DiscretePost.jl")

export PostOperator,
       ContinuousPost,
       DiscretePost,
       init, init!,
       post,
       tube⋂inv!

# ==========================
# Continuous post operators
# ==========================
include("ContinuousPost/BFFPSV18/BFFPSV18.jl")
include("ContinuousPost/BFFPSV18/reach.jl")
include("ContinuousPost/BFFPSV18/reach_blocks.jl")
include("ContinuousPost/BFFPSV18/reach_blocks_wrapping_effect.jl")

include("ContinuousPost/GLGM06/GLGM06.jl")

include("ContinuousPost/TMJets/TMJets.jl")

# ========================
# Reachability Algorithms
# ========================
import Reachability.check_aliases_and_add_default_value!

# dictionary of registered algorithms
available_algorithms = Dict{String, Dict{String, Any}}()

# "explicit" backends
push!(available_algorithms, "explicit_blocks"=>Dict("func"=>reach_blocks!,
                                                    "is_explicit"=>true))

push!(available_algorithms,
      "wrap"=>Dict("func"=>reach_blocks_wrapping_effect!,
                   "is_explicit"=>true))

export available_algorithms

# ==========================
# Discrete post operators
# ==========================
include("DiscretePost/LazyDiscretePost.jl")
include("DiscretePost/ConcreteDiscretePost.jl")

# ==============================================
# Projection of the reach set in two dimensions
# ==============================================
include("project_reach.jl")

export project_reach,
       project

end # module ReachSets


text:
"""
OAuth 2.0 Google Sessions
"""
module session

export GoogleSession, authorize

import Base: string, print, show
using Dates

using HTTP
using Base64
import JSON
import MbedTLS

using ..error
using ..credentials
using ..root

"""
    SHA256withRSA(message, key)

Sign message using private key with RSASSA-PKCS1-V1_5-SIGN algorithm.
"""
SHA256withRSA(message, key::MbedTLS.PKContext) = MbedTLS.sign(key, MbedTLS.MD_SHA256,
    MbedTLS.digest(MbedTLS.MD_SHA256, message), MbedTLS.MersenneTwister(0)
)

"""
GoogleSession(...)

OAuth 2.0 session for Google using provided credentials.

Caches authorisation tokens up to expiry.

```julia
sess = GoogleSession(JSONCredentials(expanduser("~/auth.json")), ["devstorage.full_control"])
```
"""
mutable struct GoogleSession{T <: Credentials}
    credentials::T
    scopes::Vector{String}
    authorization::Dict{Symbol, Any}
    expiry::DateTime
    """
        GoogleSession(credentials, scopes)

    Set up session with credentials and OAuth scopes.
    """
    function GoogleSession(credentials::T, scopes::AbstractVector{<: AbstractString}) where {T <: Credentials}
        scopes = [isurl(scope) ? scope : "$SCOPE_ROOT/$scope" for scope in scopes]
        new{T}(credentials, scopes, Dict{String, String}(), DateTime(1))
    end
end
function GoogleSession(credentials::Union{AbstractString, Nothing},
                       scopes::AbstractVector{<: AbstractString}=String[])
    if credentials === nothing
        credentials = ""
    end
    if isfile(credentials)
        credentials = JSONCredentials(credentials)
    else
        if isempty(credentials)
            credentials = MetadataCredentials()
        elseif isurl(credentials)
            credentials = MetadataCredentials(; url=credentials)
        else
            throw(CredentialError("Invalid credentials: ", credentials))
        end
        scopes = intersect(scopes, credentials.scopes)
    end
    GoogleSession(credentials, scopes)
end
function GoogleSession(scopes::AbstractVector{<: AbstractString})
    GoogleSession(get(ENV, "GOOGLE_APPLICATION_CREDENTIALS", ""), scopes)
end

function print(io::IO, x::GoogleSession)
    println(io, "scopes: $(x.scopes)")
    println(io, "authorization: $(x.authorization)")
    println(io, "expiry: $(x.expiry)")
end
show(io::IO, x::GoogleSession) = print(io, x)

"""
    unixseconds(x)

Convert date-time into unix seconds.
"""
unixseconds(x::DateTime) = trunc(Int, datetime2unix(x))

"""
    JWTHeader

JSON Web Token header.
"""
struct JWTHeader
    algorithm::String
end
function string(x::JWTHeader)
    base64encode(JSON.json(Dict(:alg => x.algorithm, :typ => "JWT")))
end
print(io::IO, x::JWTHeader) = print(io, string(x))

"""
    JWTClaimSet

JSON Web Token claim-set.
"""
struct JWTClaimSet
    issuer::String
    scopes::Vector{String}
    assertion::DateTime
    expiry::DateTime
    function JWTClaimSet(issuer::AbstractString, scopes::AbstractVector{<: AbstractString},
                         assertion::DateTime=now(UTC), expiry::DateTime=now(UTC) + Hour(1))
        new(issuer, scopes, assertion, expiry)
    end
end
function string(x::JWTClaimSet)
    base64encode(JSON.json(Dict(
        :iss => x.issuer, :scope => join(x.scopes, " "), :aud => AUD_ROOT,
        :iat => unixseconds(x.assertion),
        :exp => unixseconds(x.expiry)
    )))
end
print(io::IO, x::JWTClaimSet) = print(io, string(x))

"""
    JWS(credentials, claimset)

Construct the Base64-encoded JSON Web Signature based on the JWT header, claimset
and signed using the private key provided in the Google JSON service-account key.
"""
function JWS(credentials::JSONCredentials, claimset::JWTClaimSet, header::JWTHeader=JWTHeader("RS256"))
    payload = "$header.$claimset"
    signature = base64encode(SHA256withRSA(payload, credentials.private_key))
    "$payload.$signature"
end

function token(credentials::JSONCredentials, scopes::AbstractVector{<: AbstractString})
    # construct claim-set from service account email and requested scopes
    claimset = JWTClaimSet(credentials.client_email, scopes)
    data = HTTP.URIs.escapeuri(Dict{Symbol, String}(
        :grant_type => "urn:ietf:params:oauth:grant-type:jwt-bearer",
        :assertion => JWS(credentials, claimset)
    ))
    headers = Dict{String, String}("Content-Type" => "application/x-www-form-urlencoded")
    res = HTTP.post("$AUD_ROOT", headers, data)
    if res.status != 200
        throw(SessionError("Unable to obtain authorization: $(String(res.body))"))
    end
    authorization = JSON.parse(String(res.body); dicttype=Dict{Symbol, Any})
    authorization, claimset.assertion
end

function token(credentials::MetadataCredentials, ::AbstractVector{<: AbstractString})
    assertion = now(UTC)
    authorization = JSON.parse(get(credentials, "token"); dicttype=Dict{Symbol, Any})
    authorization, assertion
end

"""
    authorize(session; cache=true)

Get OAuth 2.0 authorisation token from Google.

If `cache` set to `true`, get a new token only if the existing token has expired.
"""
function authorize(session::GoogleSession; cache::Bool=true)
    # don't get a new token if a non-expired one exists
    if cache && (session.expiry >= now(UTC)) && !isempty(session.authorization)
        return session.authorization
    end

    authorization, assertion = try token(session.credentials, session.scopes) catch e
        session.expiry = DateTime(1)
        empty!(session.authorization)
        rethrow(e)
    end

    # cache authorization if required
    if cache
        session.expiry = assertion + Second(authorization[:expires_in] - 30)
        session.authorization = authorization
    else
        session.expiry = DateTime(1)
        empty!(session.authorization)
    end
    authorization
end

end


text:
module VisualParticipationAnalyticsServers

using DataFrames
using HTTP
using JSON3
using JSONTables
using SQLite
using Tables

export start,
    startrest,
    startws

include("rest.jl")
include("websocket.jl")

"""
    start()
"""
function start()
    startrest()
    startws()
end

end


text:
module pigpiod_if2

using Libdl
using pigpiod_if2_jll

include("./constants.jl")
include("./i2c.jl")
include("./serial.jl")
include("./spi.jl")
include("./waves.jl")
export INPUT, OUTPUT
export PUD_UP, PUD_OFF, PUD_DOWN
export HIGH, LOW

export pigpio_start, pigpio_stop, set_mode, gpio_write, get_mode, set_pull_up_down, gpio_read

# ESSENTIAL
function pigpio_start(addr::String = "", port::String = "")
    return ccall((:pigpio_start, libpigpiod_if2), Int32, (Cstring, Cstring), addr, port)
end

function pigpio_stop(p::Int32)
    return ccall((:pigpio_stop, libpigpiod_if2), Cvoid, (Int32,), p)
end

# BASIC
function set_mode(p::Int32, gpio::Int, mode::Int)
    return ccall((:set_mode, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, mode)
end

function get_mode(p::Int32, gpio::Int)
    return ccall((:get_mode, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

function set_pull_up_down(p::Int32, gpio::Int, pud::Int)
    return ccall((:set_pull_up_down, libpigpiod_if2), Int32, (Int32,  UInt32, UInt32), p, gpio, pud)
end

function gpio_read(p::Int32, gpio::Int)
    return ccall((:gpio_read, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

function gpio_write(p::Int32, gpio::Int, level::Int)
    return ccall((:gpio_write, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, level)
end

# PWM
export set_PWM_range, set_PWM_frequency, set_PWM_dutycycle, get_PWM_range, get_PWM_dutycycle, get_PWM_frequency, get_PWM_real_range
# set_PWM_dutycycle
function set_PWM_dutycycle(p::Int32, gpio::Int, dutycycle::Int)
    return ccall((:set_PWM_dutycycle, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, dutycycle)
end

# set_PWM_frequency
function set_PWM_frequency(p::Int32, gpio::Int, frequency::Int)
    return ccall((:set_PWM_frequency, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, frequency)
end

# set_PWM_range
function set_PWM_range(p::Int32, gpio::Int, r::Int)
    return ccall((:set_PWM_range, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, r)
end

# get_PWM_dutycycle
function get_PWM_dutycycle(p::Int32, gpio::Int)
    return ccall((:get_PWM_dutycycle, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

# get_PWM_frequency
function get_PWM_frequency(p::Int32, gpio::Int)
    return ccall((:get_PWM_frequency, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

# get_PWM_range
function get_PWM_range(p::Int32, gpio::Int)
    return ccall((:get_PWM_range, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

# get_PWM_real_range
function get_PWM_real_range(p::Int32, gpio::Int)
    return ccall((:get_PWM_real_range, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

# Servo
export set_servo_pulsewidth, get_servo_pulsewidth
# set_servo_pulsewidth
function set_servo_pulsewidth(p::Int32, gpio::Int, pulsewidth::Int)
    return ccall((:set_servo_pulsewidth, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, pulsewidth)
end

# get_servo_pulsewidth
function get_servo_pulsewidth(p::Int32, gpio::Int)
    return ccall((:get_servo_pulsewidth, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

# INTERMEDIATE
export gpio_trigger, set_watchdog, read_bank_1, read_bank_2, clear_bank_1, clear_bank_2, set_bank_1, set_bank_2, callback, callback_ex, callback_cancel, wait_for_edge
# gpio_trigger
function gpio_trigger(p::Int32, gpio::Int, pulseLen::Int, level::Int)
    return ccall((:gpio_trigger, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, gpio, pulseLen, level)
end

# set_watchdog
function set_watchdog(p::Int32, gpio::Int, timeout::Int)
    return ccall((:set_watchdog, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, timeout)
end

# read_bank_1
function read_bank_1(p::Int32)
    return ccall((:read_bank_1, libpigpiod_if2), UInt32, (Int32,), p)
end

# read_bank_2
function read_bank_2(p::Int32)
    return ccall((:read_bank_2, libpigpiod_if2), UInt32, (Int32,), p)
end

# clear_bank_1
function clear_bank_1(p::Int32, bits::UInt32)
    return ccall((:clear_bank_1, libpigpiod_if2), Int32, (Int32, UInt32), p, bits)
end

# set_bank_1
function set_bank_1(p::Int32, bits::UInt32)
    return ccall((:set_bank_1, libpigpiod_if2), Int32, (Int32, UInt32), p, bits)
end

# set_bank_2
function set_bank_2(p::Int32, bits::UInt32)
    return ccall((:set_bank_2, libpigpiod_if2), Int32, (Int32, UInt32), p, bits)
end

# clear_bank_2
function clear_bank_2(p::Int32, bits::UInt32)
    return ccall((:clear_bank_2, libpigpiod_if2), Int32, (Int32, UInt32), p, bits)
end

# callback
function callback(p::Int32, gpio::Int, edge::Int, f::Ptr{Nothing})
    return ccall((:callback, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ptr{Cvoid}), p, gpio, edge, f)
end

# callback_ex
function callback_ex(p::Int32, gpio::Int, edge::Int, f::Ptr{Nothing}, user_data)
    return ccall((:callback_ex, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ptr{Cvoid}, Any), p, gpio, edge, f, user_data)
end

# callback_cancel
function callback_cancel(callback_id::Int32)
    return ccall((:callback_cancel, libpigpiod_if2), Int32, (UInt32,), callback_id)
end

# wait_for_edge
function wait_for_edge(p::Int32, gpio::Int, edge::Int, timeout::Float32)
    return ccall((:wait_for_edge, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Float32), p, gpio, edge, timeout)
end

# Let's not do threading for now, since it depends on "pthread.h"
# start_thread
# stop_thread

# ADVANCED
# notify_open
function notify_open(p::Int32)
    return ccall((:notify_open, libpigpiod_if2), Int32, (Int32,), p)
end

# notify_begin
function notify_begin(p::Int32, handle::Int32, bits::UInt32)
    return ccall((:notify_begin, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, handle, bits)
end

# notify_pause
function notify_pause(p::Int32, handle::Int32)
    return ccall((:notify_pause, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

# notify_close
function notify_close(p::Int32, handle::Int32)
    return ccall((:notify_close, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

# hardware_clock
function hardware_clock(p::Int32, gpio::Int, clkfreq::Int)
    ccall((:hardware_clock, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, clkfreq)
end

# hardware_PWM
function hardware_PWM(p::Int32, gpio::Int, PWMfreq::Int, PWMduty::Int)
    ccall((:hardware_PWM, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, gpio, PWMfreq, PWMduty)
end

# set_glitch_filter
function set_glitch_filter(p::Int32, gpio::Int, steady::Int)
    return ccall((:set_glitch_filter, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, steady)
end

# set_noise_filter
function set_noise_filter(p::Int32, gpio::Int, steady::Int, active::Int)
    return ccall((:set_noise_filter, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, gpio, steady, active)
end

# set_pad_strength
function set_pad_strength(p::Int32, pad::Int, padStrength::Int)
    return ccall((:set_pad_strength, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, pad, padStrength)
end

# get_pad_strength
function get_pad_strength(p::Int32, pad::Int)
    return ccall((:get_pad_strength, libpigpiod_if2), Int32, (Int32, UInt32), p, pad)
end

# shell_ # Not necessary
# Custom (skipping)
# Events
export event_callback, event_callback_ex, event_callback_cancel, event_trigger, wait_for_event
# event_callback
function event_callback(p::Int32, event::Int, f)
    return ccall((:event_callback, libpigpiod_if2), Int32, (Int32, UInt32, Ptr{Cvoid}), p, event, f)
end

# event_callback_ex
function event_callback_ex(p::Int32, event::Int, f, user_data)
    return ccall((:event_callback_ex, libpigpiod_if2), Int32, (Int32, UInt32, Ptr{Cvoid}, Any), p, event, f, user_data)
end

# event_callback_cancel
function event_callback_cancel(callback_id::Int)
    return ccall((:event_callback_cancel, libpigpiod_if2), Int32, (UInt32,), callback_id)
end

# event_trigger
function event_trigger(p::Int32, event::Int)
    return ccall((:event_trigger, libpigpiod_if2), Int32, (Int32, UInt32), p, event)
end

# wait_for_event
function wait_for_event(p::Int32, event::Int, timeout::Float32)
    return ccall((:wait_for_event, libpigpiod_if2), Int32, (Int32, UInt32, Float32), p, event, timeout)
end

# Scripts
export store_script, run_script, update_script, script_status, stop_script, delete_script
# store_script
# int store_script(int pi, char *script)
function store_script(p::Int32, script::String)
    return ccall((:store_script, libpigpiod_if2), Int32, (Int32, Cstring), p, script)
end

# run_script
# int run_script(int pi, unsigned script id, unsigned numPar, uint32_t *param)
function run_script(p::Int32, script_id::Int32, numPar::Int, param::Array{Int})
    return ccall((:run_script, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ptr{Cint}), p, script_id, numPar, param)
end

# update_script
# int update_script(int pi, unsigned script_id, unsigned numPar, uint32_t *param)
function update_script(p::Int32, script_id::Int32, numPar::Int, param::Array{Int})
    return ccall((:update_script, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ptr{Cint}), p, script_id, numPar, param)
end

# script_status
# int script_status(int pi, unsigned script_id, uint32_t *param)
function script_status(p::Int32, script_id::Int32, param::Array{Int})
    return ccall((:script_status, libpigpiod_if2), Int32, (Int32, UInt32, Ptr{Cint}), p, script_id, param)
end

# stop_script
# int stop_script(int pi, unsigned script_id)
function stop_script(p::Int32, script_id::Int32)
    return ccall((:stop_script, libpigpiod_if2), Int32, (Int32, UInt32), p, script_id)
end

# delete_script
# int delete_script(int pi, unsigned script_id)
function delete_script(p::Int32, script_id::Int32)
    return ccall((:delete_script, libpigpiod_if2), Int32, (Int32, UInt32), p, script_id)
end

# FILES
#  file_open	Opens a file
#  int file_open(int pi, char *file, unsigned mode)
#  This function returns a handle to a file opened in a specified mode.
function file_open(p::Int32, file::String, mode::Int)
    return ccall((:file_open, libpigpiod_if2), Int32, (Int32, Cstring, UInt32), p, file, mode)
end

#  file_close	Closes a file
#  int file_close(int pi, unsigned handle)
#  This function closes the file associated with handle.
function file_close(p::Int32, handle::Int32)
    return ccall((:file_close, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

#  file_read	Reads bytes from a file
#  int file_read(int pi, unsigned handle, char *buf, unsigned count)
#  This function reads up to count bytes from the the file associated with handle and writes them to buf.
function file_read(p::Int32, handle::Int32, buf::Vector{Int}, count::Int)
    return ccall((:file_read, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32), p, handle, buf, count)
end

#  file_write	Writes bytes to a file
#  int file_write(int pi, unsigned handle, char *buf, unsigned count)
#  This function writes count bytes from buf to the the file associated with handle.
#  Figure out if I need to change the vector type to a character instead of an Int
function file_write(p::Int32, handle::Int32, buf::Vector{Int}, count::Int)
    return ccall((:file_write, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32), p, handle, buf, count)
end

#  file_seek	Seeks to a position within a file
#  int file_seek(int pi, unsigned handle, int32_t seekOffset, int seekFrom)
#  This function seeks to a position within the file associated with handle.
function file_seek(p::Int32, handle::Int32, seekOffset::Int, seekFrom::Int)
    return ccall((:file_seek, libpigpiod_if2), Int32, (Int32, UInt32, Int32, Int32), p, handle, seekOffset, seekFrom)
end

#  file_list	List files which match a pattern
#  int file_list(int pi, char *fpat, char *buf, unsigned count)
#  This function returns a list of files which match a pattern.
function file_list(p::Int32, fpat::String, buf::String, count::Int)
    return ccall((:file_list, libpigpiod_if2), Int32, (Int32, Cstring, Cstring, UInt32), p, fpat, buf, count)
end

# UTILITIES
export time_time, time_sleep, pigpio_error
function time_time()
    return ccall((:time_time, libpigpiod_if2), Cdouble, ())
end

function time_sleep(seconds::Float64)
    return ccall((:time_sleep, libpigpiod_if2), Cvoid, (Cdouble,), seconds)
end

function pigpio_error(errnum::Int)
    error_str = ccall((:pigpio_error, libpigpiod_if2), Cstring, (Int32,), errnum)
    return unsafe_string(error_str)
end

end # module


text:
using Dates: DateTime
using Registrator.WebUI: isauthorized, AuthFailure, AuthSuccess, User
using GitForge: GitForge, GitHub, GitLab, Bitbucket
using HTTP: stacktrace

using Mocking

Mocking.activate()

function patch_gitforge(body::Function; is_collaborator=false, is_member=false)
    patches = [
        @patch GitForge.is_collaborator(args...) = is_collaborator
        @patch GitForge.is_member(args...) = is_member
    ]

    apply(patches) do
        return body()
    end
end

@testset "gitutils" begin

@testset "isauthorized" begin
    @test isauthorized("username", "reponame") == AuthFailure("Unkown user type or repo type")

    @testset "GitHub" begin

        user = GitHub.User(login="user123")
        org = GitHub.User(login="JuliaLang")
        private_repo = GitHub.Repo(name="Example.jl", private=true, owner=user)
        public_repo_of_user = GitHub.Repo(name="Example.jl", private=false, owner=user, organization=nothing)
        public_repo_of_org = GitHub.Repo(name="Example.jl", private=false, owner=org, organization=org)
        u = User(user, GitHub.GitHubAPI())

        @testset "private repo" begin
            # Assuming CONFIG["allow_private"] is false
            @test isauthorized(u, private_repo) == AuthFailure("Repo Example.jl is private")
        end

        @testset "public repo of user" begin
            # authorized if user is a collaborator on the repo
            patch_gitforge(is_collaborator=true) do
                @test isauthorized(u, public_repo_of_user) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false) do
                @test isauthorized(u, public_repo_of_user) == AuthFailure("User user123 is not a collaborator on repo Example.jl")
            end
        end

        @testset "public repo of org" begin
            # authorized if user is either a collaborator on the repo or member of the org
            patch_gitforge(is_collaborator=true, is_member=true) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=true, is_member=false) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=true) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=false) do
                @test isauthorized(u, public_repo_of_org) == AuthFailure("User user123 is not a member of the org JuliaLang and not a collaborator on repo Example.jl")
            end
        end
    end

    @testset "GitLab" begin

        user = GitLab.User(name="user123", username="user123", id=111)
        org = GitLab.User(name="org123", username="org123", id=222)
        private_project = GitLab.Project(name="Example.jl", visibility="private", owner=user)
        public_project_of_user = GitLab.Project(name="Example.jl", visibility="public", owner=user, namespace=GitLab.Namespace(kind="user"))
        public_project_of_group = GitLab.Project(name="Example.jl", visibility="public", owner=org, namespace=GitLab.Namespace(kind="group", full_path="org123/subgroup/Example.jl"))
        u = User(user, GitLab.GitLabAPI())

        @testset "private project" begin
            # Assuming CONFIG["allow_private"] is false
            @test isauthorized(u, private_project) == AuthFailure("Project Example.jl is private")
        end

        @testset "public project of user" begin
            # authorized if user is a collaborator on the project
            patch_gitforge(is_collaborator=true) do
                @test isauthorized(u, public_project_of_user) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false) do
                @test isauthorized(u, public_project_of_user) == AuthFailure("User user123 is not a member of project Example.jl")
            end
        end

        @testset "public project of group" begin
            # authorized if user is a collaborator on the project or member of the group/subgroups
            patch_gitforge(is_collaborator=true, is_member=true) do
                @test isauthorized(u, public_project_of_group) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=true) do
                @test isauthorized(u, public_project_of_group) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=true, is_member=false) do
                @test isauthorized(u, public_project_of_group) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=false) do
                @test isauthorized(u, public_project_of_group) == AuthFailure("Project Example.jl belongs to the group org123/subgroup/Example.jl, and user user123 is not a member of that group or its parent group(s)")
            end
        end
    end

    @testset "Bitbucket" begin

        #user = @gf get_user(UI.PROVIDERS["bitbucket"].client, "wrburdick")
        user = Bitbucket.User(; nickname = "wrburdick")
        uworkspace = Bitbucket.Workspace(slug="wrb-julia-test")
        workspace = Bitbucket.Workspace(slug="wrburdick")
        private_repo = Bitbucket.Repo(; slug="Example.jl", is_private=true, owner=user, uworkspace)
        public_repo_of_user = Bitbucket.Repo(; slug="Example.jl", is_private=false, owner=user, workspace)
        public_repo_of_org = Bitbucket.Repo(; slug="Example.jl", is_private=false, owner=user, workspace)
        u = User(user, Bitbucket.BitbucketAPI())

        @testset "private repo" begin
            # Assuming CONFIG["allow_private"] is false
            @test isauthorized(u, private_repo) == AuthFailure("Repo $(private_repo.slug) is private")
        end

        @testset "public repo of user" begin
            # authorized if user is a collaborator on the repo
            patch_gitforge(is_collaborator=true) do
                @test isauthorized(u, public_repo_of_user) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false) do
                @test isauthorized(u, public_repo_of_user) == AuthFailure("User $(user.nickname) is not a member of the workspace $(workspace.slug) or a collaborator on repo $(public_repo_of_user.slug)")
            end
        end

        @testset "public repo of org" begin
            # authorized if user is either a collaborator on the repo or member of the org
            patch_gitforge(is_collaborator=true, is_member=true) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=true, is_member=false) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=true) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=false) do
                @test isauthorized(u, public_repo_of_org) == AuthFailure("User $(user.nickname) is not a member of the workspace $(workspace.slug) or a collaborator on repo $(public_repo_of_org.slug)")
            end
        end
    end
end

end