text:
"""
Abstract stiffness topology optimization problem. All subtypes must have the following fields:
    ch::ConstraintHandler
    black::BitVector
    white::BitVector
    varind::AbstractVector{Int}
    metadata::Metadata
"""
abstract type StiffnessTopOptProblem{dim, T} <: AbstractTopOptProblem end

whichdevice(p::StiffnessTopOptProblem) = whichdevice(p.ch)
whichdevice(ch::ConstraintHandler) = whichdevice(ch.dh)
whichdevice(dh::DofHandler) = whichdevice(dh.grid)
whichdevice(g::JuAFEM.Grid) = whichdevice(g.cells)

# Fallbacks
getdim(::StiffnessTopOptProblem{dim, T}) where {dim, T} = dim
floattype(::StiffnessTopOptProblem{dim, T}) where {dim, T} = T
getE(p::StiffnessTopOptProblem) = p.E
getν(p::StiffnessTopOptProblem) = p.ν
getgeomorder(p::StiffnessTopOptProblem) = nnodespercell(p) == 9 ? 2 : 1
getdensity(::StiffnessTopOptProblem{dim, T}) where {dim, T} = T(0)
getmetadata(p::StiffnessTopOptProblem) = p.metadata
getdh(p::StiffnessTopOptProblem) = p.ch.dh
getcloaddict(p::StiffnessTopOptProblem{dim, T}) where {dim, T} = Dict{String, Vector{T}}()
getpressuredict(p::StiffnessTopOptProblem{dim, T}) where {dim, T} = Dict{String, T}()
getfacesets(p::StiffnessTopOptProblem{dim, T}) where {dim, T} = Dict{String, Tuple{Int, T}}()
JuAFEM.getncells(problem::StiffnessTopOptProblem) = JuAFEM.getncells(getdh(problem).grid)

"""
Stiffness problem imported from a .inp file.
"""
struct InpStiffness{dim, N, TF, M, TI, TBool, GO, TInds<:AbstractVector{TI}, TMeta<:Metadata} <: StiffnessTopOptProblem{dim, TF}
    inp_content::InpParser.InpContent{dim, TF, N, TI}
    geom_order::Type{Val{GO}}
    ch::ConstraintHandler{DofHandler{dim, N, TF, M}, TF}
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end

"""
Imports stiffness problem from a .inp file.
"""
function InpStiffness(filepath_with_ext::AbstractString)
    problem = InpParser.extract_inp(filepath_with_ext)
    return InpStiffness(problem)
end
function InpStiffness(problem::InpParser.InpContent)
    ch = InpParser.inp_to_juafem(problem)
    black, white = find_black_and_white(ch.dh)
    varind = find_varind(black, white)
    metadata = Metadata(ch.dh)
    geom_order = JuAFEM.getorder(ch.dh.field_interpolations[1])
    return InpStiffness(problem, Val{geom_order}, ch, black, white, varind, metadata)
end

getE(p::InpStiffness) = p.inp_content.E
getν(p::InpStiffness) = p..inp_content.ν
nnodespercell(::InpStiffness{dim, N}) where {dim, N} = N
getgeomorder(p::InpStiffness{dim, N, TF, M, TI, GO}) where {dim, N, TF, M, TI, GO} = GO
getdensity(p::InpStiffness) = p.inp_content.density
getpressuredict(p::InpStiffness) = p.inp_content.dloads
getcloaddict(p::InpStiffness) = p.inp_content.cloads
getfacesets(p::InpStiffness) = p.inp_content.facesets

"""
```
///**********************************
///*                                *
///*                                * |
///*                                * |
///********************************** v

struct PointLoadCantilever{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{DofHandler{dim, N, T, M}, T}
    force::T
    force_dof::Int
    metadata::Metadata
end
```

`dim`: dimension of the problem

`T`: number type for computations and coordinates

`N`: number of nodes in a cell of the grid

`M`: number of faces in a cell of the grid

`rect_grid`: a RectilinearGrid struct

`E`: Young's modulus

`ν`: Poisson's ration

`ch`: a JuAFEM.ConstraintHandler struct

`force`: force at the center right of the cantilever beam (positive is downward)

`force_dof`: dof number at which the force is applied

`metadata`:: Metadata having various cell-node-dof relationships

API:
```
    PointLoadCantilever(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
```

Example:
```

nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = 1.0;
problem = PointLoadCantilever(nels, sizes, E, ν, force)
```
"""
struct PointLoadCantilever{dim, T, N, M, TBool, TInds<:AbstractVector{Int}, TG <: RectilinearGrid{dim, T, N, M}, TMeta<:Metadata, CH<:ConstraintHandler{<:DofHandler{dim, N, T, M}, T}} <: StiffnessTopOptProblem{dim, T}
    rect_grid::TG
    E::T
    ν::T
    ch::CH
    force::T
    force_dof::Int
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end

function PointLoadCantilever(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim, CellType}
    iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) || throw("Grid does not have an even number of elements along the y and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(rect_grid.grid, "down_force", x -> right(rect_grid, x) && middley(rect_grid, x));

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim) # Add a displacement field
    else
        ip = Lagrange{2, RefCube, 2}()
        push!(dh, :u, dim, ip) # Add a displacement field
    end
    close!(dh)

    ch = ConstraintHandler(dh)

    #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    dbc = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    add!(ch, dbc)
    close!(ch)
    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return PointLoadCantilever(rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata)
end

"""
```
 |
 |
 v
O*********************************
O*                               *
O*                               *
O*                               *
O*********************************
                                 O

struct HalfMBB{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{DofHandler{dim, N, T, M}, T}
    force::T
    force_dof::Int
    metadata::Metadata
end
```

`dim`: dimension of the problem

`T`: number type for computations and coordinates

`N`: number of nodes in a cell of the grid

`M`: number of faces in a cell of the grid

`rect_grid`: a RectilinearGrid struct

`E`: Young's modulus

`ν`: Poisson's ration

`ch`: a JuAFEM.ConstraintHandler struct

`force`: force at the top left of half the MBB (positive is downward)

`force_dof`: dof number at which the force is applied

`metadata`:: Metadata having various cell-node-dof relationships

API:
```
    HalfMBB(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
```

Example:
```

nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = -1.0;
problem = HalfMBB(nels, sizes, E, ν, force)
```
"""
struct HalfMBB{dim, T, N, M, TBool, TInds<:AbstractVector{Int}, TG<:RectilinearGrid{dim, T, N, M}, TMeta<:Metadata, CH<:ConstraintHandler{<:DofHandler{dim, N, T, M}, T}} <: StiffnessTopOptProblem{dim, T}
    rect_grid::TG
    E::T
    ν::T
    ch::CH
    force::T
    force_dof::Int
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end
function HalfMBB(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim, CellType}
    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_u1")
        pop!(rect_grid.grid.facesets, "fixed_u1")
    end
    #addfaceset!(rect_grid.grid, "fixed_u1", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_u1", x -> left(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "fixed_u2")
        pop!(rect_grid.grid.nodesets, "fixed_u2")
    end
    addnodeset!(rect_grid.grid, "fixed_u2", x -> bottom(rect_grid, x) && right(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(rect_grid.grid, "down_force", x -> top(rect_grid, x) && left(rect_grid, x));

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim)
    else
        ip = Lagrange{2, RefCube, 2}()
        push!(dh, :u, dim, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    #dbc1 = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_u1"), (x,t)->T[0], [1])
    dbc1 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_u1"), (x,t)->T[0], [1])
    add!(ch, dbc1)
    dbc2 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_u2"), (x,t)->T[0], [2])
    add!(ch, dbc2)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return HalfMBB(rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata)
end

nnodespercell(p::Union{PointLoadCantilever, HalfMBB}) = nnodespercell(p.rect_grid)
function getcloaddict(p::Union{PointLoadCantilever{dim, T}, HalfMBB{dim, T}}) where {dim, T}
    f = T[0, -p.force, 0]
    fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1]
    return Dict{Int, Vector{T}}(fnode => f)
end

"""
```
////////////
............
.          .
.          .
.          .
.          .
.          ......................
.                               .
.                               .
.                               . |
................................. v
                                force
```
"""
struct LBeam{T, N, M, TBool, TInds<:AbstractVector{Int}, TMeta<:Metadata, CH<:ConstraintHandler{<:DofHandler{2, N, T, M}, T}} <: StiffnessTopOptProblem{2, T}
    E::T
    ν::T
    ch::CH
    force::T
    force_dof::Int
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end
function LBeam(::Type{Val{CellType}}, ::Type{T}=Float64; length = 100, height = 100, upperslab = 50, lowerslab = 50, E = 1.0, ν = 0.3, force = 1.0) where {T, CellType}
    # Create displacement field u
    grid = LGrid(Val{CellType}, T, length=length, height=height, upperslab=upperslab,
        lowerslab=lowerslab)

    dh = DofHandler(grid)
    if CellType === :Linear
        push!(dh, :u, 2)
    else
        ip = Lagrange{2, RefCube, 2}()
        push!(dh, :u, 2, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc = Dirichlet(:u, getfaceset(grid, "top"), (x,t)->T[0, 0], [1, 2])
    add!(ch, dbc)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(grid, "load"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    TInds = typeof(varind)
    TMeta = typeof(metadata)
    return LBeam(E, ν, ch, force, force_dof, black, white, varind, metadata)
end

function boundingbox(grid::JuAFEM.Grid{dim}) where dim
    xmin1 = minimum(n->n.x[1], grid.nodes)
    xmax1 = maximum(n->n.x[1], grid.nodes)
    xmin2 = minimum(n->n.x[2], grid.nodes)
    xmax2 = maximum(n->n.x[2], grid.nodes)
    if dim === 2
        return ((xmin1, xmin2), (xmax1, xmax2))
    else
        xmin3 = minimum(n->n.x[3], grid.nodes)
        xmax3 = maximum(n->n.x[3], grid.nodes)
        return ((xmin1, xmin2, xmin3), (xmax1, xmax2, xmax3))
    end
end

function RectilinearTopology(b, topology = ones(getncells(getdh(b).grid)))
    bb = boundingbox(getdh(b).grid)
    go = getgeomorder(b)
    nels = Int.(round.(bb[2] .- bb[1]))
    dim = length(nels)
    if go === 1
        rectgrid = generate_grid(Quadrilateral, nels, Vec{dim}(bb[1]), Vec{dim}(bb[2]))
    elseif go === 2
        rectgrid = generate_grid(QuadraticQuadrilateral, nels, Vec{dim}(bb[1]), Vec{dim}(bb[2]))
    else
        throw("Unsupported geometry.")
    end
    new_topology = zeros(prod(nels))
    for (i, cell) in enumerate(CellIterator(getdh(b).grid))
        sub = Int.(round.((cell.coords[1]...,))) .+ (1, 1)
        ind = sub2ind(nels, sub...,)
        new_topology[ind] = topology[i]
    end
    return reshape(new_topology, nels)'
end

nnodespercell(p::LBeam{T, N}) where {T, N} = N
getdim(::LBeam) = 2
function getcloaddict(p::LBeam{T}) where {T}
    f = T[0, -p.force]
    fnode = Tuple(getnodeset(getdh(p).grid, "load"))[1]
    return Dict{Int, Vector{T}}(fnode => f)
end

"""
```
                                                               1

                                                              OOO
                                                              ...
                                                              . .
                                                           4  . .
                                30                            . .
/ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . <-
/ .                                                                 . <- 2 f
/ .    3                                                            . <-
/ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . <-
                                                              ^^^
                                                              |||
                                                              1 f
```
"""
struct TieBeam{T, N, M, TBool, TInds<:AbstractVector{Int}, TMeta<:Metadata, CH<:ConstraintHandler{<:DofHandler{2, N, T, M}, T}} <: StiffnessTopOptProblem{2, T}
    E::T
    ν::T
    force::T
    ch::CH
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end
function TieBeam(::Type{Val{CellType}}, ::Type{T}=Float64, refine = 1, force=T(1); E=T(1), ν=T(0)) where {T, CellType}
    grid = TieBeamGrid(Val{CellType}, T, refine)
    dh = DofHandler(grid)
    if CellType === :Linear
        push!(dh, :u, 2)
    else
        ip = Lagrange{2, RefCube, 2}()
        push!(dh, :u, 2, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getfaceset(grid, "leftfixed"), (x,t)->T[0, 0], [1, 2])
    add!(ch, dbc1)
    dbc2 = Dirichlet(:u, getfaceset(grid, "toproller"), (x,t)->T[0], [2])
    add!(ch, dbc2)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return TieBeam(E, ν, force, ch, black, white, varind, metadata)
end

getdim(::TieBeam) = 2
nnodespercell(::TieBeam{T, N}) where {T, N} = N
getpressuredict(p::TieBeam{T}) where {T} = Dict{String, T}("rightload"=>2*p.force, "bottomload"=>-p.force)
getfacesets(p::TieBeam) = getdh(p).grid.facesets


text:
### Processing GEOJSON file: INPUT is path_to_json_data and OUTPUT is distance_matrix_km, number_zones

function processgeodata(path_to_json_data, path_to_data, csv_dataset_list, path_to_results)

    open(path_to_json_data, "r") do json_file
        global julia_file
        julia_file = JSON.parse(json_file);
    end

    a = julia_file["features"]
    first_zone_id = parse(Int,a[1]["properties"]["MOVEMENT_ID"]);
    if (first_zone_id == 0) # catches the case that the first zone ID is a zero
        number_zones = parse(Int,a[length(a)]["properties"]["MOVEMENT_ID"]) + 1;
    else
        number_zones =  parse(Int,a[length(a)]["properties"]["MOVEMENT_ID"]);
    end
    longitude_matrix = zeros(number_zones, 100000);
    latitude_matrix = zeros(number_zones, 100000);

    for i=1:length(a)
        entry = a[i];
        id = parse(Int,entry["properties"]["MOVEMENT_ID"]);
        if (id == 0) # catches the case that the first zone ID is a zero
            id = number_zones;
        end
        coordinates = entry["geometry"]["coordinates"];
        counter = 0;

        if(typeof(coordinates) == Float64)
            counter = counter + 1;
            longitude_matrix[id, counter] = coordinates[1];
            latitude_matrix[id, counter] = coordinates[2];
        else
            for j=1:length(coordinates)
                nested_coordinates = coordinates[j]
                if(typeof(nested_coordinates) == Float64)
                    counter = counter + 1;
                    longitude_matrix[id, counter] = coordinates[1];
                    latitude_matrix[id, counter] = coordinates[2];
                else
                    for k=1:length(nested_coordinates)
                        nested_coordinates = coordinates[j][k];
                        if(typeof(nested_coordinates) == Float64)
                            counter = counter + 1;
                            longitude_matrix[id, counter] = coordinates[j][1];
                            latitude_matrix[id, counter] = coordinates[j][2];
                        else
                            for l=1:length(nested_coordinates)
                                nested_coordinates = coordinates[j][k][l];
                                if(typeof(nested_coordinates) == Float64)
                                    counter = counter + 1;
                                    longitude_matrix[id, counter] = coordinates[j][k][1];
                                    latitude_matrix[id, counter] = coordinates[j][k][2];
                                else
                                    for m=1:length(nested_coordinates)
                                        nested_coordinates = coordinates[j][k][l][m];
                                        if(typeof(nested_coordinates) == Float64)
                                            counter = counter + 1;
                                            longitude_matrix[id, counter] = coordinates[j][k][l][1];
                                            latitude_matrix[id, counter] = coordinates[j][k][l][2];
                                        else
                                            for n=1:length(nested_coordinates)
                                                nested_coordinates = coordinates[j][k][l][m][n];
                                                if(typeof(nested_coordinates) == Float64)
                                                    counter = counter + 1;
                                                    longitude_matrix[id, counter] = coordinates[j][k][l][m][1];
                                                    latitude_matrix[id, counter] = coordinates[j][k][l][m][2];
                                                else
                                                    for o=1:length(nested_coordinates)
                                                        nested_coordinates = coordinates[j][k][l][m][n][o];
                                                        if(typeof(coordinates) == Float64)
                                                            counter = counter + 1;
                                                            longitude_matrix[id, counter] = coordinates[j][k][l][m][n][1];
                                                            latitude_matrix[id, counter] = coordinates[j][k][l][m][n][2];
                                                        else
                                                            println("Double-check how often this json file is nested at:")
                                                            println("i = ", i)
                                                            println("j = ", j)
                                                            println("k = ", k)
                                                            println("l = ", l)
                                                            println("m = ", m)
                                                            println("n = ", n)
                                                            println("o = ", o)
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    centroid_lat = zeros(number_zones,1);
    centroid_long = zeros(number_zones,1);
    area_polygon = zeros(number_zones,1);

    ### Calculating the area of each polygon
    for i=1:number_zones
        if(longitude_matrix[i,1] != 0)
            summation_term = 0;
            for j=1:10000
                if(longitude_matrix[i,j] == 0)
                    ### Correction of summation term
                    summation_term = summation_term - (latitude_matrix[i,j-1] * longitude_matrix[i,j] - latitude_matrix[i,j] * longitude_matrix[i,j-1])

                    ### Addition of first vertex as last one
                    longitude_matrix[i,j] = longitude_matrix[i,1];
                    latitude_matrix[i,j] = latitude_matrix[i,1];

                    ### Recalculation of last summation term
                    summation_term = summation_term + (latitude_matrix[i,j-1] * longitude_matrix[i,j] - latitude_matrix[i,j] * longitude_matrix[i,j-1])
                    break
                else
                    summation_term = summation_term + (latitude_matrix[i,j] * longitude_matrix[i,j+1] - latitude_matrix[i,j+1] * longitude_matrix[i,j])
                end
            end
            area_polygon[i] = summation_term/2;
        end
    end

    ### Calculating the polygon centroid coordinates
    for i=1:number_zones
        if(longitude_matrix[i,1] != 0)
            summation_long = 0;
            summation_lat = 0;
            for j=1:10000
                if(longitude_matrix[i,j] == 0)
                    ### Correction of summation term
                    summation_long = summation_long - (longitude_matrix[i,j-1] + longitude_matrix[i,j]) * (longitude_matrix[i,j-1] * latitude_matrix[i,j] - longitude_matrix[i,j] * latitude_matrix[i,j-1]);
                    summation_lat = summation_lat - (latitude_matrix[i,j-1] + latitude_matrix[i,j]) * (longitude_matrix[i,j-1] * latitude_matrix[i,j] - longitude_matrix[i,j] * latitude_matrix[i,j-1]);
                else
                    summation_long = summation_long + (longitude_matrix[i,j] + longitude_matrix[i,j+1]) * (longitude_matrix[i,j] * latitude_matrix[i,j+1] - longitude_matrix[i,j+1] * latitude_matrix[i,j]);
                    summation_lat = summation_lat + (latitude_matrix[i,j] + latitude_matrix[i,j+1]) * (longitude_matrix[i,j] * latitude_matrix[i,j+1] - longitude_matrix[i,j+1] * latitude_matrix[i,j]);
                end
            end
            centroid_long[i] = -summation_long/(6 * area_polygon[i]);
            centroid_lat[i] = -summation_lat/(6 * area_polygon[i]);
        end
    end

    ## calculating geographic distances between zones
    c_lat_long = 111.3; # this is a natural constant
    conv_deg_rad = 0.01745; # this is a constant for the conversion from degrees into radians. It serves for input into trig. functions.
    distance_matrix_km = zeros(number_zones, number_zones);
    for i=1:(number_zones-1)
        long1 = centroid_long[i];
        lat1 = centroid_lat[i];
        for j=i:number_zones
            long2 = centroid_long[j];
            lat2 = centroid_lat[j];
            distance_km = c_lat_long * sqrt((cos((lat1+lat2)/2 * conv_deg_rad))^2 * (long1-long2)^2 + (lat1-lat2)^2);
            distance_matrix_km[i,j] = distance_km;
            distance_matrix_km[j,i] = distance_km;
        end
    end

    ## Filling the diagonal zeros that stand for a travel distance within the same zone with a distance of 1 km.
    for i=1:number_zones
        distance_matrix_km[i,i] = 1;
    end

    ## save the results
    results = zeros(number_zones, 2)
    for i = 1:number_zones
        results[i,1] = centroid_lat[i];
        results[i,2] = centroid_long[i];
    end
    results = convert(DataFrame, results);
    header_vector = ["latitude", "longitude"]; # creates an empty string vector in the desired dimension
    results_string = string(path_to_results,"/zoneID_coordinates.csv");
    CSV.write(results_string, results, header=header_vector);

    return distance_matrix_km, number_zones

end

ext: jl

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
export LPHRep

MOI.Utilities.@model(_MOIModel,
                     (), (MOI.EqualTo, MOI.LessThan,), (), (),
                     (), (MOI.ScalarAffineFunction,), (), ())
# We need the `SingleVariable` constraints to be bridged so we should say that
# they are not supported. We notably exclude `Integer` as we just ignore
# integrality constraints. Binary constraint should be bridged to integrality
# once https://github.com/JuliaOpt/MathOptInterface.jl/issues/704 is done.
function MOI.supports_constraint(
    ::_MOIModel{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T},
                   MOI.Interval{T}, MOI.ZeroOne}}) where T
    return false
end

mutable struct LPHRep{T} <: HRepresentation{T}
    model::_MOIModel{T}
    hyperplane_indices::Union{Nothing,
                              Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                                         MOI.EqualTo{T}}}}
    halfspace_indices::Union{Nothing,
                             Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                                        MOI.LessThan{T}}}}
end
function LPHRep(model::_MOIModel{T}) where T
    return LPHRep(model, nothing, nothing)
end
function LPHRep(model::MOI.ModelLike, T::Type = Float64)
    _model = _MOIModel{T}()
    bridged = MOI.Bridges.LazyBridgeOptimizer(_model)
    # Only enable constraint bridges that don't create variables and don't add
    # any variable bridge so that there is an identity mapping betwenen
    # variables of `model` and polyhedra dimensions.
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.GreaterToLessBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.LessToGreaterBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NonnegToNonposBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NonposToNonnegBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.VectorizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarFunctionizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.VectorFunctionizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.SplitIntervalBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NormInfinityBridge{T})
    # This one creates variables so the user need to consider the polyhedra as the
    # feasible set of the extended formulation.
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NormOneBridge{T})
    MOI.copy_to(bridged, model)
    return LPHRep(_model)
end
# returns `Int64` so need to convert for 32-bit system
FullDim(rep::LPHRep) = convert(Int, MOI.get(rep.model, MOI.NumberOfVariables()))

hvectortype(::Type{LPHRep{T}}) where {T} = SparseVector{T, Int}
similar_type(::Type{LPHRep{S}}, ::FullDim, ::Type{T}) where {S, T} = LPHRep{T}
fulltype(::Type{LPHRep{T}}) where {T} = LPHRep{T}

LPHRep(h::HRep{T}) where {T} = LPHRep{T}(h)
LPHRep{T}(h::HRep) where {T} = convert(LPHRep{T}, h)

function _saf(a::AbstractVector{T}, vars) where T
    func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], zero(T))
    for i in eachindex(a)
        coef = a[i]
        if !iszero(coef)
            push!(func.terms, MOI.ScalarAffineTerm(coef, vars[i]))
        end
    end
    return func
end

function LPHRep{T}(d::FullDim,
                   hyperplanes::ElemIt{<:HyperPlane{T}},
                   halfspaces::ElemIt{<:HalfSpace{T}}) where {T}
    model = _MOIModel{T}()
    vars = MOI.add_variables(model, fulldim(d))
    for hyperplane in hyperplanes
        func = _saf(hyperplane.a, vars)
        set = MOI.EqualTo(hyperplane.β)
        MOI.add_constraint(model, func, set)
    end
    for halfspace in halfspaces
        func = _saf(halfspace.a, vars)
        set = MOI.LessThan(halfspace.β)
        MOI.add_constraint(model, func, set)
    end
    return LPHRep(model)
end

function Base.copy(lp::LPHRep{T}) where {T}
    model = _MOIModel{T}()
    MOI.copy_to(model, lp.model)
    return LPHRep(model)
end

function constraint_indices(rep::LPHRep{T},
                            ::Union{HyperPlaneIndex, HyperPlaneIndices}) where T
    if rep.hyperplane_indices === nothing
        rep.hyperplane_indices = MOI.get(
            rep.model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},
                                        MOI.EqualTo{T}}())
    end
    return rep.hyperplane_indices
end
function MOI.add_constraint(rep::LPHRep{T},
                            func::MOI.ScalarAffineFunction{T},
                            set::MOI.EqualTo{T}) where T
    rep.hyperplane_indices = nothing
    return MOI.add_constraint(rep.model, func, set)
end

function constraint_indices(rep::LPHRep{T},
                            ::Union{HalfSpaceIndex, HalfSpaceIndices}) where T
    if rep.halfspace_indices === nothing
        rep.halfspace_indices = MOI.get(
            rep.model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},
                                        MOI.LessThan{T}}())
    end
    return rep.halfspace_indices
end
function MOI.add_constraint(rep::LPHRep{T},
                            func::MOI.ScalarAffineFunction{T},
                            set::MOI.LessThan{T}) where T
    rep.halfspace_indices = nothing
    return MOI.add_constraint(rep.model, func, set)
end

function Base.length(idxs::HIndices{T, LPHRep{T}}) where {T}
    return length(constraint_indices(idxs.rep, idxs))
end
function Base.isempty(idxs::HIndices{T, LPHRep{T}}) where {T}
    return isempty(constraint_indices(idxs.rep, idxs))
end
function startindex(idxs::HIndices{T, LPHRep{T}}) where {T}
    if isempty(idxs)
        return nothing
    else
        return eltype(idxs)(1)
    end
end
function Base.get(rep::LPHRep{T}, idx::HIndex{T}) where {T}
    ci = constraint_indices(rep, idx)[idx.value]
    func = MOI.get(rep.model, MOI.ConstraintFunction(), ci)::MOI.ScalarAffineFunction
    indices = [t.variable_index.value for t in func.terms]
    values = [t.coefficient for t in func.terms]
    a = sparsevec(indices, values, FullDim(rep))
    set = MOI.get(rep.model, MOI.ConstraintSet(), ci)
    β = MOI.constant(set) - func.constant
    if idx isa HyperPlaneIndex
        return HyperPlane(a, β)
    else
        return HalfSpace(a, β)
    end
end
function nextindex(rep::LPHRep{T}, idx::HIndex{T}) where {T}
    if idx.value >= length(constraint_indices(rep, idx))
        return nothing
    else
        return typeof(idx)(idx.value + 1)
    end
end

function Base.isvalid(lp::LPHRep{T}, idx::HIndex{T}) where {T}
    return 1 ≤ idx.value ≤ length(constraint_indices(rep, idx))
end

dualtype(::Type{LPHRep{T}}, ::Type{AT}) where {T, AT} = dualtype(Intersection{T, AT, Int}, AT)
dualfullspace(h::LPHRep, d::FullDim, ::Type{T}, ::Type{AT}) where {T, AT} = dualfullspace(Intersection{T, AT, Int}, d, T, AT)


text:
using IntervalWavelets
using Wavelets
using Base.Test

println("Testing dyadic rationals...")

#=
Translation between indices and values for dyadic rationals
=#

p = 2
S = DaubSupport(-p+1, p)

# Level 0 = integers
supportx = dyadic_rationals( S, 0 )
for x in round(Int, supportx)
	@test supportx[ x2index(x,S) ] == x
end

for idx in 1:length(supportx)
	@test supportx[idx] == index2x(idx,S)
end

# Higher level
R = rand( 1:10 )

supportx = dyadic_rationals( S, R )
for x in supportx
	@test supportx[ x2index(x,S,R) ] == x
end

for idx in 1:length(supportx)
	@test supportx[idx] == index2x(idx,S,R)
end


text:
# All the types native to ogl, wgl and vulkan shaders

const number_types = (Float32, Cint, Cuint, Cdouble)
const small_vecs = (((StaticVector{N, T} for T in number_types, N in (2, 3, 4)))...,)
const small_mats = (((StaticArray{Tuple{i, j}, T, 2} for T in ShaderAbstractions.number_types, i in (2, 3, 4), j in (2, 3, 4)))...,)
const small_arrays = (small_vecs..., small_mats...)
const native_types = (number_types..., small_arrays...)

const NativeNumbers = Union{number_types...}
const SmallVecs = Union{small_vecs...}
const SmallMats = Union{small_mats...}
const SmallArrays = Union{small_arrays...}
const NativeTypes = Union{native_types...}

"""
    native_type(context, Int128)
Returns a native type for a non native type.
E.g. native_type(context, Int128) -> Cint
"""
native_type(context::AbstractContext, ::Type{T}) where T <: NativeTypes = T

native_type(context::AbstractContext, x::Type{T}) where {T <: Integer} = Cint
native_type(context::AbstractContext, x::Type{Union{Int16, Int8}}) = x

native_type(context::AbstractContext, x::Type{T}) where {T <: Unsigned} = Cuint
native_type(context::AbstractContext, x::Type{Union{UInt16, UInt8}})  = x

native_type(context::AbstractContext, x::Type{T}) where {T <: AbstractFloat} = Float32
native_type(context::AbstractContext, x::Type{Float16})               = x

native_type(context::AbstractContext, x::Type{T}) where T <: Normed = N0f32
native_type(context::AbstractContext, x::Type{N0f16}) = x
native_type(context::AbstractContext, x::Type{N0f8}) = x

native_type(context::AbstractContext, x::Type{<: StaticArray{S, T, N}}) where {S, T, N} = similar_type(x, native_type(context, T))

map_t(f, tuple) = map_t(f, (), tuple)
map_t(f, result, ::Type{Tuple{}}) = Tuple{result...}
function map_t(f, result, T::Type{<: Tuple})
    map_t(
        f,
        (result..., f(Base.tuple_type_head(T))),
        Base.tuple_type_tail(T)
    )
end

native_type(context::AbstractContext) = x-> native_type(context, x)

function native_type(context::AbstractContext, ::Type{NamedTuple{Names, Types}}) where {Names, Types}
    return NamedTuple{Names, map_t(native_type(context), Types)}
end

"""
Curry form to easier pass this via e.g. map
"""
convert_uniform(context::AbstractContext) = x-> convert_uniform(context, x)

"""
All native types don't need conversion
"""
convert_uniform(context::AbstractContext, x::NativeTypes) = x

"""
Vector of native types, e.g `vec3 [4]`
"""
convert_uniform(context::AbstractContext, x::StaticVector{N, T}) where {N, T <: NativeTypes} = x

"""
Static Array with non native uniform type
"""
function convert_uniform(context::AbstractContext, x::StaticVector{N, T}) where {N, T}
    convert(similar_type(x, native_type(context, T)), x)
end

"""
Colors get special treatment to get them as vecs in shaders
"""
convert_uniform(context::AbstractContext, x::Colorant{T}) where T <: NativeNumbers = x

convert_uniform(context::AbstractContext, x::Colorant{T}) where T = mapc(native_type(context, T), x)

convert_uniform(context::AbstractContext, x::AbstractSampler{T, N}) where {T, N} = x

function convert_uniform(context::AbstractContext, x::AbstractVector{T}) where T
    return convert(Vector{native_type(context, T)}, x)
end

function convert_uniform(context::AbstractContext, x::NamedTuple{Names, Types}) where {Names, Types}
    return map(convert_uniform(context), x)
end

function convert_uniform(context::AbstractContext, x::Observable)
    return map(convert_uniform(context), x)
end
function convert_uniform(context::AbstractContext, x::T) where T
    all(t-> isbits(t), fieldtypes(T)) || error("All field types need to be isbits. Found: $(T) with $(fieldtypes(T))")
    return x
end

type_prefix(context::AbstractContext, x::Type{T}) where {T <: Union{FixedPoint, Float32, Float16}} = ""
type_prefix(context::AbstractContext, x::Type{T}) where {T <: Float64} = "d"
type_prefix(context::AbstractContext, x::Type{Cint}) = "i"
type_prefix(context::AbstractContext, x::Type{T}) where {T <: Union{Cuint, UInt8, UInt16}} = "u"

type_postfix(context::AbstractContext, x::Type{Float64}) = "dv"
type_postfix(context::AbstractContext, x::Type{Float32}) = "fv"
type_postfix(context::AbstractContext, x::Type{Cint})    = "iv"
type_postfix(context::AbstractContext, x::Type{Cuint})   = "uiv"

type_string(context::AbstractContext, x::T) where {T} = type_string(context, T)
type_string(context::AbstractContext, t::Type{Float32}) = "float"
type_string(context::AbstractContext, t::Type{Float64}) = "double"
type_string(context::AbstractContext, t::Type{Cuint}) = "uint"
type_string(context::AbstractContext, t::Type{Cint}) = "int"

function type_string(context::AbstractContext, t::Type{T}) where {T <: Union{StaticVector, Colorant}}
    return string(type_prefix(context, eltype(T)), "vec", length(T))
end

function type_string(context::AbstractContext, t::Type{<: AbstractSamplerBuffer{T}}) where T
    string(type_prefix(context, eltype(T)), "samplerBuffer")
end
is_arraysampler(t) = false
function type_string(context::AbstractContext, t::Type{<: AbstractSampler{T, D}}) where {T, D}
    str = string(type_prefix(context, eltype(T)), "sampler", D == 1 ? 2 : D, "D")
    is_arraysampler(t) && (str *= "Array")
    return str
end

function type_string(context::AbstractContext, t::Type{<: StaticMatrix})
    M, N = size(t)
    string(type_prefix(context, eltype(t)), "mat", M == N ? M : string(M, "x", N))
end
type_string(context::AbstractContext, t::Bool) = "bool"

type_string(context::AbstractContext, t::Observable) = type_string(context, t[])

type_string(context::AbstractContext, t::Type) = error("Type $t not supported")


text:
module PhyDCA

using GaussDCA
using PlmDCA
using PyPlot
using Printf

export
	#phylogenetic distance class
	Hamming,
	Correlation,
	pValue,
	mfDCA,
	plmDCA,
	PhyloOut,

	#function
    phydca,
	make_phylo_profile,
	evaluate_distance,
	compute_distance_matrix

include("io.jl")
include("types.jl")
include("fisher.jl")
include("commonDCA.jl") #tools for DCA inference, both for mfDCA and plmDCA

########################################
# MAIN FUNCTIONS
########################################
#method 1,  directly the phylogenetic matrix (MxN) file containing knwon functional relations is given
function phylodca_matrix(filename_phyloMatrix::AbstractString,
                         dist::UnionDistances)

    P = readdlm(filename_phyloMatrix)
    P = Array{Int8}(P)
    M,N = size(P)

    sss = collect(1:M)
    species = [string(sss[i]) for i =1:M]
    ddd= collect(1:N)
    domains= [string(ddd[i]) for i =1:N]

	@printf("%s", "Computing distance matrix...")
	dist_matr=compute_distance_matrix(P,dist)
	@printf("%s\n", "done")

	@printf("%s\n", "Sorting domain-pairs by their phylogenetic distance")
    final_unsort, final_sorted=compute_final_matrix(dist_matr,domains, dist)

	#print_result(final_sorted)

    x=PhyloOut(P,domains,species,dist_matr,final_unsort, final_sorted)
	return x
end

#method 2, no file containing known functional relations is given
function phylodca(filename_phylo::AbstractString,
                  dist::UnionDistances)

	@printf("%s", "Computing phylogenetic matrix...")
	P,species,domains=make_phylo_profile(filename_phylo)
	@printf("%s\n", "done")

	@printf("%s", "Computing distance matrix...")
	dist_matr=compute_distance_matrix(P,dist)
	@printf("%s\n", "done")

	@printf("%s\n", "Sorting domain-pairs by their phylogenetic distance")
    final_unsort, final_sorted=compute_final_matrix(dist_matr,domains, dist)

	#print_result(final_sorted)

    x=PhyloOut(P,domains,species,dist_matr,final_unsort, final_sorted)
	return x
end

#method 3: one (or more) file(s) containing known functional relation is given
function phylodca(filename_phylo::AbstractString,
		  dist::UnionDistances,
          known_relations::Array{String,1}) #VarArg allows variable number of inputs

	@printf("%s", "Computing phylogenetic matrix...")
	P,species,domains = make_phylo_profile(filename_phylo)
	@printf("%s\n", "done")

	@printf("%s\n", "Computing distance matrix...")
	dist_matr = compute_distance_matrix(P,dist)

	@printf("%s\n", "Sorting domain-pairs by their phylogenetic couplings")
    final_unsort, final_sorted = compute_final_matrix(dist_matr,domains, dist,known_relations)

	#print_result(final_sorted)

    x = PhyloOut(P,domains,species,dist_matr,final_unsort, final_sorted)
	return x
end

##############################
#1) construct phylogenetic matrix
##############################
function make_phylo_profile(filename::AbstractString)

	list_species,list_domains = readData(filename)
	P = makePhyloMatrix(filename, list_species,list_domains)

	#PLOT phylogenetic matrix.... not a good idea...

	#plot_matrix(P,list_domains,list_species)
	#plot_hist(P)
	#implementa plot_distr_domains (per vedere se mettere cut_off)

	return P,list_species,list_domains

end

##############################
#2) construct distance matrix
##############################
function compute_distance_matrix(P::Matrix{Int8},
				 dist::UnionDistances)

	!isa(P,Matrix) && error("Phlogenetic Matrix not given")

	dist_matrix = evaluate_distance(dist,P)
end

############################################################
#3) construct result matrix sorted by phylogenetic distances
############################################################
#method 1, no file containing knwown functional relations is given
function compute_final_matrix(distance_matrix::Matrix{Float64},
			      list_domains::Array{String},
			      dist::UnionDistances)

	res_matr=make_result_matrix(distance_matrix,list_domains)
	sort_res_matr=sort_matrix(res_matr,dist)
	return res_matr,sort_res_matr
end

#method 2: one (or more) file(s) containing kwnown functional relations is given
function compute_final_matrix(distance_matrix::Matrix{Float64},
			      list_domains::Array{String},
			      dist::UnionDistances,
                  bench_filenames::Array{String,1}) #Vararg allows variable number of inputs

	res_matr=make_result_matrix(distance_matrix,list_domains)

	for i=1:length(bench_filenames)
		println(bench_filenames[i])
		res_matr=add_benchmark(res_matr,list_domains,bench_filenames[i])
	end

	sort_res_matr=sort_matrix(res_matr,dist)
	return res_matr,sort_res_matr
end

######################################################################
# plot the PPV (positive predictive value)
######################################################################
function  plot_PPV(final_matrix::Matrix{Any},
		   ppv_index::Vararg{Int})
	ss=sum(final_matrix[:,k] for k in ppv_index)
	num_pred,ppv_rate=compute_tp_rate(Array{Int}(ss))
	xscale("log")
	xlabel("Number of predictions")
	ylabel("PPV")
	plot(num_pred,ppv_rate)
end

####################################################################################################
#plotting the matrix
####################################################################################################
function plot_matrix(P::Matrix{Int8},
		     x_label,
		     y_label)
    fig = figure("Phylogenetic Profile", figsize=(10,10))
    ax = gca()
    x_label=union(0,x_label)
    y_label=union(0,y_label)
    ax[:set_xticklabels](x_label)
    ax[:set_yticklabels](y_label)
    ax[:matshow](P) #cmap=plt.cm.gray
end

function plot_hist(P::Matrix{Int8})
    fig = figure("Perc of batt", figsize=(10,10))
    ax = gca()
    s=sum(P,1)
    ax[:hist](s) #cmap=plt.cm.gray

end

####################################################################################################

###############
#compare with set of known relations
##############

#for j>i, the pair {i,j} in the list {{1,2}, {1,3}, {1,4},...,{1,N},{2,3},... } is located in offset(i,j,N)
function offset(i::Int,j::Int,N::Int)
	off=(i-1)*N - i*(i-1)/2 +(j-i)
	return Int(off)
end

function make_result_matrix(Jij::Matrix{Float64},
			    list_domains::Array{String})
	N=length(list_domains)
	num=Int(N*(N-1)/2)
	result_matrix=Matrix{Any}(num,3)

	cont=1
	for i=1:N
		for j=(i+1):N
			result_matrix[cont,1]=list_domains[i]
			result_matrix[cont,2]=list_domains[j]
			result_matrix[cont,3]=Jij[i,j]
			cont+=1
		end
	end
	return result_matrix
end

function add_benchmark(result_matrix::Matrix,
		       list_domains::Array{String},
		       filename::AbstractString)

		bench=read_benchmark(filename)
		tot=size(result_matrix)[1]
		vec_bench=zeros(Int,tot)
		Nb=size(bench)[1]

        if(size(bench)[2]!=2)
            error("knwon relations must be pairwise")
        end

		num_domains=size(list_domains)[1]
        dict_dom = Dict(list_domains[i] => i for i=1:length(list_domains))

		for i=1:Nb
			dom1=bench[i,1]
			dom2=bench[i,2]
            if ( haskey(dict_dom, dom1) && haskey(dict_dom, dom2))
                pos1 = dict_dom[dom1]
                pos2 = dict_dom[dom2]
                if(pos1<pos2)
                    vec_bench[offset(pos1,pos2,num_domains)]=1
                else
                    vec_bench[offset(pos2,pos1,num_domains)]=1
                end
            end
		end

		result_matrix=hcat(result_matrix,vec_bench)

	return result_matrix
end

#N.b sort matrix must depend on the phylogenetic distance:
#Hamming, pValue distance-> minimum
#correlation, mfDCA,plmDCA-> maximum
#
function sort_matrix(final_matrix::Matrix,
		     ::Union{Hamming,pValue})
	sorted_matr=sortrows(final_matrix, by=x->x[3])
	return sorted_matr
end
function sort_matrix(final_matrix::Matrix,
		     ::Union{Correlation,mfDCA,plmDCA,MutualInfo})
	sorted_matr=sortrows(final_matrix, by=x->x[3],rev=true)
	return sorted_matr
end

#################################
##plot tp rate
###################

function compute_tp_rate(tp::Array{Int})
    m=length(tp)
    tp_check=zeros(m)
    for a=1:m
        tp_check[a]=(tp[a]!=0)
    end
    tp_cumsum=cumsum(tp_check)

    num_pred=zeros(m)
    ppv_rate=zeros(m)
    for a=1:m
	    num_pred[a]=a
	    ppv_rate[a]=tp_cumsum[a]/a
    end
    return num_pred,ppv_rate
end

end

max_stars_repo_name: clayfos/PhyDCA

max_issues_repo_name: clayfos/PhyDCA

max_forks_repo_name: clayfos/PhyDCA


text:
include("structure.jl")
include("initialize.jl")
include("shuffle_config.jl")
include("sro_calc.jl")
include("plotting.jl")

"""
WELCOME TO MY SRO CALCULATOR.
This script calculates Warren-Cowley short range order parameters for pairs
and triplets in arbitrary alloys. The code works for an arbitrary number of
different elements. This file shall guide you through the necessary steps.
"""

# INPUT
# Please declare your input YAML-file:

file = "testfile.yaml"

# This will create all structures for a SRO calculation and stored it in an
# object of type 'Backpack'.

my_backpack = initialize(file)

# Now, the SRO parameters are being calculated while the atomic configuration is
# shuffled in between. The input number below determines the number of
# non-repetitive shuffles:

shuffles = 100
perms = shuffle_atoms(my_backpack, shuffles)
#xyz_animation(my_backpack, unique_perms[2])

# Alternatively, if you want to generate aLL permutations, you can use this function
# below. You find the documentations for this functions in 'shuffle_sonfig.jl' as well.
# Please keep in mind that this can be very CPU-intensive and is therefore not
# suitable for too large systems. But feel free to try it out. ;-)
# @time begin
    # perms = all_permutations(sequence(my_backpack))
# end

# Now, we start SRO calculation. With the function 'evaluate', SROs are calculated for
# all created permutations and the minima and maxima for pairs, triplets and a
# sum of both are returned separately. The value in the function is the triplet weight factor,
# which defines the contribution of the triplets
results = evaluate(my_backpack, perms, 0.2)

# Finally, see what you did! :-)
# You can either plot some reults, like here the pairs...
plot_configs(results.pairs)
#...or you plot the whole structure, for example the minimum for triplets...
#xyzplot(my_backpack, perms[results.triplet_extrema[1][2]], 15.0)
#...you can even create a little gif, if you want:
#xyz_animation(my_backpack, perms[results.triplet_extrema[1][2]])

max_stars_repo_name: ulipoto/WCSRO

max_issues_repo_name: ulipoto/WCSRO

max_forks_repo_name: ulipoto/WCSRO


text:
@testset "Knapsack" begin

@testset "Knapsack: ZeroOnePack!" begin
    dp=zeros(Int,30)
    @test ZeroOnePack!(20,[1,3,11],[2,5,30],dp) == 37
    dp=zeros(Int,30)
    @test ZeroOnePack!(10,[1,3,11],[20,5,80],dp) == 25
end

@testset "Knapsack: CompletePack!" begin
    dp=zeros(Int,30)
    @test CompletePack!(20,[1,2,9],[1,3,20],dp) == 43
    dp=zeros(Int,30)
    @test CompletePack!(10,[1,3,11],[20,5,80],dp) == 200
end

end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
div1(x, y) = div(x - 1, y) + 1

w_eachindex(m::Matrix) = eachindex(m)

struct TbBlock{T <: AbstractFloat, LT, M <: AbstractMatrix{Complex{T}}}
    R_cryst ::Vec3{Int}
    R_cart  ::Vec3{LT}
    block::M #has contributions from all the ws shifts and divided by the degeneracies 1/(degen * nshifts), use for k-point interpolation
    tb_block::M #this is the correct tight binding block, with ws shifts applied such that entries represent hopping between closest two wfs
end

block(x::TbBlock) = x.block

for f in (:getindex, :size, :similar)
	@eval Base.$f(h::TbBlock, args...) = $f(block(h), args...)
end

LinearAlgebra.eigen(h::TbBlock) =
	eigen(block(h))

const TbHami{T, LT, M}  = Vector{TbBlock{T, LT, M}}

getindex(h::TbHami, R::Vec3{Int}) =
	getfirst(x -> x.R_cryst == R, h)

#some small type piracy?
Base.zeros(m::AbstractArray{T}) where {T} =
	fill!(similar(m), zero(T))

zeros_block(h::TbHami) =
	zeros(block(h[1]))

similar_block(h::TbHami) =
	similar(block(h[1]))

blocksize(h::TbHami, args...) =
	size(block(h[1]), args...)

for op in (:+, :-, :*, :/)
	@eval $op(t::TbBlock{T}, v::T) where {T} = TbBlock(t.R_cart, t.R_cryst, $op(block(t), v), $op(t.tb_block, v))
	@eval $op(v::T, t::TbBlock{T}) where {T} = TbBlock(t.R_cart, t.R_cryst, $op(v, block(t)), $op(v, t.tb_block))
	@eval $op(t::TbBlock{T,M}, v::M) where {T,M} = TbBlock(t.R_cart, t.R_cryst, $op(block(t), v), $op(t.tb_block, v))
	@eval $op(v::M, t::TbBlock{T,M}) where {T,M} = TbBlock(t.R_cart, t.R_cryst, $op(v, block(t)), $op(v, t.tb_block))
	@eval $op(t::TbBlock{T,M}, v::TbBlock{T,M}) where {T,M} = TbBlock(t.R_cart, t.R_cryst, $op(block(t), block(v)), $op(t.tb_block, v.tb_block))
end

make_noncolin(tb::TbBlock) =
    TbBlock(tb.R_cryst, tb.R_cart, convert(NonColinMatrix, tb.block), convert(NonColinMatrix, tb.tb_block))

make_noncolin(tb::TbBlock{T, LT, ColinMatrix{Complex{T}, Matrix{Complex{T}}}}) where {T<:AbstractFloat, LT<:Length{T}}=
    TbBlock(tb.R_cryst, tb.R_cart, NonColinMatrix(tb.block[Up()], tb.block[Down()]), NonColinMatrix(tb.tb_block[Up()], tb.tb_block[Down()]))

make_noncolin(v::Vector) =
    [v[1:2:end];v[2:2:end]]

struct HamiltonianKGrid{T,MT<:AbstractMatrix{Complex{T}}, VT<:AbstractVector{T}} <: AbstractKGrid{T}
    core::CoreKGrid{T}
    Hk::Vector{MT}
    eigvals::Vector{VT}
    eigvecs::Vector{MT}
end
HamiltonianKGrid(kpoints::Vector{<:Vec3}, args...) = HamiltonianKGrid(CoreKGrid(kpoints), args...)

@doc raw"""
	HamiltonianKGrid(hami::TbHami{T}, nk, H_function_k::Function = x -> nothing) where T
	HamiltonianKGrid(hami::TbHami{T}, k_grid, H_function_k::Function = x -> nothing) where T

Takes a k grid, calculates Hk for each of them and diagonalizes. Only the eigenvectors and eigenvalues of Hk are stored,
the H_function_k function is called on the intermediate Hk.
"""
function HamiltonianKGrid(hami::TbHami{T}, kpoints::Vector{<:Vec3}, Hk_function::Function = x -> nothing) where {T}
	# kpoints = [KPoint(k, blocksize(hami), R, zeros_block(hami)) for k in k_grid]
	n_eigvals = max(blocksize(hami)...)
	eigvals = hami[1].block isa AbstractMagneticMatrix ? [MagneticVector(zeros(T, n_eigvals)) for k in kpoints] : [zeros(T, n_eigvals) for k in kpoints]
	kgrid = HamiltonianKGrid(kpoints, [zeros_block(hami) for k in kpoints], eigvals, [zeros_block(hami) for k in kpoints])
	nk    = length(kpoints)
	calc_caches = [EigCache(block(hami[1])) for i=1:nthreads()]
	p = Progress(nk, 1, "Calculating H(k)...")
    @threads for i=1:nk
	    tid = threadid()
	    Hk!(kgrid.eigvecs[i], hami, k_cryst(kgrid)[i])
        kgrid.Hk[i] = copy(kgrid.eigvecs[i])
	    Hk_function(kgrid.Hk[i])
	    eigen!(kgrid.eigvals[i], kgrid.eigvecs[i], calc_caches[Threads.threadid()])
	    next!(p)
    end
    return kgrid
end

function Hk!(out::AbstractMatrix, tbhami::TbHami, kpoint::Vec3)
    fill!(out, zero(eltype(out)))
    fourier_transform(tbhami, kpoint) do i, iR, R_cart, b, fac
        @inbounds out[i] += fac * b.block[i]
    end
end

function Hk(tbhami::TbHami, kpoint::Vec3)
    out = similar(tbhami[1].block)
    Hk!(out, tbhami, kpoint)
    return out
end

Hk(g::HamiltonianKGrid) = g.Hk
eigvecs(g::HamiltonianKGrid) = g.eigvecs
eigvals(g::HamiltonianKGrid) = g.eigvals

"Fourier transforms the tight binding hamiltonian and calls the R_function with the current index and the phase."
function fourier_transform(R_function::Function, tb_hami::TbHami{T}, kpoint::Vec3) where {T}
    for (iR, b) in enumerate(tb_hami)
        fac = ℯ^(2im*π*(b.R_cryst ⋅ kpoint))
        for i in eachindex(block(b))
            R_function(i, iR, b.R_cart, b, fac)
        end
    end
end

# "Constructs the total spin-orbit-coupled Hamiltonian out of supplied angular momentums between the Wannier functions and uses the l_soc of the atoms."
# function construct_soc_hami(hami, structure::WanStructure{T})::Matrix{Complex{T}} where T
#     dim = getwandim(structure)
#     Lx_soc = zeros(Complex{T}, dim, dim)
#     Ly_soc = zeros(Complex{T}, dim, dim)
#     Lz_soc = zeros(Complex{T}, dim, dim)
#     i = 1
#     for at in atoms(structure)
#         soc = lsoc(at)
#         mom = angmom(at)
#         len = length(wfcs(at))-1
#         Lx_soc[i:i+len, i:i+len] = 0.5 * soc * getindex.(mom, 1)
#         Ly_soc[i:i+len, i:i+len] = 0.5 * soc * getindex.(mom, 2)
#         Lz_soc[i:i+len, i:i+len] = 0.5 * soc * getindex.(mom, 3)
#         i += len + 1
#     end
#     Lx_soc = (Lx_soc+Lx_soc')/2
#     Ly_soc = (Ly_soc+Ly_soc')/2
#     Lz_soc = (Lz_soc+Lz_soc')/2
#     out = [hami+Lz_soc Lx_soc-1im*Ly_soc;Lx_soc+1im*Ly_soc hami-Lz_soc]
#     return out
# end

"Holds all the calculated values from a wannier model."
@with_kw mutable struct WannierBand{T<:AbstractFloat, VT<:AbstractVector} <: DFC.AbstractBand
    kpoints_cryst ::Vector{Vec3{T}}
    eigvals        ::Vector{T}
    eigvec         ::Vector{VT}
    cms      ::Vector{Point3{T}} = Point3{T}[]
    angmoms  ::Vector{Vector{Point3{T}}} = Vector{Point3{T}}[]
    spins    ::Vector{Vector{Point3{T}}} = Vector{Point3{T}}[]
end

DFControl.eigvals(b::WannierBand) = b.eigvals

function wannierbands(tbhamis::TbHami{T}, kpoints::Vector{<:Vec3}) where {T}
    matdim   = blocksize(tbhamis, 2)
    kgrid    = HamiltonianKGrid(tbhamis, kpoints)
    nbnd = size(tbhamis[1].block, 2)
    evals = [zeros(length(kpoints)) for i=1:nbnd]
    evecs = [[similar(kgrid.eigvecs[1][:, 1]) for i=1:length(kpoints)] for i=1:nbnd]
    @threads for i = 1:length(kpoints)
        eigvals, eigvecs = kgrid.eigvals[i], kgrid.eigvecs[i]
        for e=1:length(eigvals)
            evals[e][i] = eigvals[e]
            evecs[e][i] = eigvecs[:,e]
        end
    end
    return [WannierBand{T, eltype(evecs[1])}(kpoints_cryst = kpoints, eigvals=evals[i], eigvec=evecs[i]) for i =1:nbnd]
end
wannierbands(tbhamis, dfbands::Vector{<:DFC.AbstractBand}) =
	wannierbands(tbhamis, dfbands[1].k_points_cryst)
wannierbands(tbhamis, dfbands::Union{NamedTuple, Tuple}) =
	wannierbands(tbhamis, dfbands[1][1].k_points_cryst)

function energy_bins(binfunc::Function, wbands::Vector{<:WannierBand}, E_range, normalize_bins=false)
    nbins   = length(E_range)-1
    bins    = zeros(typeof(binfunc(wbands[1].eigvec[1])), nbins)
    # sum = zero(typeof(binfunc(bands[1].eigvec[1])))
    nperbin = zeros(Int, nbins)
    for b in wbands
        for (e, v) in zip(b.eigvals, b.eigvec)
            ie = findfirst(i -> E_range[i] <= e <= E_range[i+1], 1:nbins)
            if ie === nothing
                continue
            end
            bins[ie] += binfunc(v)
            nperbin[ie] += 1
        end
    end
    if normalize_bins #like taking the mean
        for i=1:nbins
            if nperbin[i] > 0
                bins[i] /= nperbin[i]
            end
        end
    end
    return bins
end

function character_contribution(wband::WannierBand, atoms::Vector{Atom})
    contributions = zeros(length(wband.kpoints_cryst))
    for (i, v) in enumerate(wband.eigvec)
        for a in atoms
            contributions[i] += norm(v[a])^2
        end
    end
    return contributions
end

function DFControl.FileIO.pdos(wbands::Vector{<:WannierBand}, atoms::Vector{Atom}, dE = 0.02)
    Emin = minimum(wbands[1].eigvals)
    Emax = maximum(wbands[end].eigvals)
    E_range = range(Emin, Emax, step=dE)

    bins = energy_bins(wbands, E_range, false) do v
        tot = 0.0
        for a in atoms
            tot += norm(v[a])^2/dE
        end
        return tot
    end
    return (E=E_range, pdos=bins./length(wbands[1].kpoints_cryst))
end

kpdos(bands::Vector{<:WannierBand}, atoms::Vector{Atom}) = map(x -> character_contribution(x, atoms), bands)


text:
"""
extract a fully sampled calibration area of size `M` by `N` from the Matrix `A`

# Arguments
* `A::Array{Int64,2}` - Array to be sampled
* `M::Int64`          - size of A along the first dimension
* `N::Int64`          - size of N along the second dimension
* `redFac::Float64`   - reduction factor of the target sampling pattern (used for a sanity check)
* `calsize::Int64`    - size of the calibration area
"""
function setcalibrationarea(A::Array{Int64,2},M::Int64,N::Int64,redFac::Float64,calsize::Int64)
  cal_x, cal_y = min(calsize,M), min(calsize,N)
  ChosenSamples = 0
  NumberSamples = floor(Int64,M*N/redFac)
  if calsize > 0
    if cal_x*cal_y > NumberSamples
      @error "Number of calibration samples exceeds number of required total samples"
    end
    ChosenSamples = cal_x*cal_y
    A[floor(Int64,end/2)-floor(Int64,cal_x/2)+1:floor(Int64,end/2)-floor(Int64,cal_x/2)+cal_x,floor(Int64,end/2)-floor(Int64,cal_y/2)+1:floor(Int64,end/2)-floor(Int64,cal_y/2)+cal_y].=1
  end

  return A, ChosenSamples

end


text:
using AbstractMCMC
using DynamicPPL
using Random
using Turing

export AutoConditional

"""
    AutoConditional(sym)

A "pseudo-sampler" to use automatically extracted conditionals within `Gibbs`.  The conditionals are
extracted from the model trace at sampling step and allow dynamically changing models.  If you don't
need that, use `StaticConditional` instead.

`AutoConditional(:x)` will sample the Gibbs conditional of variable `x`, or equivalently, all
variables subsumed by `x`.

# Examples

```julia

@model test(x) = begin
    w ~ Dirichlet(2, 1.0)
    p ~ DiscreteNonParametric([0.3, 0.7], w)
    x ~ Bernoulli(p)
end

m = test(false)

sample(m, Gibbs(AutoConditional(:p), MH(:w)), 10)
```
"""
struct AutoConditional{S}
    AutoConditional(sym::Symbol) = new{sym}()
end

DynamicPPL.getspace(::AutoConditional{S}) where {S} = (S,)
DynamicPPL.alg_str(::AutoConditional) = "AutoConditional"
Turing.isgibbscomponent(::AutoConditional) = true

function Turing.Sampler(
    alg::AutoConditional,
    model::Turing.Model,
    s::DynamicPPL.Selector=DynamicPPL.Selector()
)
    state = Turing.Inference.SamplerState(Turing.VarInfo(model))
    return Turing.Sampler(alg, Dict{Symbol, Any}(), s, state)
end

function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::Turing.Model,
    spl::Turing.Sampler{<:AutoConditional{S}},
    N::Integer,
    transition;
    kwargs...
) where {S}
    graph = trackdependencies(model, spl.state.vi)
    conditioned_vn = DynamicPPL.VarName(S)
    values = sampled_values(graph)
    for (vn, cond) in conditionals(graph, conditioned_vn)
        updated = rand(rng, cond(values))
        spl.state.vi[vn] = [updated;]
    end

    return transition
end


text:
import Plotly
function showfn(space,i)
    geo = geometry(space)
    T = coordtype(geo)
    X = Dict{Int,T}()
    Y = Dict{Int,T}()
    Z = Dict{Int,T}()
    U = Dict{Int,T}()
    V = Dict{Int,T}()
    W = Dict{Int,T}()
    for sh in space.fns[i]
        chrt = chart(geo, cells(geo)[sh.cellid])
        nbd = CompScienceMeshes.center(chrt)
        vals = refspace(space)(nbd)
        x,y,z = cartesian(nbd)
        @show vals[sh.refid].value
        u,v,w = vals[sh.refid].value
        # @show x, y, z
        # @show u, v, w
        X[sh.cellid] = x
        Y[sh.cellid] = y
        Z[sh.cellid] = z
        U[sh.cellid] = get(U,sh.cellid,zero(T)) + sh.coeff * u
        V[sh.cellid] = get(V,sh.cellid,zero(T)) + sh.coeff * v
        W[sh.cellid] = get(W,sh.cellid,zero(T)) + sh.coeff * w
    end
    X = collect(values(X))
    Y = collect(values(Y))
    Z = collect(values(Z))
    U = collect(values(U))
    V = collect(values(V))
    W = collect(values(W))
    Plotly.cone(x=X,y=Y,z=Z,u=U,v=V,w=W)
end


text:
## general interface:
## - specify minimum size of sample and whether rolling window should
##   be applied
## - specify function to estimate and auxilliary data required for
##   estimator
##
## each estimator function signature should be:
## (discRet, begInd, endInd, auxData)

#######################
## general interface ##
#######################

function produceMoments(discRet::Timematr,
                        estimator::Function,
                        auxData::Tuple;
                        minObs=200, rolling=false)
    ## produce estimated moments for subsample, and return them as
    ## estimate for next day after subsample
    ##
    ## Inputs:
    ## 	discRet 		nObs x nAss Timematr of discrete returns
    ## 	estimator	estimator function that shall be applied to each
    ## 					subsample
    ## 	auxData 		1 x nAux tuple of auxiliary data used to impose
    ## 					structure within estimator
    ## 	minObs 		define (minimum) number of observations per
    ## 					estimate: first day that estimate is produced,
    ## 					and subsample size for rolling window approach
    ## 	rolling		Boolean value to indicate whether rolling window
    ## 					approach will be used or not
    ##
    ## Outputs:
    ## 	tuple (mu, covMatr, nextDay) with estimated moments for given
    ## 	subsample and scalar index of next day, where estimated
    ## 	moments will be used for to build a strategy at the next day

    ## general interface to moments estimation
    nObs = size(discRet, 1)

    ## specify how subsamples should be created
    produceSubsamples =
        Task(() -> AssetMgmt.subsampleProducer(nObs, minObs))

    dataArr = core(discRet)

    for subsample in produceSubsamples
        indices = subsample
        (mu, covMatr) = estimator(dataArr, indices..., auxData...)
        nextDay = indices[2]+1
        produce(mu, covMatr, nextDay)
    end
end

#######################
## empirical moments ##
#######################

function empiricalEstimator(discRet::Array{Float64, 2},
                            begInd::Int,
                            endInd::Int)
    ## for given data, apply estimator

    mus = mean(discRet[begInd:endInd, :], 1)
    covMatr = cov(discRet[begInd:endInd, :])

    return (mus, covMatr)
end

###################
## factor models ##
###################


text:
using QuasiMonteCarlo, Distributions, Test
using Test

#1D
lb = 0.0
ub = 5.0
n = 5
d = 1
QuasiMonteCarlo.sample(n,lb,ub,GridSample(0.1))
QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
QuasiMonteCarlo.sample(n,lb,ub,LatticeRuleSample())
QuasiMonteCarlo.sample(5,d,Cauchy())
QuasiMonteCarlo.sample(5,d,Normal(0,4))

@testset "1D" begin
    @testset "LowDiscrepancySample" begin
        s = QuasiMonteCarlo.sample(n, 0.0, 1.0, LowDiscrepancySample(2))
        @test isa(s, Vector{Float64})
        @test size(s) == (n,)
        @test s ≈ [0.5, 0.25, 0.75, 0.125, 0.625]

        s = QuasiMonteCarlo.sample(n, 0, 1, LowDiscrepancySample(2))
        @test isa(s, Vector{Float64})
        @test size(s) == (n,)
        @test s ≈ [0.5, 0.25, 0.75, 0.125, 0.625]

        s = QuasiMonteCarlo.sample(n, zero(Float32), one(Float32), LowDiscrepancySample(2))
        @test isa(s, Vector{Float32})
        @test size(s) == (n,)
        @test s ≈ [0.5, 0.25, 0.75, 0.125, 0.625] rtol = 1e-7
    end
end

#ND
lb = [0.0,0.0]
ub = [1.0,1.0]
n = 5
d = 2

@testset "GridSample" begin
    #GridSample{T}
    s = QuasiMonteCarlo.sample(n,lb,ub,GridSample([0.1,0.5]))
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "UniformSample" begin
    #UniformSample()
    s = QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "SobolSample" begin
    #SobolSample()
    s = QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "LHS" begin
    #LHS
    s = QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "LatticeRuleSample" begin
    #LatticeRuleSample()
    s = QuasiMonteCarlo.sample(n,lb,ub,LatticeRuleSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "LDS" begin
    #LDS
    s = QuasiMonteCarlo.sample(n,lb,ub,LowDiscrepancySample([2,3]))
    @test isa(s,Matrix{Float64})
    @test size(s) == (d, n)
    @test s[1, :] ≈ [0.5, 0.25, 0.75, 0.125, 0.625]
    @test s[2, :] ≈ [1/3, 2/3, 1/9, 4/9, 7/9]

    s = QuasiMonteCarlo.sample(n,Int.(lb),Int.(ub),LowDiscrepancySample([2,3]))
    @test isa(s,Matrix{Float64})
    @test size(s) == (d, n)
    @test s[1, :] ≈ [0.5, 0.25, 0.75, 0.125, 0.625]
    @test s[2, :] ≈ [1/3, 2/3, 1/9, 4/9, 7/9]

    s = QuasiMonteCarlo.sample(n,zeros(Float32, 2),ones(Float32, 2),LowDiscrepancySample([2,3]))
    @test isa(s,Matrix{Float32})
    @test size(s) == (d, n)
    @test s[1, :] ≈ [0.5, 0.25, 0.75, 0.125, 0.625] rtol = 1e-7
    @test s[2, :] ≈ [1/3, 2/3, 1/9, 4/9, 7/9] rtol = 1e-7
end

@testset "Distribution 1" begin
    #Distribution 1
    s = QuasiMonteCarlo.sample(n,d,Cauchy())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "Distribution 2" begin
    #Distribution 2
    s = QuasiMonteCarlo.sample(n,d,Normal(3,5))
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "generate_design_matrices" begin
    num_mat = 3
    Ms = QuasiMonteCarlo.generate_design_matrices(n,lb,ub,UniformSample(), num_mat)
    @test length(Ms) == num_mat
    A = Ms[1]
    B = Ms[2]
    @test isa(A, Matrix{typeof(A[1][1])}) == true
    @test size(A) == (d, n)
    @test isa(B, Matrix{typeof(B[1][1])}) == true
    @test size(B) == (d, n)
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 4af6dbd0-19b1-11eb-059b-c7ac5d8c4e0d
using FailureRepresentation

# ╔═╡ 5c706e30-19b1-11eb-18ad-5190eca820b1
begin
	using BSON
	using BSON: load, @load
	using CUDA
	using FailureRepresentation.SystemUnderTest
	using FailureRepresentation.Autoencoder
	using Flux
	using Flux: onecold
	using Flux.Data: DataLoader
	using NNlib
	Core.eval(Main, :(import NNlib)) # required for .load
	Core.eval(Main, :(import Flux)) # required for .load
end

# ╔═╡ 330f1d40-282d-11eb-29e0-f3affbba2711
using PlutoUI

# ╔═╡ d4fcfdae-1a14-11eb-0a76-5dfe76d1754a
using Plots; plotlyjs()

# ╔═╡ 2c1e5cc0-27c5-11eb-1e52-6121d656b973
using PyPlot, Seaborn

# ╔═╡ 27c76940-27c6-11eb-0ba2-c9213d619fb4
using MLBase

# ╔═╡ 5455e46e-27d3-11eb-1478-0f9879fc2215
using ScikitLearn

# ╔═╡ 4147f9a0-27d7-11eb-32de-a9d04fc45211
using ColorSchemes

# ╔═╡ 54164520-19b6-11eb-376c-b5246d05c3a6
MODEL = BSON.load("../models/sut.bson")[:m]

# ╔═╡ 30bd5680-19b2-11eb-3ed2-81593376e639
sut = SystemUnderTest.load("../models/sut.bson")

# ╔═╡ 1732d5e0-19b3-11eb-3b9d-8b5b61fc38ae
traindata, testdata = SystemUnderTest.getdata();

# ╔═╡ 1e4525a0-27a4-11eb-0422-9d0d7dc3dd9c
testdata

# ╔═╡ 24ef3930-19b3-11eb-0d9d-e3a0733ba1da
X, Y = rand(testdata, 3)

# ╔═╡ c3996450-1a14-11eb-2b50-934336a06d5a
Xc, Yc = deepcopy(X), deepcopy(Y);

# ╔═╡ b3418cb2-19b3-11eb-25b4-4b48a12ad82e
hcat([Autoencoder.img(X[:,i])' for i in 1:3])

# ╔═╡ c2f9dec0-27a3-11eb-0b14-75176c151264
hasfailure = onecold(sut(X)) .!= onecold(Y)

# ╔═╡ 8becbb82-19b3-11eb-2c35-b17e8dee1ed2
onecold(sut(X)) .- 1, onecold(Y) .- 1 # notice off-by-one

# ╔═╡ 708c22b0-19b1-11eb-3ee7-f32bf48ebf1f
autoencoder = Autoencoder.load("../models/autoencoder.bson")

# ╔═╡ 1cef6c80-19b5-11eb-1f69-7356f410b1fc
hcat([Autoencoder.img(autoencoder(X[:,i]))' for i in 1:3])

# ╔═╡ 3cf13db0-19b5-11eb-081b-ed8141b6b8d7
onecold(sut(autoencoder(X))) .- 1, onecold(Y) .- 1 # notice off-by-one

# ╔═╡ f4fbfc40-19b2-11eb-0ee1-0b3e9ece6008
Autoencoder.sample(autoencoder, 5)

# ╔═╡ 8cde0240-1a14-11eb-0274-090f1dcdeb0d
softmax(sut(X[:,1]))

# ╔═╡ b6901270-1a16-11eb-09a2-c1073bff3f8f
md"""
$$\Delta\operatorname{softmax}(\mathbf y, \mathbf{\hat y}) = \max\left(\operatorname{softmax}(\mathbf{\hat y})\right) - \operatorname{softmax}(\mathbf{\hat y})^{\left[ \operatorname{argmax} \mathbf y \right]}$$
"""

# ╔═╡ 5d045420-2833-11eb-1f67-59da9ee8ef92
md"""
-  $$\Delta\operatorname{softmax} \in [0, 1]$$
- When $$\Delta\operatorname{softmax}$$ is large, severe failure
- When $$\Delta\operatorname{softmax}$$ is small, likely failure
- When $$\Delta\operatorname{softmax}$$ is zero, no failure
"""

# ╔═╡ 95bd6740-1a17-11eb-3372-21bac6614f53
function Δsoftmax(model, x::Vector, y)
	ŷ = model(x)
	return maximum(softmax(ŷ)) - softmax(ŷ)[argmax(y)]
end;

# ╔═╡ a9de0db0-1a2b-11eb-2144-e5025a1a9e5e
function Δsoftmax(model, x::Matrix, y)
	ŷ = model(x)
	softmaxŷ = softmax(ŷ)
	max_softmax = mapslices(maximum, softmaxŷ, dims=1)'
	true_yᵢ = mapslices(argmax, y, dims=1)
	true_softmax = [softmaxŷ[r,c] for (r,c) in zip(true_yᵢ, 1:size(y,2))]
	return vec(max_softmax - true_softmax)
end;

# ╔═╡ 3da097b0-2833-11eb-3566-11862e85da01
Δsoftmax(sut, X, Y)

# ╔═╡ 1f747b72-2834-11eb-0015-dd02035ce6bc
FailureRepresentation.Adversary.bce(0.5, 1) .+ log.(Δsoftmax(sut, X, Y) .+ eps(1.0))

# ╔═╡ 8fd9b780-2835-11eb-2698-ed54bb00d2c7
dl = Flux.Data.DataLoader([1, 1], [2, 2], [3, 3])

# ╔═╡ 974bb132-2835-11eb-29fa-d342c6b4712b
for (x,y) in dl
	@info x,y
end

# ╔═╡ 009237b0-1a15-11eb-3e6c-9ff29ea9a32a
begin
	idx = 1
	Plots.bar(0:9, Y[:,idx], xticks=0:9, fillalpha=0.6, label="true")
	Plots.bar!(0:9, softmax(sut(X[:,idx])), fillalpha=0.6, label="prediction")
end

# ╔═╡ e083cc00-1ac7-11eb-103a-13fb6fb3d586
Autoencoder.img(X[:,idx])'

# ╔═╡ 0e50baa0-1ad5-11eb-1603-cba64a9c41e7
Δsoftmax(sut, X[:, idx], Y[:, idx])

# ╔═╡ 99924160-1a35-11eb-0df0-ff718845bd4b
Xi, Yi = rand(testdata, 10)

# ╔═╡ 7b8cf5c0-1a30-11eb-149c-5fbf3a4d9a12
Δsoftmax(sut, Xi[:, 1], Yi[:, 1])

# ╔═╡ 769264b0-1a30-11eb-3742-7f8badcea17f
Δsoftmax(sut, Xi, Yi)

# ╔═╡ 9f53f2fe-1a35-11eb-241a-15d6dc131650
Δsoftmax(sut, Xi, Yi)

# ╔═╡ 2059e2d0-1ad5-11eb-21c5-abf90763e4b9
md"""
---
"""

# ╔═╡ 805fa7c0-1a15-11eb-2200-a166cafe9431
# function Δsoftmax_expanded(model, x, y)
# 	probs = softmax(model(x))
# 	ŷᵢₙₓ = onecold(probs)
# 	yᵢₙₓ = onecold(y)
# 	max_softmax = [probs[r,c] for (r,c) in zip(onecold(probs), 1:size(probs,2))]
# 	true_softmax = [probs[r,c] for (r,c) in zip(onecold(Yi), 1:size(probs,2))]
# 	return max_softmax - true_softmax
# end;

# ╔═╡ 4b5dc850-27c5-11eb-015b-a91f5e61a3cc
function confusion(testdata, model)
	confY = Int[]
	confY′ = Int[]
	for (x,y) in testdata
		push!(confY, onecold(y)...)
		push!(confY′, onecold(softmax(model(x)))...)
	end
	return confY, confY′
end

# ╔═╡ 5c81b510-27cf-11eb-2081-e7eb2f26ee09
function confusion_adversary(testdata, model, threshold=0.5)
	confY = Int[]
	confY′ = Int[]
	confY′_prob = Float64[]
	for (x,y) in testdata
		push!(confY, y...)
		push!(confY′, FailureRepresentation.Adversary.predict(model, x, threshold)...)
		push!(confY′_prob, model(x)...)
	end
	return confY, confY′, confY′_prob
end

# ╔═╡ 9038c9d0-2869-11eb-21e5-b5fe26e0435f
function confusion_random(testdata)
	confY = Int[]
	confY′ = Int[]
	for (x,y) in testdata
		push!(confY, y...)
		push!(confY′, rand() .>= 0.5)
	end
	return confY, confY′
end

# ╔═╡ 607c1fc0-27c5-11eb-2b0c-db7dd00ada2c
# confY, confY′ = confusion(x->sut(autoencoder(x)))
confY, confY′ = confusion(testdata, sut)

# ╔═╡ 41d1f9e0-27c6-11eb-22eb-cb7988e1c216
begin
	conf = confusmat(10, confY, confY′)
	# conf = conf ./ sum(conf, dims=2) # normalize per class
end

# ╔═╡ 345edf20-27c7-11eb-14c7-ed2a1f653528
correctrate(confY, confY′)

# ╔═╡ 74419690-27c8-11eb-2bef-95b760034b9f
FailureRepresentation.SystemUnderTest.accuracy(testdata, sut)

# ╔═╡ 87c13220-27c8-11eb-33a9-7b92300eb881
M = conf

# ╔═╡ 02e5174e-27c9-11eb-0f16-ffa235e1c87c
precision(M, i) = M[i,i] / sum(M[:,i])

# ╔═╡ 6e904330-27c9-11eb-1ba2-afae7afad48b
recall(M, i) = M[i,i] / sum(M[i,:])

# ╔═╡ 8369c970-27c9-11eb-202d-cbc34e5eb985
mean(precision(M, i) for i in 1:10), mean(recall(M, i) for i in 1:10)

# ╔═╡ 3bad2432-27c7-11eb-26bb-e51b66fa7698
errorrate(confY, confY′)

# ╔═╡ 2fffb5a0-27c5-11eb-28bc-35c520a05b32
begin
	close("all")
	figure(figsize=(4,3))
	Seaborn.heatmap(conf, annot=true, cmap="viridis",
		            xticklabels=0:9, yticklabels=0:9, fmt="d",
		            annot_kws=Dict("size"=>8))
	title("Confusion matrix for system \$\\mathcal{S}\$")
	xlabel("predicted labels")
	ylabel("true labels")
	gcf()
end

# ╔═╡ adc62590-27da-11eb-1c83-c3065a5b1c54
PyPlot.savefig("confusion_sut.pgf", bbox_inches="tight")

# ╔═╡ c6a09110-27ce-11eb-196e-4365b178198a
adversary = BSON.load("../models/adversary_10.bson")[:adversary]

# ╔═╡ 906c4e70-27d0-11eb-16d3-15b462cf88e8
encoder = BSON.load("../models/encoder.bson")[:encoder]

# ╔═╡ fb4d27a0-27df-11eb-11dc-ddbfffd9fc84
decoder = BSON.load("../models/decoder.bson")[:decoder]

# ╔═╡ f5a06d50-27ce-11eb-17b4-c1599c1d93db
function convert_data(testdata, sut, encoder)
    x̃ = reduce(hcat, [encoder(x) for (x,_) in testdata])
    ỹ = reduce(vcat, [onecold(sut(x)) .!= onecold(y) for (x,y) in testdata])
    return DataLoader(x̃, ỹ)
end

# ╔═╡ e5c0b4d0-27ce-11eb-330b-2dfe47120ede
testdata_tilde = convert_data(testdata, sut, encoder)

# ╔═╡ 20b3f42e-27cf-11eb-00e7-73c83aa66e00
conf𝒜, conf𝒜′, conf𝒜′_prob = confusion_adversary(testdata_tilde, adversary)

# ╔═╡ 4bf84f10-27cf-11eb-0435-a12463478ae0
M𝒜 = confusmat(2, conf𝒜 .+ 1, conf𝒜′ .+ 1)

# ╔═╡ 15da21f0-27d0-11eb-025b-e7f9374b020b
begin
	close("all")
	PyPlot.figure(figsize=(4,3))
	Seaborn.heatmap(M𝒜, annot=true, cmap="viridis",
		            xticklabels=["not failure", "failure"],
		            yticklabels=["not failure", "failure"],
		            fmt="d", annot_kws=Dict("size"=>12))
	title("Confusion matrix for adversary \$\\mathcal{A}\$")
	xlabel("predicted failures")
	ylabel("true failures")
	gcf()
end

# ╔═╡ a1d0c69e-27da-11eb-163b-ab92a7859f62
PyPlot.savefig("confusion_adversary.pgf", bbox_inches="tight")

# ╔═╡ 2f42f450-27d0-11eb-1f41-cf1df27badc2
mean(precision(M𝒜, i) for i in 1:2), mean(recall(M𝒜, i) for i in 1:2)

# ╔═╡ bae39552-27d0-11eb-1c3b-cb0bfcacdb81
precision(M𝒜, 2), recall(M𝒜, 2)

# ╔═╡ 82a942e0-2869-11eb-17f2-5550a7d223a0
md"""
## Random baseline
"""

# ╔═╡ 867c6aee-2869-11eb-15bb-2b5934eaba9f
conf𝒜_r, conf𝒜′_r = confusion_random(testdata_tilde)

# ╔═╡ b8f05f50-2869-11eb-158b-47205e98fd88
M𝒜_r = confusmat(2, conf𝒜_r .+ 1, conf𝒜′_r .+ 1)

# ╔═╡ 65a19ee0-2869-11eb-0f9d-6d7759d79d9f
begin
	close("all")
	PyPlot.figure(figsize=(4,3))
	Seaborn.heatmap(M𝒜_r, annot=true, cmap="viridis",
		            xticklabels=["not failure", "failure"],
		            yticklabels=["not failure", "failure"],
		            fmt="d", annot_kws=Dict("size"=>8))
	title("Confusion matrix for adversary \$\\mathcal{A}\$")
	xlabel("predicted failures")
	ylabel("true failures")
	gcf()
end

# ╔═╡ 789ea780-286a-11eb-28dc-2185259857bd
precision(M𝒜_r, 2), recall(M𝒜_r, 2)

# ╔═╡ 798a4000-286a-11eb-1f12-e30338eb7ead
md"""
## ROC curve scikitlearn
"""

# ╔═╡ 268df680-27d4-11eb-175c-f79f746fbab7
@sk_import metrics: roc_curve

# ╔═╡ 5bc8ac00-27d4-11eb-31aa-edb430062b1b
@sk_import metrics: auc

# ╔═╡ 7e894f10-27d4-11eb-2b28-dbed5929b12b
begin
	rX, rY = [], []
	for (x,y) in testdata_tilde
		push!(rX, x...)
		push!(rY, y...)
	end
end

# ╔═╡ 96d3e030-27d4-11eb-39ea-db39558701a6
rY

# ╔═╡ 469ee2e2-27d4-11eb-3e68-39a1c8de84e3
fpr, tpr, _ = roc_curve(conf𝒜, conf𝒜′_prob)

# ╔═╡ 460879d0-27d5-11eb-2495-7920d55f677c
auc(fpr, tpr)

# ╔═╡ 15281992-27dc-11eb-2e1e-1b4d5028142b
md"""
## Threshold sweep
"""

# ╔═╡ 1985ce60-27dc-11eb-3aa6-f55d06695c2d
begin
	precision_t = []
	recall_t = []
	thresholds = 0:0.05:1
	for thresh in thresholds
		conf𝒜_t, conf𝒜′_t, conf𝒜′_prob_t =
			confusion_adversary(testdata_tilde, adversary, thresh)
		M𝒜_t = confusmat(2, conf𝒜_t .+ 1, conf𝒜′_t .+ 1)
		push!(precision_t, precision(M𝒜_t, 2))
		push!(recall_t, recall(M𝒜_t, 2))
	end
end

# ╔═╡ e4bf4200-27dc-11eb-051f-8b50f7f4bdd8
begin
	close("all")
	figure(figsize=(4,3))
	PyPlot.plot(thresholds, precision_t, color="blue", marker=".")
	PyPlot.plot(thresholds, recall_t, color="red", marker=".")
	legend(["Precision", "Recall"])
	title("Prediction threshold sweep: precision and recall")
	xlabel("prediction threshold: \$\\hat y \\ge {threshold}\$")
	ylabel("precision and recall")
	PyPlot.plot([0.5005, 0.5005], [-0.01, recall_t[11]], linestyle=":", color="gray")
	ylim([-0.01, 1.05])
	gcf()
end

# ╔═╡ ae24ea8e-27de-11eb-216e-55bf2ed07117
PyPlot.savefig("threshold_sweep.pgf", bbox_inches="tight")

# ╔═╡ 1cad742e-27dc-11eb-1d62-ed5ef14f01a0
md"""
## ROC curves
"""

# ╔═╡ 430dc580-27d7-11eb-2edd-81e1104ab57b
ColorSchemes.Blues_9

# ╔═╡ de03e060-27d7-11eb-19e4-f303bb9078f3
rgb(c) = (red(c), green(c), blue(c))

# ╔═╡ edbb0240-27d7-11eb-1c26-e3a457165160
rgb(ColorSchemes.Reds_9[1])

# ╔═╡ 384e63a0-27d9-11eb-1540-7731666f9aac
# using PGFPlots

# ╔═╡ dde1c780-27d9-11eb-2345-b9067ba467b7
begin
	matplotlib.rc("font", family=["serif"]) # sans-serif keyword not supported
	# matplotlib.rc("font", serif=["Helvetica"]) # matplotlib.rc("font", serif=["Palatino"])
	matplotlib.rc("text", usetex=true)
	matplotlib.rc("pgf", rcfonts=false)
end

# ╔═╡ c0cf6762-27d4-11eb-3a68-e1fc72404bd7
begin
	close("all")
	figure(figsize=(4,3))
	legends = []
	for i in 1:10
		color = i >= 10 ? "black" : rgb(ColorSchemes.Blues_9[i])
		adversary_i = BSON.load("../models/adversary_$i.bson")[:adversary]
		conf𝒜_i, conf𝒜′_i, conf𝒜′_prob_i =
			confusion_adversary(testdata_tilde, adversary_i)
		fprᵢ, tprᵢ, _ = roc_curve(conf𝒜_i, conf𝒜′_prob_i)
		PyPlot.plot(fprᵢ, tprᵢ, color=color)
		push!(legends, "AUC\$_{$i}\$ = \$$(round(auc(fprᵢ, tprᵢ), digits=3))\$")
	end
	ylabel("true positive rate")
	xlabel("false positive rate")
	title("Receiver operating characteristic (ROC) curve per iteration \$t \\in T\$")
	legend(legends, loc="center left", bbox_to_anchor=(1.04, 0.5))
	PyPlot.plot(0:1, 0:1, linestyle="--", color="gray")
	gcf()
end

# ╔═╡ bc9ee040-27d8-11eb-2ae2-fd3ad839dffb
PyPlot.savefig("roc.pgf", bbox_inches="tight")

# ╔═╡ 880aede0-27df-11eb-342c-af78db5ae8bb
md"""
## Qualitative analysis
Most likely failure via the prediction, but wasn't a failure.
"""

# ╔═╡ 9207b11e-27df-11eb-262b-83039547db91
arg = argmax([maximum(adversary(x)) for (x,y) in testdata_tilde])

# ╔═╡ c06496a0-27df-11eb-0557-35927114712c
arg_x, arg_y = collect(testdata_tilde)[arg]

# ╔═╡ e9f2b600-27df-11eb-3190-d3a359b8b5cb
Autoencoder.img(decoder(vec(arg_x)))'

# ╔═╡ 314cdad0-27e0-11eb-261b-359761d4268b
Flux.flatten(collect(testdata))

# ╔═╡ 1a62c410-27e0-11eb-14ab-339edb9da195
function select_top_k(model, testdata, k=5, rev=true)
    return partialsortperm([model(x)[1] for (x,y) in testdata], 1:k, rev=rev)
end

# ╔═╡ b41442a0-2830-11eb-03d0-ed6c96e9dde9
function select_bottom_k_false_negatives(model, testdata, k=5, rev=false)
	idx = []
	predicted_value = []
	for (i,(x,y)) in enumerate(testdata)
		# if false negative
		if y[1] && !FailureRepresentation.Adversary.predict(model, x)[1]
			push!(idx, i)
			push!(predicted_value, model(x)[1])
		end
	end
	perm = partialsortperm(predicted_value, 1:k, rev=rev)
	return idx[perm], predicted_value[perm]
end

# ╔═╡ 2117c860-287a-11eb-11f2-5304606bcdf3
function select_top_k_true_positives(model, testdata, k=5, rev=true)
	idx = []
	predicted_value = []
	for (i,(x,y)) in enumerate(testdata)
		# if true positive
		if y[1] && FailureRepresentation.Adversary.predict(model, x)[1]
			push!(idx, i)
			push!(predicted_value, model(x)[1])
		end
	end
	perm = partialsortperm(predicted_value, 1:k, rev=rev)
	return idx[perm], predicted_value[perm]
end

# ╔═╡ ab2aeb90-287a-11eb-0a90-ff7e62c0f950
md"""
## False negatives
"""

# ╔═╡ 4b99eee0-2831-11eb-20e4-cbae0c67e4fb
bot_false_negs = select_bottom_k_false_negatives(adversary, testdata_tilde, 10)

# ╔═╡ ce96ccf0-2831-11eb-3f91-4bbe163284b4
false_neg_decoded_images = Autoencoder.img.([decoder(vec(args_x)) for (args_x, args_y) in collect(testdata_tilde)[bot_false_negs[1]]])'

# ╔═╡ 4499ad50-2832-11eb-1316-a3370fe12816
false_neg_true_failures = [args_y[1] for (args_x, args_y) in collect(testdata_tilde)[bot_false_negs[1]]]

# ╔═╡ 16d1a010-287a-11eb-1ef2-eb3e698fe85f
md"""
## True positives
"""

# ╔═╡ 33caec30-287a-11eb-099b-5b6571982daf
top_true_pos = select_top_k_true_positives(adversary, testdata_tilde, 10)

# ╔═╡ 4af4cb10-287a-11eb-3718-f7d59a2c0eed
true_pos_decoded_images = Autoencoder.img.([decoder(vec(args_x)) for (args_x, args_y) in collect(testdata_tilde)[top_true_pos[1]]])'

# ╔═╡ cdf44130-287a-11eb-13dd-796a7a81dd8a
md"""
## Likely failures
"""

# ╔═╡ e02a6c70-2830-11eb-08d9-631b8f5c8166
testdata_tilde

# ╔═╡ 9c05b770-27e0-11eb-2116-55c62bb8e032
args = select_top_k(adversary, testdata_tilde, 10)

# ╔═╡ 89daa100-2830-11eb-0ab0-09728faf861b
args_bot = select_top_k(adversary, testdata_tilde, 10, false)

# ╔═╡ 9ab0ee80-2830-11eb-15eb-6dcee557eadd
[minimum(adversary(x)) for (x,y) in testdata_tilde][args_bot]

# ╔═╡ a4de85c0-2830-11eb-181a-079c0c560cb6
true_failures_bot = [args_y[1] for (args_x, args_y) in collect(testdata_tilde)[args_bot]]

# ╔═╡ b4aabd40-27e3-11eb-1d99-4bb53010e540
[maximum(adversary(x)) for (x,y) in testdata_tilde][args]

# ╔═╡ b31de312-27e0-11eb-0a4e-cde578822d68
decoded_images = Autoencoder.img.([decoder(vec(args_x)) for (args_x, args_y) in collect(testdata_tilde)[args]])'

# ╔═╡ df09c890-27e0-11eb-1509-afdf248cd717
true_failures = [args_y[1] for (args_x, args_y) in collect(testdata_tilde)[args]]

# ╔═╡ 7a323e12-27e1-11eb-3720-a35f5fff7948
function flatten_testdata(testdata)
	X, Y = [], []
	for (x,y) in testdata
		@show size(x), size(y)
		for i in 1:size(y,2)
			push!(X, x[:,i])
			push!(Y, y[:,i])
		end
	end
	@info size(X), size(Y)
	return Flux.Data.DataLoader(X, Y)
end

# ╔═╡ 90073b4e-27e1-11eb-1b41-9d02c393df7c
ftestdata = flatten_testdata(testdata)

# ╔═╡ f7aada50-2831-11eb-2318-0f64039d365d
false_neg_true_inputs = [x[1] for x in first.(collect(ftestdata))[bot_false_negs[1]]]

# ╔═╡ eb076110-2831-11eb-2151-6ba1beee687c
false_neg_true_images = Autoencoder.img.([x for x in false_neg_true_inputs])'

# ╔═╡ 83009a90-287d-11eb-364a-014efae85a3a
fn_true_labels = [onecold(y[1]) - 1 for y in last.(collect(ftestdata))[bot_false_negs[1]]]

# ╔═╡ 89381872-287d-11eb-0991-bdd241001a4c
fn_sut_labels = [onecold(sut(x[1])) - 1 for x in first.(collect(ftestdata))[bot_false_negs[1]]]

# ╔═╡ 5e6a0660-287a-11eb-0308-874829b878aa
true_pos_true_inputs = [x[1] for x in first.(collect(ftestdata))[top_true_pos[1]]]

# ╔═╡ 678d3280-287a-11eb-2c43-bd919979d395
true_pos_true_images = Autoencoder.img.([x for x in true_pos_true_inputs])'

# ╔═╡ 675fb282-287d-11eb-0c1d-7d21c31b3ae6
tp_true_labels = [onecold(y[1]) - 1 for y in last.(collect(ftestdata))[top_true_pos[1]]]

# ╔═╡ 72e42730-287d-11eb-26d4-e9cbe9b7b2cc
tp_sut_labels = [onecold(sut(x[1])) - 1 for x in first.(collect(ftestdata))[top_true_pos[1]]]

# ╔═╡ f18d3460-27e1-11eb-0351-f9b44125b4bc
true_labels = [onecold(y[1]) - 1 for y in last.(collect(ftestdata))[args]]

# ╔═╡ 885a17f0-27e2-11eb-2720-1d13935ac861
sut_labels = [onecold(sut(x[1])) - 1 for x in first.(collect(ftestdata))[args]]

# ╔═╡ 563506e0-27e2-11eb-154d-cb3af1504823
true_inputs = [x[1] for x in first.(collect(ftestdata))[args]]

# ╔═╡ 502164b0-27e2-11eb-1201-b93ded002d96
true_images = Autoencoder.img.([x for x in true_inputs])'

# ╔═╡ 90ae7b00-27e5-11eb-0991-532c3e8e3941
redimg(x::Vector) = reshape([RGB(clamp(i, 0, 1), 0, 0) for i in x], 28, 28)

# ╔═╡ f6f1ae00-27e5-11eb-2538-6dcf631a43b3
greenimg(x::Vector) = reshape([RGB(0, clamp(i, 0, 1), 0) for i in x], 28, 28)

# ╔═╡ 34b11b30-2832-11eb-0cb5-775f817a2d7c
false_neg_failure_indicators = [y ? greenimg(fill(0.4, 28^2)) : redimg(fill(0.4, 28^2)) for y in false_neg_true_failures]'

# ╔═╡ 291ceb80-27e4-11eb-32d3-69542e48e4e7
failure_indicators = [y ? greenimg(fill(0.4, 28^2)) : redimg(fill(0.4, 28^2)) for y in true_failures]'

# ╔═╡ a51cbc50-27e5-11eb-11d0-317e3a30e5ea
redimg(fill(0.4, 28^2)), greenimg(fill(0.4, 28^2))

# ╔═╡ 8c6f5540-27e6-11eb-226b-3d61418cc764
hcat([Gray.(reshape(args_x, 8, 8)) for (args_x, args_y) in collect(testdata_tilde)[args]]...)

# ╔═╡ a8bb0e30-27e4-11eb-1578-53806dd0f0e9
Autoencoder.img(decoder(randn(64)))'

# ╔═╡ 32323f80-27e5-11eb-1581-3537f89b6549
hcat((true_images .- decoded_images)...)

# ╔═╡ 54351b60-27e6-11eb-3898-15504bbff380
Gray.(reshape(collect(testdata_tilde)[6][1], 8, 8))

# ╔═╡ 5c03be90-27e7-11eb-0d07-6529c6633b4f
function onehotvector(i, n, v=1)
	X = zeros(n)
	X[i] = v
	return X
end

# ╔═╡ 0a565d00-2832-11eb-3110-a9c7d74cb153
false_neg_feature_images = Autoencoder.img.([decoder(onehotvector(argmax(softmax(vec(args_x))), 64, 10)) for (args_x, args_y) in collect(testdata_tilde)[bot_false_negs[1]]])'

# ╔═╡ 15abd26e-2832-11eb-1cae-818f96fce3f9
false_negatives = hcat(vcat.(false_neg_true_images, false_neg_feature_images, false_neg_decoded_images)...)

# ╔═╡ 70830d10-287a-11eb-269c-c3fa1eaeb9e2
true_pos_feature_images = Autoencoder.img.([decoder(onehotvector(argmax(softmax(vec(args_x))), 64, 10)) for (args_x, args_y) in collect(testdata_tilde)[top_true_pos[1]]])'

# ╔═╡ 05517e02-287a-11eb-167c-0935bd9a2a70
true_positives = hcat(vcat.(true_pos_true_images, true_pos_feature_images, true_pos_decoded_images)...)

# ╔═╡ 4999e350-27e8-11eb-143b-dd46d2a2a8d9
feature_images = Autoencoder.img.([decoder(onehotvector(argmax(softmax(vec(args_x))), 64, 10)) for (args_x, args_y) in collect(testdata_tilde)[args]])'

# ╔═╡ eddba8a0-27e2-11eb-3c2b-456abfb6d91c
likely_failures = hcat(vcat.(true_images, feature_images, decoded_images, failure_indicators)...)

# ╔═╡ dcb73e50-27e6-11eb-3f42-35c9f67c96da
Autoencoder.img(decoder(onehotvector(1,64,10)))'

# ╔═╡ 9bd64200-27eb-11eb-2c81-f7e1bb552e84
function twohotvector(i, j, n, v=1)
	X = zeros(n)
	X[i] = v
	X[j] = v
	return X
end

# ╔═╡ 7a474bb2-27ec-11eb-3a1b-4d6f11cc8042
function threehotvector(i, j, k, n, v=1)
	X = zeros(n)
	X[i] = v
	X[j] = v
	X[k] = v
	return X
end

# ╔═╡ 6ff44460-27e7-11eb-371e-9dddb1d9a76c
onehotvector(1,10)

# ╔═╡ 8adfbf20-27e7-11eb-0cec-91f326f42e55
ex = onehotvector(1,64)

# ╔═╡ 0ac1d330-27e9-11eb-148e-156e7d670757
hcat(vcat.([[Autoencoder.img(decoder(onehotvector(i, 64, 10)))' for i in (8j - 7):8j] for j in 1:8]...)...)

# ╔═╡ 734e20f0-27eb-11eb-11fb-d51498c1592f
Autoencoder.img(decoder(twohotvector(1, 6, 64, 10)))'

# ╔═╡ c1f48f00-27eb-11eb-13af-599491b28241
argmax([[adversary(twohotvector(i, j, 64, 10)) for i in 1:64] for j in 1:64])

# ╔═╡ 0bd585c0-27ec-11eb-32ea-f33ef2084890
combinations = vec([(i,j) for i in 1:64, j in 1:64])

# ╔═╡ 81a56f40-27ec-11eb-1bf8-b30e0f259dc4
combinations3 = vec([(i,j,k) for i in 1:64, j in 1:64, k in 1:64])

# ╔═╡ 2e574b10-27ec-11eb-26ea-db2bfc5278fe
max2hot = argmax([adversary(twohotvector(i, j, 64, 10)) for (i,j) in combinations])

# ╔═╡ 67663a10-27ec-11eb-0173-27217c95d913
adversary(twohotvector(combinations[max2hot]..., 64, 1))

# ╔═╡ eb2cc630-27eb-11eb-0496-831311ac09e9
Autoencoder.img(decoder(twohotvector(combinations[max2hot]..., 64, 10)))'

# ╔═╡ 8a935a92-27ec-11eb-3508-d9ca334ea4b5
max3hot = argmax([adversary(threehotvector(i, j, k, 64, 10)) for (i,j,k) in combinations3])

# ╔═╡ 96b791b0-27ec-11eb-27a4-bf55cb102ee2
Autoencoder.img(decoder(threehotvector(combinations3[max3hot]..., 64, 10)))'

# ╔═╡ 292d3070-27ea-11eb-088b-efc17f03f787
maximum(adversary(onehotvector(i, 64, 1)) for i in 1:64)

# ╔═╡ 273869e0-282d-11eb-12ae-e90ddca72641
@bind intensity Slider(1:100, default=10, show_value=true)

# ╔═╡ 47fee610-27ea-11eb-0a2c-c1b0c93f4446
encoding_map = hcat(vcat.([[Autoencoder.img(decoder(onehotvector(i, 64, intensity*adversary(onehotvector(i,64))[1])))' ./ intensity for i in (8j - 7):8j] for j in 1:8]...)...)

# ╔═╡ 29974c00-282e-11eb-044e-a972ed21be77
begin
	close("all")
	figure(figsize=(4,3))
	ticks8 = range(28/2, stop=8*28-(28/2), length=8)
	Seaborn.heatmap(convert(Array{Float32, 2}, encoding_map),
		            cmap="viridis", xticklabels=1:8, yticklabels=1:8)
	xticks(ticks8)
	yticks(ticks8)
	title("Likelihood sweep over low-dimensional features")
	xlabel("feature index \$i\$")
	ylabel("feature index \$j\$")
	gcf()
end

# ╔═╡ 0e986a50-2870-11eb-2bdc-ef8bd3938307
PyPlot.savefig("feature-likelihood.pgf", bbox_inches="tight")

# ╔═╡ 5936b910-282d-11eb-107f-69a8742adfa0
@bind likelihood_intensity Slider(1:1)

# ╔═╡ 632eb3a0-2828-11eb-1933-7b74bf88bcdc
begin
	close("all")
	figure(figsize=(4,3))
	Seaborn.heatmap(reshape([adversary(encoder(onehotvector(i, 28^2, likelihood_intensity)))[1] for i in 1:28^2], 28, 28)', cmap="viridis", xticklabels=0:7:28, yticklabels=0:7:28)
	title("Failure likelihood per pixel")
	xticks(0:7:28)
	yticks(0:7:28)
	xlabel("\$x\$-pixel")
	ylabel("\$y\$-label")
	gcf()
end

# ╔═╡ 3584b5b0-2870-11eb-3d47-7762282618a9
PyPlot.savefig("failure-likelihood.pgf", bbox_inches="tight")

# ╔═╡ ab1efe50-2872-11eb-0a14-f900c5907e05
md"""
## SUT metrics
"""

# ╔═╡ Cell order:
# ╠═4af6dbd0-19b1-11eb-059b-c7ac5d8c4e0d
# ╠═5c706e30-19b1-11eb-18ad-5190eca820b1
# ╠═330f1d40-282d-11eb-29e0-f3affbba2711
# ╠═54164520-19b6-11eb-376c-b5246d05c3a6
# ╠═30bd5680-19b2-11eb-3ed2-81593376e639
# ╠═1732d5e0-19b3-11eb-3b9d-8b5b61fc38ae
# ╠═1e4525a0-27a4-11eb-0422-9d0d7dc3dd9c
# ╠═24ef3930-19b3-11eb-0d9d-e3a0733ba1da
# ╠═c3996450-1a14-11eb-2b50-934336a06d5a
# ╠═b3418cb2-19b3-11eb-25b4-4b48a12ad82e
# ╠═c2f9dec0-27a3-11eb-0b14-75176c151264
# ╠═8becbb82-19b3-11eb-2c35-b17e8dee1ed2
# ╠═708c22b0-19b1-11eb-3ee7-f32bf48ebf1f
# ╠═1cef6c80-19b5-11eb-1f69-7356f410b1fc
# ╠═3cf13db0-19b5-11eb-081b-ed8141b6b8d7
# ╠═f4fbfc40-19b2-11eb-0ee1-0b3e9ece6008
# ╠═8cde0240-1a14-11eb-0274-090f1dcdeb0d
# ╠═d4fcfdae-1a14-11eb-0a76-5dfe76d1754a
# ╟─b6901270-1a16-11eb-09a2-c1073bff3f8f
# ╟─5d045420-2833-11eb-1f67-59da9ee8ef92
# ╠═95bd6740-1a17-11eb-3372-21bac6614f53
# ╠═7b8cf5c0-1a30-11eb-149c-5fbf3a4d9a12
# ╠═a9de0db0-1a2b-11eb-2144-e5025a1a9e5e
# ╠═769264b0-1a30-11eb-3742-7f8badcea17f
# ╠═3da097b0-2833-11eb-3566-11862e85da01
# ╠═1f747b72-2834-11eb-0015-dd02035ce6bc
# ╠═8fd9b780-2835-11eb-2698-ed54bb00d2c7
# ╠═974bb132-2835-11eb-29fa-d342c6b4712b
# ╠═009237b0-1a15-11eb-3e6c-9ff29ea9a32a
# ╠═e083cc00-1ac7-11eb-103a-13fb6fb3d586
# ╠═0e50baa0-1ad5-11eb-1603-cba64a9c41e7
# ╠═99924160-1a35-11eb-0df0-ff718845bd4b
# ╠═9f53f2fe-1a35-11eb-241a-15d6dc131650
# ╟─2059e2d0-1ad5-11eb-21c5-abf90763e4b9
# ╠═805fa7c0-1a15-11eb-2200-a166cafe9431
# ╠═2c1e5cc0-27c5-11eb-1e52-6121d656b973
# ╠═4b5dc850-27c5-11eb-015b-a91f5e61a3cc
# ╠═5c81b510-27cf-11eb-2081-e7eb2f26ee09
# ╠═9038c9d0-2869-11eb-21e5-b5fe26e0435f
# ╠═607c1fc0-27c5-11eb-2b0c-db7dd00ada2c
# ╠═27c76940-27c6-11eb-0ba2-c9213d619fb4
# ╠═41d1f9e0-27c6-11eb-22eb-cb7988e1c216
# ╠═345edf20-27c7-11eb-14c7-ed2a1f653528
# ╠═74419690-27c8-11eb-2bef-95b760034b9f
# ╠═87c13220-27c8-11eb-33a9-7b92300eb881
# ╠═02e5174e-27c9-11eb-0f16-ffa235e1c87c
# ╠═6e904330-27c9-11eb-1ba2-afae7afad48b
# ╠═8369c970-27c9-11eb-202d-cbc34e5eb985
# ╠═3bad2432-27c7-11eb-26bb-e51b66fa7698
# ╠═2fffb5a0-27c5-11eb-28bc-35c520a05b32
# ╠═adc62590-27da-11eb-1c83-c3065a5b1c54
# ╠═c6a09110-27ce-11eb-196e-4365b178198a
# ╠═906c4e70-27d0-11eb-16d3-15b462cf88e8
# ╠═fb4d27a0-27df-11eb-11dc-ddbfffd9fc84
# ╠═e5c0b4d0-27ce-11eb-330b-2dfe47120ede
# ╠═f5a06d50-27ce-11eb-17b4-c1599c1d93db
# ╠═20b3f42e-27cf-11eb-00e7-73c83aa66e00
# ╠═4bf84f10-27cf-11eb-0435-a12463478ae0
# ╠═15da21f0-27d0-11eb-025b-e7f9374b020b
# ╠═a1d0c69e-27da-11eb-163b-ab92a7859f62
# ╠═2f42f450-27d0-11eb-1f41-cf1df27badc2
# ╠═bae39552-27d0-11eb-1c3b-cb0bfcacdb81
# ╠═82a942e0-2869-11eb-17f2-5550a7d223a0
# ╠═867c6aee-2869-11eb-15bb-2b5934eaba9f
# ╠═b8f05f50-2869-11eb-158b-47205e98fd88
# ╠═65a19ee0-2869-11eb-0f9d-6d7759d79d9f
# ╠═789ea780-286a-11eb-28dc-2185259857bd
# ╟─798a4000-286a-11eb-1f12-e30338eb7ead
# ╠═5455e46e-27d3-11eb-1478-0f9879fc2215
# ╠═268df680-27d4-11eb-175c-f79f746fbab7
# ╠═5bc8ac00-27d4-11eb-31aa-edb430062b1b
# ╠═7e894f10-27d4-11eb-2b28-dbed5929b12b
# ╠═96d3e030-27d4-11eb-39ea-db39558701a6
# ╠═469ee2e2-27d4-11eb-3e68-39a1c8de84e3
# ╠═460879d0-27d5-11eb-2495-7920d55f677c
# ╟─15281992-27dc-11eb-2e1e-1b4d5028142b
# ╠═1985ce60-27dc-11eb-3aa6-f55d06695c2d
# ╠═e4bf4200-27dc-11eb-051f-8b50f7f4bdd8
# ╠═ae24ea8e-27de-11eb-216e-55bf2ed07117
# ╟─1cad742e-27dc-11eb-1d62-ed5ef14f01a0
# ╠═4147f9a0-27d7-11eb-32de-a9d04fc45211
# ╠═430dc580-27d7-11eb-2edd-81e1104ab57b
# ╠═de03e060-27d7-11eb-19e4-f303bb9078f3
# ╠═edbb0240-27d7-11eb-1c26-e3a457165160
# ╠═384e63a0-27d9-11eb-1540-7731666f9aac
# ╠═dde1c780-27d9-11eb-2345-b9067ba467b7
# ╠═c0cf6762-27d4-11eb-3a68-e1fc72404bd7
# ╠═bc9ee040-27d8-11eb-2ae2-fd3ad839dffb
# ╟─880aede0-27df-11eb-342c-af78db5ae8bb
# ╠═9207b11e-27df-11eb-262b-83039547db91
# ╠═c06496a0-27df-11eb-0557-35927114712c
# ╠═e9f2b600-27df-11eb-3190-d3a359b8b5cb
# ╠═314cdad0-27e0-11eb-261b-359761d4268b
# ╠═1a62c410-27e0-11eb-14ab-339edb9da195
# ╠═b41442a0-2830-11eb-03d0-ed6c96e9dde9
# ╠═2117c860-287a-11eb-11f2-5304606bcdf3
# ╟─ab2aeb90-287a-11eb-0a90-ff7e62c0f950
# ╠═4b99eee0-2831-11eb-20e4-cbae0c67e4fb
# ╠═ce96ccf0-2831-11eb-3f91-4bbe163284b4
# ╠═f7aada50-2831-11eb-2318-0f64039d365d
# ╠═eb076110-2831-11eb-2151-6ba1beee687c
# ╠═0a565d00-2832-11eb-3110-a9c7d74cb153
# ╠═4499ad50-2832-11eb-1316-a3370fe12816
# ╠═34b11b30-2832-11eb-0cb5-775f817a2d7c
# ╠═15abd26e-2832-11eb-1cae-818f96fce3f9
# ╠═83009a90-287d-11eb-364a-014efae85a3a
# ╠═89381872-287d-11eb-0991-bdd241001a4c
# ╟─16d1a010-287a-11eb-1ef2-eb3e698fe85f
# ╠═33caec30-287a-11eb-099b-5b6571982daf
# ╠═4af4cb10-287a-11eb-3718-f7d59a2c0eed
# ╠═5e6a0660-287a-11eb-0308-874829b878aa
# ╠═678d3280-287a-11eb-2c43-bd919979d395
# ╠═70830d10-287a-11eb-269c-c3fa1eaeb9e2
# ╠═05517e02-287a-11eb-167c-0935bd9a2a70
# ╠═675fb282-287d-11eb-0c1d-7d21c31b3ae6
# ╠═72e42730-287d-11eb-26d4-e9cbe9b7b2cc
# ╟─cdf44130-287a-11eb-13dd-796a7a81dd8a
# ╠═e02a6c70-2830-11eb-08d9-631b8f5c8166
# ╠═9c05b770-27e0-11eb-2116-55c62bb8e032
# ╠═89daa100-2830-11eb-0ab0-09728faf861b
# ╠═9ab0ee80-2830-11eb-15eb-6dcee557eadd
# ╠═a4de85c0-2830-11eb-181a-079c0c560cb6
# ╠═b4aabd40-27e3-11eb-1d99-4bb53010e540
# ╠═b31de312-27e0-11eb-0a4e-cde578822d68
# ╠═4999e350-27e8-11eb-143b-dd46d2a2a8d9
# ╠═df09c890-27e0-11eb-1509-afdf248cd717
# ╠═7a323e12-27e1-11eb-3720-a35f5fff7948
# ╠═90073b4e-27e1-11eb-1b41-9d02c393df7c
# ╠═f18d3460-27e1-11eb-0351-f9b44125b4bc
# ╠═885a17f0-27e2-11eb-2720-1d13935ac861
# ╠═563506e0-27e2-11eb-154d-cb3af1504823
# ╠═502164b0-27e2-11eb-1201-b93ded002d96
# ╠═291ceb80-27e4-11eb-32d3-69542e48e4e7
# ╠═90ae7b00-27e5-11eb-0991-532c3e8e3941
# ╠═f6f1ae00-27e5-11eb-2538-6dcf631a43b3
# ╠═a51cbc50-27e5-11eb-11d0-317e3a30e5ea
# ╠═eddba8a0-27e2-11eb-3c2b-456abfb6d91c
# ╠═8c6f5540-27e6-11eb-226b-3d61418cc764
# ╠═a8bb0e30-27e4-11eb-1578-53806dd0f0e9
# ╠═32323f80-27e5-11eb-1581-3537f89b6549
# ╠═54351b60-27e6-11eb-3898-15504bbff380
# ╠═dcb73e50-27e6-11eb-3f42-35c9f67c96da
# ╠═5c03be90-27e7-11eb-0d07-6529c6633b4f
# ╠═9bd64200-27eb-11eb-2c81-f7e1bb552e84
# ╠═7a474bb2-27ec-11eb-3a1b-4d6f11cc8042
# ╠═6ff44460-27e7-11eb-371e-9dddb1d9a76c
# ╠═8adfbf20-27e7-11eb-0cec-91f326f42e55
# ╠═0ac1d330-27e9-11eb-148e-156e7d670757
# ╠═734e20f0-27eb-11eb-11fb-d51498c1592f
# ╠═c1f48f00-27eb-11eb-13af-599491b28241
# ╠═0bd585c0-27ec-11eb-32ea-f33ef2084890
# ╠═81a56f40-27ec-11eb-1bf8-b30e0f259dc4
# ╠═2e574b10-27ec-11eb-26ea-db2bfc5278fe
# ╠═67663a10-27ec-11eb-0173-27217c95d913
# ╠═eb2cc630-27eb-11eb-0496-831311ac09e9
# ╠═8a935a92-27ec-11eb-3508-d9ca334ea4b5
# ╠═96b791b0-27ec-11eb-27a4-bf55cb102ee2
# ╠═292d3070-27ea-11eb-088b-efc17f03f787
# ╠═273869e0-282d-11eb-12ae-e90ddca72641
# ╠═47fee610-27ea-11eb-0a2c-c1b0c93f4446
# ╠═29974c00-282e-11eb-044e-a972ed21be77
# ╠═0e986a50-2870-11eb-2bdc-ef8bd3938307
# ╠═5936b910-282d-11eb-107f-69a8742adfa0
# ╠═632eb3a0-2828-11eb-1933-7b74bf88bcdc
# ╠═3584b5b0-2870-11eb-3d47-7762282618a9
# ╠═ab1efe50-2872-11eb-0a14-f900c5907e05


text:
function is_valid(x, ::Val{:domain_wall})
    for (i, b) in enumerate(x)
        iszero(b) && (return iszero(x[i+1:end]))
    end
    return true
end

function binarize(x, d::D, ::Val{:domain_wall}) where {T <: Number, D <: DiscreteDomain{T}}
    b0 = typeof(d) <: RangeDomain && first(get_domain(d)) == 0
    y = zeros(T, length(d) - b0)
    is_in = false
    for (i, v) in enumerate(get_domain(d))
        i > b0 && (y[i - b0] = 1)
        (is_in = x == v) && break
    end
    @assert is_in "The value $x is not in the domain $d"
    return y
end

function debinarize(x, d::D, ::Val{:domain_wall}) where {T <: Number, D <: DiscreteDomain{T}}
    b0::Bool = typeof(d) <: RangeDomain && first(get_domain(d)) == 0
    at_one = true
    val = typemax(T)
    for (i, v) in zip((b0 ? 0 : 1):length(x), get_domain(d))
        b = i > 0 ? x[i] : 1
        if at_one
            b == 1 ? (val = v) : (at_one = false)
        else
            b == 1 && return typemax(T)
        end
    end
    return val
end

function debinarize(x, domains::Vector{D}, ::Val{:domain_wall}
) where {T <: Number, D <: DiscreteDomain{T}}
    start = 0
    stop = 0
    function aux(d)
        b0 = typeof(d) <: RangeDomain && first(get_domain(d)) == 0
        start = stop + 1
        stop += length(d) - b0

        return debinarize(x[start:stop], d, Val(:domain_wall))
    end
    return map(aux, domains)
end