text:
export get_cstrname_lower, get_cstrname_upper, get_cstrname_eq
export get_momentcstrname, get_blockname, get_cstrname
export get_normalizedpoly, get_pbcstrname, get_ccmultvar, format_string, shortname_moment
export change_eq_to_ineq!

get_cstrname_lower(cstrname::String) = cstrname*"_lo"
get_cstrname_upper(cstrname::String) = cstrname*"_hi"
get_cstrname_eq(cstrname::String) = cstrname*"_eq"

get_momentcstrname() = "moment_cstr"

function get_blockname(cstrname, cliquename, mmtrelax_pb)
    cstrcliques = Set(Iterators.filter(x->x[1] == cstrname, keys(mmtrelax_pb.constraints)))
    if length(cstrcliques) == 1
        return cstrname
    else
        cstrname == get_momentcstrname() || warn(LOGGER, "get_blockname(): several cliques found for $cstrname : $cstrcliques")
        return cstrname*"_"*cliquename
    end
end

function get_cstrname(cstrname::String, cstrtype::Symbol)
    if cstrtype == :ineqdouble
        return get_cstrname_lower(cstrname), get_cstrname_upper(cstrname)
    elseif cstrtype == :eq
        return get_cstrname_eq(cstrname)
    elseif cstrtype == :ineqlo
        return cstrname
        # return get_cstrname_lower(cstrname)
    elseif cstrtype == :ineqhi
        return cstrname
        # return get_cstrname_upper(cstrname)
    else
        error(LOGGER, "get_cstrname(): Unknown type $cstrtype.")
    end
end

function get_normalizedpoly(cstr::Constraint, cstrtype::Symbol)
    if cstrtype == :eq
        return cstr.p - cstr.lb
    elseif cstrtype == :ineqlo
        return cstr.p - cstr.lb
    elseif cstrtype == :ineqhi
        return cstr.ub - cstr.p
    else
        error(LOGGER, "get_normalizedpoly(): Unhandeld type $cstrtype.")
    end
end

function get_pbcstrname(cstrname::String)
    if ismatch(r".+(\_lo|\_hi|\_eq)", cstrname)
        return cstrname[1:end-3]
    else
        return cstrname
    end
end

function get_ccmultvar(relaxctx::RelaxationContext, moment::Exponent, clique1::String, clique2::String)
    if relaxctx.hierarchykind == :Real
        return Variable("lagmult_cc_$(format_string(moment))_$(clique1)_$(clique2)", Real)
    elseif relaxctx.hierarchykind == :Complex
        return Variable("lagmult_cc_$(format_string(moment))_$(clique1)_$(clique2)", Complex)
    else
        error(LOGGER, "get_ccmultvar(): Unhandled hierarchy kind $(relaxctx.hierarchykind) ")
    end
end

function format_string(α::Exponent, β::Exponent)
    s = "$α,$β"
    return replace(s, " " => "")
end

function format_string(s1::String, s2::String)
    s = "$s1,$s2"
    return replace(s, " " => "")
end

function format_string(α::Exponent)
    s = "$α"
    return replace(s, " " => "")
end

function format_string(α::Exponent, block::String)
    s = "$(α)_$block"
    return format_string(α)*"_$block"
end

shortname_moment(n::Int) = "mmt_$n"

function change_eq_to_ineq!(problem::Problem)
    for (ctrname, ctr) in problem.constraints
        if get_cstrtype(ctr) == :eq
            rm_constraint!(problem, ctrname)
            add_constraint!(problem, get_cstrname_lower(ctrname), 0 << (ctr.p - ctr.lb))
            add_constraint!(problem, get_cstrname_upper(ctrname), 0 << (ctr.ub - ctr.p))
        end
    end
end

## Misc
function print_cmat(mat::AbstractArray, round = 1e-3)
    for i=1:size(mat, 1)
        for j=1:size(mat, 2)
            re, im = real(mat[i, j]), imag(mat[i, j])
            @printf("% 5.4f", re)
            @printf(" ")
            @printf("%+5.4fim", im)
            @printf("   ")
        end
        @printf("\n")
    end
end

import Base: hashindex, isslotempty, isslotmissing, isslotfilled, _setindex!, rehash!

### Efficient dict manipulation
"""
    addindex!(h::Dict{K,V}, v0, key::K) where V<:Number where K

    Add the value `v0` to `h[key]` if `key` is already a key, else add the pair `key=>v0`.
"""
function addindex!(h::Dict{K,V}, v0, key::K) where V<:Number where K
    v = convert(V, v0)
    index = ht_keyindex2!(h, key)

    if index > 0
        h.age += 1
        @inbounds h.keys[index] = key
        @inbounds h.vals[index] += v
    else
        @inbounds _setindex!(h, v, key, -index)
    end

    return h
end

#####################################################################################
###  Utils functions, duplicate from base/dict.jl ... (TODO: fix this ducplicate)
#####################################################################################

# get the index where a key is stored, or -pos if not present
# and the key would be inserted at pos
# This version is for use by setindex! and get!
function ht_keyindex2!(h::Dict{K,V}, key) where V where K
    maxallowedprobe = 16
    maxprobeshift   = 6

    age0 = h.age
    sz = length(h.keys)
    iter = 0
    maxprobe = h.maxprobe
    index = hashindex(key, sz)
    avail = 0
    keys = h.keys

    @inbounds while true
        if isslotempty(h,index)
            if avail < 0
                return avail
            end
            return -index
        end

        if isslotmissing(h,index)
            if avail == 0
                # found an available slot, but need to keep scanning
                # in case "key" already exists in a later collided slot.
                avail = -index
            end
        elseif key === keys[index] || isequal(key, keys[index])
            return index
        end

        index = (index & (sz-1)) + 1
        iter += 1
        iter > maxprobe && break
    end

    avail < 0 && return avail

    maxallowed = max(maxallowedprobe, sz>>maxprobeshift)
    # Check if key is not present, may need to keep searching to find slot
    @inbounds while iter < maxallowed
        if !isslotfilled(h,index)
            h.maxprobe = iter
            return -index
        end
        index = (index & (sz-1)) + 1
        iter += 1
    end

    rehash!(h, h.count > 64000 ? sz*2 : sz*4)

    return ht_keyindex2!(h, key)
end


text:
using Base.Test

include("collatz-conjecture.jl")

# canonical data
@test collatz_steps(1) == 0
@test collatz_steps(16) == 4
@test collatz_steps(12) == 9
@test collatz_steps(1000000) == 152
@test_throws DomainError collatz_steps(0)
@test_throws DomainError collatz_steps(-15)

ext: jl


text:
export vandpoly
Reduce.Preload() # required otherwise Reduce has a fit

# The Vandermonde matrix method

macro genvandpoly(index,leny,x,y) # generate poly functions using the Vandermonde matrix method
    leny = :($leny)
    symX = []
    symY = []
    A = Array{Expr}(undef, leny, leny)
    for i = 1:leny
        push!(symX, Meta.parse("x$i"))
        push!(symY, Meta.parse("y$i"))
    end
    for i = 1:leny # Assemble A matrix
        for j = 1:leny
            pwr = leny-j
            var = symX[i]
            A[i,j] = :($var^$pwr)
        end
    end
    coefs = Algebra.:*(Algebra.inv(A),symY) # invert to find coefs of polynomial
    for i = 1:leny
        str = repr(coefs[i])
        for j = 1:leny
            str = replace(str,"y"*string(j)=>"$y"*"["*string(j)*"]")
            str = replace(str,"x"*string(j)=>"$x"*"["*string(j)*"]")
        end
        coefs[i] = Meta.parse(chop(str, head=2, tail=1))
    end

    return coefs[:($index)]
end

function vandpoly(x,y)
    xLength = length(x)
    if xLength == 2
        c1 = @genvandpoly(1, 2,:($x),:($y))
        c2 = @genvandpoly(2, 2,:($x),:($y))
        return Polynomial([c2,c1])
    elseif xLength == 3
        c1 = @genvandpoly(1, 3,:($x),:($y))
        c2 = @genvandpoly(2, 3,:($x),:($y))
        c3 = @genvandpoly(3, 3,:($x),:($y))
        return Polynomial([c3,c2,c1])
    elseif xLength == 4
        c1 = @genvandpoly(1, 4,:($x),:($y))
        c2 = @genvandpoly(2, 4,:($x),:($y))
        c3 = @genvandpoly(3, 4,:($x),:($y))
        c4 = @genvandpoly(4, 4,:($x),:($y))
        return Polynomial([c4,c3,c2,c1])
    elseif xLength == 5
        c1 = @genvandpoly(1, 5,:($x),:($y))
        c2 = @genvandpoly(2, 5,:($x),:($y))
        c3 = @genvandpoly(3, 5,:($x),:($y))
        c4 = @genvandpoly(4, 5,:($x),:($y))
        c5 = @genvandpoly(5, 5,:($x),:($y))
        return Polynomial([c5,c4,c3,c2,c1])
    else
        Polynomial(inv([xval^j for xval in eachindex(x), j = 0:xlength - 1]) * y) # create poly coefficient by inverting array of linear equations
        println("Warning, function behaves slow when exceeding more than 5 datapoints")
    end
end


text:
const libGR3_syms = [
:gr3_geterror
:gr3_init
:gr3_free
:gr3_terminate
:gr3_useframebuffer
:gr3_usecurrentframebuffer
:gr3_getimage
:gr3_export
:gr3_getrenderpathstring
:gr3_drawimage
:gr3_createmesh
:gr3_createindexedmesh
:gr3_drawmesh
:gr3_createheightmapmesh
:gr3_drawheightmap
:gr3_deletemesh
:gr3_setquality
:gr3_clear
:gr3_cameralookat
:gr3_setcameraprojectionparameters
:gr3_setlightdirection
:gr3_drawcylindermesh
:gr3_drawconemesh
:gr3_drawspheremesh
:gr3_drawcubemesh
:gr3_setbackgroundcolor
:gr3_createisosurfacemesh
:gr3_surface
:gr_volume
:gr3_createxslicemesh
:gr3_createyslicemesh
:gr3_createzslicemesh
]


text:
function testDCT1d(N=32)
  Random.seed!(1235)
  x = zeros(ComplexF64, N^2)
  for i=1:5
    x .+= rand()*cos.(rand(1:N^2)*collect(1:N^2)) .+ 1im*rand()*cos.(rand(1:N^2)*collect(1:N^2))
  end
  D1 = DCTOp(ComplexF64,(N^2,),2)
  D2 = sqrt(2/N^2)*[cos(pi/(N^2)*j*(k+0.5)) for j=0:N^2-1,k=0:N^2-1]
  D2[1,:] .*= 1/sqrt(2)
  D3 = DCTOp(ComplexF64,(N^2,),4)
  D4 = sqrt(2/N^2)*[cos(pi/(N^2)*(j+0.5)*(k+0.5)) for j=0:N^2-1,k=0:N^2-1]

  y1 = D1*x
  y2 = D2*x

  @test norm(y1 - y2) / norm(y1) ≈ 0 atol=0.01

  y3 = D3*x
  y4 = D4*x
  @test norm(y3 - y4) / norm(y1) ≈ 0 atol=0.01

  x1 = adjoint(D1)*y1
  @test norm(x1 - x) / norm(x) ≈ 0 atol=0.01

  x2 = adjoint(D3)*y3
  @test norm(x2 - x) / norm(x) ≈ 0 atol=0.01
end

function testFFT1d(N=32,shift=true)
  Random.seed!(1234)
  x = zeros(N^2)
  for i=1:5
    x .+= rand()*cos.(rand(1:N^2)*collect(1:N^2))
  end
  D1 = FFTOp(ComplexF64,(N^2,),shift)
  D2 =  1.0/N*[exp(-2*pi*im*j*k/N^2) for j=0:N^2-1,k=0:N^2-1]

  y1 = D1*x
  if shift
    y2 = fftshift(D2*fftshift(x))
  else
    y2 = D2*x
  end
  @test norm(y1 - y2) / norm(y1) ≈ 0 atol=0.01

  x1 = adjoint(D1) * y1
  if shift
    x2 = ifftshift(adjoint(D2)*ifftshift(y2))
  else
    x2 = adjoint(D2)*y2
  end
  @test norm(x1 - x2) / norm(x1) ≈ 0 atol=0.01
end

function testFFT2d(N=32,shift=true)
  Random.seed!(1234)
  x = zeros(N^2)
  for i=1:5
    x .+= rand()*cos.(rand(1:N^2)*collect(1:N^2))
  end
  D1 = FFTOp(ComplexF64,(N,N),shift)

  idx = CartesianIndices((N,N))[collect(1:N^2)]
  D2 =  1.0/N*[ exp(-2*pi*im*((idx[j][1]-1)*(idx[k][1]-1)+(idx[j][2]-1)*(idx[k][2]-1))/N) for j=1:N^2, k=1:N^2 ]

  y1 = D1*x
  if shift
    y2 = D2*vec(fftshift(reshape(x,N,N)))
    y2 = vec(fftshift(reshape(y2,N,N)))
  else
    y2 = D2*x
  end
  @test norm(y1 - y2) / norm(y1) ≈ 0 atol=0.01

  x1 = adjoint(D1) * y1
  if shift
    x2 = adjoint(D2)*vec(ifftshift(reshape(y2,N,N)))
    x2 = vec(ifftshift(reshape(x2,N,N)))
  else
    x2 = adjoint(D2)*y2
  end
  @test norm(x1 - x2) / norm(x1) ≈ 0 atol=0.01
end

function testWeighting(N=512)
  Random.seed!(1234)
  x1 = rand(N)
  weights = rand(N)
  W = WeightingOp(weights)
  y1 = W*x1
  y = weights .* x1

  @test norm(y1 - y) / norm(y) ≈ 0 atol=0.01

  x2 = rand(2*N)
  W2 = WeightingOp(weights,2)
  y2 = W2*x2
  y = repeat(weights,2) .* x2

  @test norm(y2 - y) / norm(y) ≈ 0 atol=0.01
end

function testGradOp1d(N=512)
  x = rand(N)
  G = GradientOp(eltype(x),size(x))
  G0 = Bidiagonal(ones(N),-ones(N-1), :U)[1:N-1,:]

  y = G*x
  y0 = G0*x
  @test norm(y - y0) / norm(y0) ≈ 0 atol=0.001

  xr = transpose(G)*y
  xr0 = transpose(G0)*y0

  @test norm(y - y0) / norm(y0) ≈ 0 atol=0.001
end

function testGradOp2d(N=64)
  x = repeat(1:N,1,N)
  G = GradientOp(eltype(x),size(x))
  G_1d = Bidiagonal(ones(N),-ones(N-1), :U)[1:N-1,:]

  y = G*vec(x)
  y0 = vcat( vec(G_1d*x), vec(x*transpose(G_1d)) )
  @test norm(y - y0) / norm(y0) ≈ 0 atol=0.001

  xr = transpose(G)*y
  y0_x = reshape(y0[1:N*(N-1)],N-1,N)
  y0_y = reshape(y0[N*(N-1)+1:end],N,N-1)
  xr0 = transpose(G_1d)*y0_x + y0_y*G_1d
  xr0 = vec(xr0)

  @test norm(xr - xr0) / norm(xr0) ≈ 0 atol=0.001
end

function testSampling(N=64)
  x = rand(ComplexF64,N,N)
  idx = shuffle(collect(1:N^2)[1:N*div(N,2)])
  SOp = SamplingOp(idx,(N,N))
  y = SOp*vec(x)
  x2 = adjoint(SOp)*y
  # references
  y_ref = vec(x[idx])
  x2_ref = zeros(ComplexF64,N^2)
  x2_ref[idx] .= y_ref
  # perform tests
  @test norm(y - y_ref) / norm(y_ref) ≈ 0 atol=0.000001
  @test norm(x2 - x2_ref) / norm(x2_ref) ≈ 0 atol=0.000001
end

function testWavelet(M=64,N=60)
  x = rand(M,N)
  WOp = WaveletOp(Float64,(M,N))
  x_wavelet = WOp*vec(x)
  x_reco = reshape( adjoint(WOp)*x_wavelet, M, N)

  @test norm(x_reco - x) / norm(x) ≈ 0 atol=0.001
end

@testset "Linear Operators" begin
  @info "test DCT-II and DCT-IV"
  for N in [2,8,16,32]
    testDCT1d(N)
  end
  @info "test FFT"
  for N in [8,16,32]
    testFFT1d(N,false)
    testFFT1d(N,true)
    testFFT2d(N,false)
    testFFT2d(N,true)
  end
  @info "test Weighting"
  testWeighting(512)
  @info "test gradientOp"
  testGradOp1d(512)
  testGradOp2d(64)
  @info "test sampling"
  testSampling(64)
  @info "test WaveletOp"
  testWavelet(64,64)
  testWavelet(64,60)
end


text:
const IMGUI_DEMO_DRAG_INT_STATES = Dict(
    "basic_drag_int" => DragInts.State(
        "drag int",
        50,
        speed=1,
    ),
    "basic_drag_int2" => DragInts.State(
        "drag int 0..100",
        42,
        speed=1,
        range=(0, 100),
        format="%d%%",
    ),
)

const IMGUI_DEMO_DRAG_FLOAT_STATES = Dict(
    "basic_drag_float" => DragFloats.State(
        "drag float",
        1.0f0,
        speed=0.005,
    ),
    "basic_drag_small_float" => DragFloats.State(
        "drag small float",
        0.0067f0,
        speed=0.0001,
        range=(0.0f0, 0.0f0),
        format="%.06f ns",
    ),
)


text:
#######################
# Linearity detection #
#######################

export dim, detecthlinearity!, detecthlinearity, detectvlinearity!, detectvlinearity, detect_new_linearities

"""
    dim(h::HRep, current=false)

Returns the dimension of the affine hull of the polyhedron.
That is the number of non-redundant hyperplanes that define it.
If `current` is `true` then it simply returns the dimension according the current number of hyperplanes, assuming that the H-linearity has already been detected.
Otherwise, it first calls [`detecthlinearity!`](@ref).
"""
function dim(h::HRep, current=false)
    if !current
        detecthlinearity!(h)
    end
    fulldim(h) - nhyperplanes(h)
end

"""
    detecthlinearity!(p::VRep)

Detects all the hyperplanes contained in the H-representation and remove all redundant hyperplanes.

## Examples
The representation
```julia
h = HalfSpace([1, 1], 1]) ∩ HalfSpace([-1, -1], -1)
```
contains the hyperplane `HyperPlane([1, 1], 1)`.
"""
detecthlinearity!(p::HRep) = error("detecthlinearity! not implemented for $(typeof(p))")
detecthlinearity!(p::Polyhedron, solver=default_solver(p); verbose=0) = sethrep!(p, detecthlinearity(hrep(p), solver; verbose=verbose))

"""
    detectvlinearity!(p::VRep)

Detects all the lines contained in the V-representation and remove all redundant lines.

## Examples
The representation
```julia
v = conichull([1, 1], [-1, -1])
```
contains the line `Line([1, 1])`.
"""
detectvlinearity!(p::VRep) = error("detectvlinearity! not implemented for $(typeof(p))")
detectvlinearity!(p::Polyhedron, solver=default_solver(p); verbose=0) = setvrep!(p, detectvlinearity(vrep(p), solver; verbose=verbose))

# No ray so no line
function detectvlinearity!(p::VPolytope) end
detectvlinearity(p::VPolytope) = p

linearize(h::HalfSpace) = hyperplane(h)
linearize(r::Ray) = line(r)
_aff_push!(aff::HAffineSpace, h) = intersect!(aff, h)
_aff_push!(aff::VLinearSpace, l) = convexhull!(aff, l)

function _detect_opposite_element(aff, non_opposite, element)
    element = remproj(element, aff)
    if !isapproxzero(element) && !any(el -> el ≈ element, non_opposite)
        lin = linearize(element)
        i = findfirst(el -> linearize(el) ≈ lin, non_opposite)
        if i === nothing
            push!(non_opposite, element)
        else
            deleteat!(non_opposite, i)
            _aff_push!(aff, lin)
        end
        return i !== nothing
    else
        return false
    end
end
function _detect_opposite_elements(aff, non_opposite, elements)
    newlin = true
    for i in 1:fulldim(aff) # could use `while newlin` but `for`-loop is safer.
        newlin || break
        newlin = false
        empty!(non_opposite)
        # Project each ray/halfspace orthogonal to the line hull/affine hull.
        # Remove rays/halfspaces that become zero and detect new lines/hyperplanes
        # with rays/halfspaces that becomes opposite to each other.
        for element in elements
            newlin |= _detect_opposite_element(aff, non_opposite, element)
        end
    end
end

"""
    detect_new_linearities(rep::HRepresentation, solver; verbose=0)

Given a polyhedron with H-representation `rep`, detect whether a new hyperplane can be generated from the halfspaces in `halfspaces` using an linear program solved by `solver`.
The method is similar to the method used for lines described as follows.
This function is automatically called by `removehredundancy` if a solver is provided.

    detect_new_linearities(rep::VRepresentation, solver; verbose=0)

Given a cone defined by the V-representation `rep` (ignoring the points in the representation if any), detect whether a new line can be generated from the rays in `rays` using an linear program solved by `solver`.
The method is as follows (suppose `lines` is empty for simplicity).
This function is automatically called by `removevredundancy` if a solver is provided.

If there was a line `l` in the cone, it would mean that there exist `μ >= 0` and `ν >= 0` such that
`Σ μ_i r_i = l` and `Σ ν_i r_i = -l`. We deduce from this that `Σ λ_i r_i = 0` where `λ = μ + ν`.

Conversely, if there are `λ >= 0` such that `Σ λ_i r_i = 0` then let `j` be the index of `λ` with largest magnitude
(to make sure it is nonzero).
We have `Σ_{i != j} λ_i/λ_j r_i = -r_j`. As both `r_j` and `-r_j` are in the cone,
`r_j` generates a line in the cone.
However, this means that we now have `Σ_{i != j} λ_i/λ_j r_i ≡ 0 (mod r_j)` so if there is another `λ_i`
with nonzero value, we can transform it to a line as well.
In summary, we have a line `r_i` for each `i` such that `λ_i != 0`.

The dual program is:
max z
    r_i'x >= z
When the primal is feasible, the dual program may still be feasible.
We know that `z = 0` by strong duality as the objective value needs to be equal to the objective of the primal which is zero.
So the constraints are `r_i'x >= 0`. If we have `r_i'x > 0` for some `i`, it means that `-r_i` does not belong to the cone
hence `r_i` can be dropped for the purpose of searching for lines.

## Note

In CDDLib, the dual program is solved, if the objective value is zero then linearity are found
by, for each `i` such that `r_i'x = 0`, solve an LP to find whether `-r_i` belongs to the cone.
CDDLib ignores the primal results provided in `λ` which directly gives linearity without the need
to solve an LP for each ray.
This method is therefore significantly more efficient as it's complexity is `O(dimension of linespace)` which is upper
bounded by `O(fulldim)` while the method of CDDLib is `O(number of rays)`.
"""
function detect_new_linearities(rep::Representation, solver; verbose=0)
    lins = _linearity(rep)
    nonlins = _nonlinearity(rep)
    isempty(nonlins) && return Int[]
    model, T = layered_optimizer(solver)
    verbose >= 2 && (_model = JuMP.direct_model(model))
    is_lin = falses(length(nonlins))
    active = trues(length(nonlins))
    # We pass `true` as we break homogeneity of the problem with `sum(λ) == 1`.
    hull = _zero_hull(rep, T)
    λ, cλ, sum_con = _hull(model, T, hull, rep, eachindex(nonlins), true)
    if verbose >= 2
        for (i, λ) in enumerate(λ)
            MOI.set(model, MOI.VariableName(), λ, "λ[$i]")
        end
    end
    if !isempty(lins)
        η, _, _ = _hull(model, T, hull, rep, eachindex(lins))
    end
    if rep isa HRepresentation
        MOI.add_constraint.(model, hull[1:end-1], MOI.EqualTo(zero(T)))
        β_con = MOI.add_constraint(model, hull[end], MOI.LessThan(zero(T)))
    else
        MOI.add_constraint.(model, hull, MOI.EqualTo(zero(T)))
    end
    for i in 1:fulldim(rep) # safer than `while ...`
        (count(is_lin) == length(is_lin) || iszero(count(active))) && break
        verbose >= 2 && println(_model)
        # FIXME stopping when we have enough hyperplanes to prove that it's empty
        #       should also resolve the issue with presolve.
        is_feasible(model, "detecting new linearity (you may need to activate presolve for some solvers, e.g. by replacing `GLPK.Optimizer` with `optimizer_with_attributes(GLPK.Optimizer, \"presolve\" => GLPK.GLP_ON)`).") || break
        if rep isa HRepresentation && !isapproxzero((β_primal = MOI.get(model, MOI.ConstraintPrimal(), β_con);))
            verbose >= 1 && @info("The polyhedron is empty as $β_primal is negative.")
            return nothing
        end
        ray_to_line = Int[]
        ray_to_drop = Int[]
        for i in eachindex(λ)
            if active[i] && !is_lin[i]
                if !isapproxzero((primal = MOI.get(model, MOI.VariablePrimal(), λ[i]);))
                    # `λ_i > 0`, we know that `-r_i` belongs to the cone so we transform the ray into a line.
                    verbose >= 1 && @info("$(i)th element is linear as $primal is positive.")
                    is_lin[i] = true
                    push!(ray_to_line, i)
                elseif !isapproxzero((dual = MOI.get(model, MOI.ConstraintDual(), cλ[i]);))
                    # `r_i'x > 0`, we know that `r_i` does not belong to the cone so we just drop the ray from the search for lines.
                    verbose >= 1 && @info("$(i)th element is nonlinear as $dual is positive.")
                    active[i] = false
                    push!(ray_to_drop, i)
                end # otherwise, we are still uncertain about this ray and we'll do a new solve.
            end
        end
        # Query `primal` and `dual` before doing any `MOI.delete` and `MOI.modify` as they
        # won't be available afterwards.
        for i in ray_to_line
            MOI.delete(model, cλ[i])
            MOI.modify(model, sum_con, MOI.ScalarCoefficientChange(λ[i], 0.0))
        end
        for i in ray_to_drop
            MOI.delete(model, λ[i])
        end
    end
    return findall(is_lin)
end

_no_nonlin_word(::HRepresentation) = "affine"
_no_nonlin_word(::VRepresentation) = "bounded"
_lin_word(::HRepresentation) = "halfspaces"
_lin_word(::VRepresentation) = "lines"
_nonlin_word(::HRepresentation) = "hyperplanes"
_nonlin_word(::VRepresentation) = "rays"
_nonlin_type(h::HRepresentation) = halfspacetype(h)
_nonlin_type(v::VRepresentation) = raytype(v)
_hasnonlinearity(h::HRepresentation) = hashalfspaces(h)
_hasnonlinearity(v::VRepresentation) = hasrays(v)
_nonlinearity(h::HRepresentation) = halfspaces(h)
_nonlinearity(v::VRepresentation) = rays(v)
_linearity_space(h::HRepresentation, current) = affinehull(h, current)
_linearity_space(v::VRepresentation, current) = linespace(v, current)

struct OppositeMockOptimizer end
function _detect_linearity(rep::Representation, solver; verbose=0)
    aff = _linearity_space(rep, true)
    if _hasnonlinearity(rep)
        if solver === nothing
            @warn("""
Cannot detect exact linearity as no solver was provided and the polyhedron is not $(_no_nonlin_word(rep)).
As fallback, we will only detect $(_lin_word(rep)) from opposite $(_nonlin_word(rep)) but that may not detect all $(_lin_word(rep)).
Set a solver if you believe that the polyhedron may have more linearity.
""" * NO_SOLVER_HELP)
            solver = OppositeMockOptimizer
        end
        if solver == OppositeMockOptimizer
            els = _nonlin_type(rep)[]
            _detect_opposite_elements(aff, els, _nonlinearity(rep))
        else
            new_lins = detect_new_linearities(rep, solver; verbose=verbose)
            if new_lins === nothing
                empty!(aff)
                els = eltype(_nonlinearity(rep))[]
            else
                els = _nonlinearity(rep)
                if !isempty(new_lins)
                    for idx in collect(eachindex(els))[new_lins]
                        _aff_push!(aff, linearize(get(rep, idx)))
                    end
                end
                els = elements_without(els, new_lins)
            end
        end
    else
        els = _nonlinearity(rep)
    end
    return removeduplicates(aff), els
end
function detecthlinearity(hr::HRepresentation, solver; verbose=0)
    aff, hs = _detect_linearity(hr, solver; verbose=verbose)
    typeof(hr)(FullDim(hr), aff.hyperplanes, hs)
end
function detectvlinearity(vr::VRepresentation, solver; verbose=0)
    aff, rays = _detect_linearity(vr, solver; verbose=verbose)
    typeof(vr)(FullDim(vr), preps(vr)..., aff.lines, rays)
end


text:
###############
###############
## partition ##
###############
###############

abstract type AbstractSign end
struct Plus <: AbstractSign end
struct Minus <: AbstractSign end
struct Zero <: AbstractSign end
struct NotPlus <: AbstractSign end
struct NotMinus <: AbstractSign end
struct NotZero <: AbstractSign end

function _sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{Plus},i::Int)
    if DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
        if ffq.rates.lwr[bndry_idx]>0.0
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{Minus},i::Int)
    if DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
        if ffq.rates.lwr[bndry_idx]<0.0
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{Zero},i::Int)
    if DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
        if ffq.rates.lwr[bndry_idx]==0.0
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{NotPlus},i::Int)
    if DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
        if !(ffq.rates.lwr[bndry_idx]>0.0)
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{NotMinus},i::Int)
    if DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
        if !(ffq.rates.lwr[bndry_idx]<0.0)
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{NotZero},i::Int)
    if DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
        if !(ffq.rates.lwr[bndry_idx]==0.0)
            return true
        else
            return false
        end
    else
        return false
    end
end

function _sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{Plus},i::Int)
    if DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
        if ffq.rates.upr[bndry_idx]>0.0
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{Minus},i::Int)
    if DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
        if ffq.rates.upr[bndry_idx]<0.0
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{Zero},i::Int)
    if DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
        if ffq.rates.upr[bndry_idx]==0.0
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{NotPlus},i::Int)
    if DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
        if !(ffq.rates.upr[bndry_idx]>0.0)
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{NotMinus},i::Int)
    if DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
        if !(ffq.rates.upr[bndry_idx]<0.0)
            return true
        else
            return false
        end
    else
        return false
    end
end
function _sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{NotZero},i::Int)
    if DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S,i)
        bndry_idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
        if !(ffq.rates.upr[bndry_idx]==0.0)
            return true
        else
            return false
        end
    else
        return false
    end
end

_sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Type{<:AbstractSign}) =
    [_sign_map_lwr_bndry(ffq,k,i) for i in 1:n_phases(ffq.dq)]
_sign_map_upr_bndry(ffq::FluidFluidQueue,k::Type{<:AbstractSign}) =
    [_sign_map_upr_bndry(ffq,k,i) for i in 1:n_phases(ffq.dq)]
_sign_map_lwr_bndry(ffq::FluidFluidQueue,k::Colon) =
    [_sign_map_lwr_bndry(ffq,NotPlus,i)||_sign_map_lwr_bndry(ffq,Plus,i) for i in 1:n_phases(ffq.dq)]
_sign_map_upr_bndry(ffq::FluidFluidQueue,k::Colon) =
    [_sign_map_upr_bndry(ffq,NotPlus,i)||_sign_map_upr_bndry(ffq,Plus,i) for i in 1:n_phases(ffq.dq)]

_sign_map_interior(ffq::FluidFluidQueue,k::Type{Plus}) = map(x->x>0.0,ffq.rates.interior)
_sign_map_interior(ffq::FluidFluidQueue,k::Type{Minus}) = map(x->x<0.0,ffq.rates.interior)
_sign_map_interior(ffq::FluidFluidQueue,k::Type{Zero}) = map(x->x==0.0,ffq.rates.interior)
_sign_map_interior(ffq::FluidFluidQueue,k::Type{NotPlus}) = map(x->!(x>0.0),ffq.rates.interior)
_sign_map_interior(ffq::FluidFluidQueue,k::Type{NotMinus}) = map(x->!(x<0.0),ffq.rates.interior)
_sign_map_interior(ffq::FluidFluidQueue,k::Type{NotZero}) = map(x->!(x==0.0),ffq.rates.interior)
_sign_map_interior(ffq::FluidFluidQueue,k::Colon) = trues(size(ffq.rates.interior))

_sign_map(ffq::FluidFluidQueue{T},k) where T<:DiscretisedFluidQueue =
    [_sign_map_lwr_bndry(ffq,k) _sign_map_interior(ffq,k) _sign_map_upr_bndry(ffq,k)]

Base.size(R::FluidFluidOperator{1}) = size(R.array)
Base.getindex(R::FluidFluidOperator{1},i::AbstractVector) = R.array[i]
Base.getindex(R::FluidFluidOperator{1},i::Int) = R.array[i]
# function Base.getindex(R::FluidFluidOperator{1},(ks,is,ps)::Tuple{Any,Any,Any})
#     return Base.getindex(R,(ks,is,ps)...)
# end
function Base.getindex(R::FluidFluidOperator{1},(ks,is,ps)::Tuple{Any,Any,Any})
    rows = DiscretisedFluidQueues._build_index_from_tuple(R.ffq.dq,(ks,is,ps))
    # (typeof(is)==Colon)&&(is=1:n_phases(R.ffq.dq))
    # (typeof(ks)==Colon)&&(ks=0:n_intervals(R.ffq.dq)+1)
    # (typeof(ps)==Colon)&&(ps=1:n_bases_per_cell(R.ffq.dq))

    # rows = Int[]
    # for k in ks
    #     for i in is
    #         for p in ps
    #             if (k==0)&&(p==1)&&DiscretisedFluidQueues._has_left_boundary(R.ffq.dq.model.S,i)
    #                 push!(rows,DiscretisedFluidQueues._map_to_index((0,i,1),R.ffq.dq))
    #             elseif (k==n_intervals(R.ffq.dq)+1)&&(p==1)&&DiscretisedFluidQueues._has_right_boundary(R.ffq.dq.model.S,i)
    #                 push!(rows,DiscretisedFluidQueues._map_to_index((n_intervals(R.ffq.dq)+1,i,1),R.ffq.dq))
    #             elseif (k!=0)&&(k!=n_intervals(R.ffq.dq)+1)
    #                 push!(rows,DiscretisedFluidQueues._map_to_index((k,i,p),R.ffq.dq))
    #             end
    #         end
    #     end
    # end
    return R.array[rows]
end
function _build_index_map_from_tuple_sign(ffq,(k,i,p))
    k_map = findall(_sign_map(ffq,k))
    (i==:)&&(i=1:n_phases(ffq.dq))
    (p==:)&&(p=1:n_bases_per_cell(ffq.dq))
    k_map = k_map[map(x->x[1]∈i,k_map)]
    return k_map
end
function _build_index_from_tuple_sign(ffq,(k,i,p))
    k_map = _build_index_map_from_tuple_sign(ffq,(k,i,p))
    idx = vcat([DiscretisedFluidQueues._build_index_from_tuple(ffq.dq,(k_map[row][2]-1,k_map[row][1],p)) for row in 1:length(k_map)]...)
    return idx
end
function Base.getindex(R::FluidFluidOperator{1},(k,i,p)::Tuple{DataType,Any,Any})
    idx = _build_index_from_tuple_sign(R.ffq,(k,i,p))
    return R[idx]
end
Base.getindex(R::FluidFluidOperator{1},(k,i)::Tuple{Any,Any}) = R[(k,i,:)]
Base.getindex(R::FluidFluidOperator{1},k::Type{<:AbstractSign}) = R[(k,:,:)]

Base.size(B::FluidFluidOperator{2}) = size(B.array)
Base.getindex(B::FluidFluidOperator{2},i,j) = B.array[i,j]

Base.getindex(B::FluidFluidOperator{2},(k,i,p)::Tuple{Any,Any,Any},(l,j,q)::Tuple{Any,Any,Any}) =
    DiscretisedFluidQueues._getindex_from_tuple(B,B.ffq.dq,(k,i,p),(l,j,q))
function _getindex_from_tuple_sign(B::FluidFluidOperator{2},(k,i,p)::Tuple,(l,j,q)::Tuple)
    # k_map = findall(_sign_map(B.ffq,k))
    # l_map = findall(_sign_map(B.ffq,l))
    # (i==:)&&(i=1:n_phases(B.ffq.dq))
    # (j==:)&&(j=1:n_phases(B.ffq.dq))
    # (p==:)&&(p=1:n_bases_per_cell(B.ffq.dq))
    # (q==:)&&(q=1:n_bases_per_cell(B.ffq.dq))
    # k_map = k_map[map(x->x[1]∈i,k_map)]
    # l_map = l_map[map(x->x[1]∈j,l_map)]
    # return vcat([hcat([B[(k_map[row][2]-1,k_map[row][1],p),(l_map[col][2]-1,l_map[col][1],q)] for col in 1:length(l_map)]...) for row in 1:length(k_map)]...)
    rows = _build_index_from_tuple_sign(B.ffq,(k,i,p))
    cols = _build_index_from_tuple_sign(B.ffq,(l,j,q))
    return B[rows,cols]
end
Base.getindex(B::FluidFluidOperator{2},(k,i,p)::Tuple{DataType,Any,Any},(l,j,q)::Tuple{Any,Any,Any}) =
    _getindex_from_tuple_sign(B,(k,i,p),(l,j,q))
Base.getindex(B::FluidFluidOperator{2},(k,i,p)::Tuple{Any,Any,Any},(l,j,q)::Tuple{DataType,Any,Any}) =
    _getindex_from_tuple_sign(B,(k,i,p),(l,j,q))
Base.getindex(B::FluidFluidOperator{2},(k,i,p)::Tuple{DataType,Any,Any},(l,j,q)::Tuple{DataType,Any,Any}) =
    _getindex_from_tuple_sign(B,(k,i,p),(l,j,q))
# function Base.getindex(B::FluidFluidOperator{2},(k,i,p)::Tuple{Any,Any,Any},(l,j,q)::Tuple{Any,Any,Any})
#     k_map = findall(_sign_map(B,k))
#     l_map = findall(_sign_map(B,l))
#     (i==:)&&(i=1:n_phases(B.dq))
#     (j==:)&&(j=1:n_phases(B.dq))
#     (p==:)&&(p=1:n_bases_per_cell(B.dq))
#     (q==:)&&(q=1:n_bases_per_cell(B.dq))
#     k_map = k_map[map(x->x[1]∈i,k_map)]
#     l_map = l_map[map(x->x[1]∈j,l_map)]
#     return vcat([hcat([B[(k_map[row][2]-1,k_map[row][1],p),(l_map[col][2]-1,l_map[col][1],q)] for col in 1:length(l_map)]...) for row in 1:length(k_map)]...)
# end

Base.getindex(B::FluidFluidOperator{2},(k,i)::Tuple{Any,Any},(l,j)::Tuple{Any,Any}) = B[(k,i,:),(l,j,:)]
Base.getindex(B::FluidFluidOperator{2},k::Type{<:AbstractSign},l) = B[(k,:,:),(l,:,:)]
Base.getindex(B::FluidFluidOperator{2},k::Colon,l::Type{<:AbstractSign}) = B[(k,:,:),(l,:,:)]

index(ffq,s::Type{<:AbstractSign}) = _build_index_from_tuple_sign(ffq,(s,:,:))


text:
using Perftests

run_perf_groups()


text:
function test_constructors(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Construct constant ", T, " from scalar, Vector, Matrix") )" begin
        scalar1 = T(1)
        scalar2 = T(1.0)
        V = T([1;2;3])
        M = T([1 2 ; 3 4 ; 5 6])

        @test casadi.is_equal(eval( Meta.parse( string("casadi.", T) ) )(1), scalar1, 1)
        @test casadi.is_equal(eval( Meta.parse( string("casadi.", T) ) )(1.0), scalar2, 1)
        @test casadi.is_equal(scalar1, scalar2, 1)
        @test casadi.is_equal(eval( Meta.parse( string("casadi.", T) ) )([1;2;3]), V, 1)
        @test casadi.is_equal(eval( Meta.parse( string("casadi.", T) ) )([1 2;3 4;5 6]), M, 1)
    end

    @testset "$( string("Construct variable ", T, " vector and matrix          ") )" begin
        v = eval( Meta.parse( string("casadi.", T, ".sym") ) )("v",3)
        V = T("V",3)
        v_val = rand(3)
        @test casadi.is_equal(casadi.substitute(v,v,v_val), casadi.substitute(V,V,v_val), 1)

        m = eval( Meta.parse( string("casadi.", T, ".sym") ) )("m",2,4)
        M = T("M",2,4)
        m_val = rand(2,4)
        @test casadi.is_equal(casadi.substitute(m,m,m_val), casadi.substitute(M,M,m_val), 1)
    end

    @testset "$( string("AbstractFloat and Integer constructors for ", T, "    ") )" begin
        af = [1.2 ; -0.0 ; -4.3]
        @test casadi.is_equal(T(af)', eval( Meta.parse( string("casadi.", T) ) )(af)', 2)

        int = [1 ; 0 ; -2]
        @test casadi.is_equal(T(int)', eval( Meta.parse( string("casadi.", T) ) )(int)', 2)
    end

    @testset "$( string("Other constructors tests for ", T, "                  ") )" begin
        x = T("x")
        @test casadi.is_equal(T(x), x)

        @test convert(T, "y").name() == "y"

        @test to_julia( T(3,2) ) == zeros(3,2)
    end
end


text:
using GeoStatsImages
using GeoStatsBase
using Test

@testset "Basic checks" begin
  # all images loaded as arrays
  for id in GeoStatsImages.available()
    @test geostatsimage(id) isa GeoData
  end

  # throws on non existing identifier
  @test_throws AssertionError geostatsimage("NonExistingIdentifier")
end


text:
module solve_eqm
# Write your package code here.
using NLsolve, LinearAlgebra
include("solve_dist.jl")
include("valueFunctions.jl")
include("helper.jl")
using .solve_dist, .valueFunctions, .helper
function logit(x, fixedParam)
    return 1 ./ (1.0 .+ exp.(-fixedParam.lambda .* x))
end

#Solve for equilibrium
function solve(par, set, fxp, init; purpose="solve")
    kk = set.Ny * set.Na
    ii = set.Nx * set.Ne

    epsf = repeat([-2, 0, 2] .* fxp.sdf .+ par.muf * fxp.mum, inner=7)
    epsm = repeat([-2, 0, 2] .* fxp.sdm .+ fxp.mum, inner=7)
    alpha = repeat(fxp.alpha, inner=7)
    #humanCapital = repeat(fxp.hc, outer = 3);
    #jobProductivityType = repeat(fxp.p,outer=3);
    #jobAmenityType = repeat(fxp.alpha,inner=7);
    output = par.K .* [(par.a * hc^fxp.rho + (1 - par.a) * p^fxp.rho)^(1 / fxp.rho) for hc in repeat(fxp.hc, outer=3), p in repeat(fxp.p, outer=3)]
    outputb = repeat([par.b * hc for hc in fxp.hc], outer=3)
    #This stages solution from nlsolve is slightly different from that from matlab's fsolve. The difference is on the order of 1e-6
    #e = Dict("stages" => nlsolve(x -> solveStages(x, fxp), 12.50 .* ones(typeof(12.50), 8)).zero, "converge" => 1, "Pineg" => 0)
    stages = [0.0714286494397354, 0.117096067390252, 0.0614754501329449, 0.249999769056968, 0.0714293405864623, 0.175877941575804, 0.00269313154484208, 0.250000013739020]
    e = Dict("stages" => stages, "converge" => 1, "Pineg" => 0)
    NC_f = e["stages"][1]
    YC_f = e["stages"][2]
    PL_f = e["stages"][3]
    D_f = e["stages"][4]
    NC_m = e["stages"][5]
    YC_m = e["stages"][6]
    PL_m = e["stages"][7]
    D_m = e["stages"][8]
    gem = par.d1m .+ par.d2m .* fxp.p
    gef = par.d1f .+ par.d2f .* fxp.p

    #unemployed workers u(x), in different stages
    umNC_0 = vec(init.umNC)
    umPL_0 = vec(init.umYC / 2)
    umYC_0 = vec(init.umYC / 2)
    umD_0 = vec(init.umD)

    ufNC_0 = vec(init.ufNC)
    ufPL_0 = vec(init.ufYC / 2)
    ufYC_0 = vec(init.ufYC / 2)
    ufD_0 = vec(init.ufD)
    #employed workers h(x, y, fxp.alpha)
    hmNC_0 = init.hmNC
    hmYC_0 = init.hmYC
    hmPL_0 = init.hmPL
    hmD_0 = init.hmD
    hfNC_0 = init.hfNC
    hfYC_0 = init.hfYC
    hfPL_0 = init.hfPL
    hfD_0 = init.hfD
    #vacancies
    v0 = vec(init.v)
    #Guess Value Functions
    # value in unemployment
    UmNC_0 = vec(init.UmNC)
    UmPL_0 = vec(init.UmYC)
    UmYC_0 = vec(init.UmYC)
    UmD_0 = vec(init.UmD)
    UfNC_0 = vec(init.UfNC)
    UfPL_0 = vec(init.UfYC)
    UfYC_0 = vec(init.UfYC)
    UfD_0 = vec(init.UfD)
    #value in vacancy Π₀(y,α)
    Pi0 = vec(init.Pi)

    #value in surplus
    #Male NC
    SmNC = init.SmNC
    SmNCPos = copy(SmNC)
    SmNCPos[SmNCPos.<0] .= 0
    #Male PL
    SmPL = init.SmPL
    SmPLPos = copy(SmPL)
    SmPLPos[SmPLPos.<0] .= 0
    #Male YC
    SmYC = init.SmYC
    SmYCPos = copy(SmYC)
    SmYCPos[SmYCPos.<0] .= 0
    #Male D
    SmD = init.SmD
    SmDPos = copy(SmD)
    SmDPos[SmDPos.<0] .= 0

    #Female NC
    SfNC = init.SfNC
    SfNCPos = copy(SfNC)
    SfNCPos[SfNCPos.<0] .= 0
    #Female PL
    SfPL = init.SfPL
    SfPLPos = copy(SfPL)
    SfPLPos[SmPLPos.<0] .= 0
    #Female YC
    SfYC = init.SfYC
    SfYCPos = copy(SfYC)
    SfYCPos[SfYCPos.<0] .= 0
    #Female D
    SfD = init.SfD
    SfDPos = copy(SfD)
    SfDPos[SfDPos.<0] .= 0

    #Joint values
    PmNC_0 = SmNC + [x + y for x in UmNC_0, y in Pi0]
    PmPL_0 = SmPL + [x + y for x in UmPL_0, y in Pi0]
    PmYC_0 = SmYC + [x + y for x in UmYC_0, y in Pi0]
    PmD_0 = SmD + [x + y for x in UmD_0, y in Pi0]

    PfNC_0 = SfNC + [x + y for x in UfNC_0, y in Pi0]
    PfPL_0 = SfPL + [x + y for x in UfPL_0, y in Pi0]
    PfYC_0 = SfYC + [x + y for x in UfYC_0, y in Pi0]
    PfD_0 = SfD + [x + y for x in UfD_0, y in Pi0]

    #initial guess for tax rate
    tau = fxp.Transfer * sum(sum(hmPL_0 + hfPL_0, dims=1) .* output, dims=1) / (sum(sum(hmNC_0 + hfNC_0 + hmYC_0 + hfYC_0 + hmD_0 + hfD_0, dims=1) .* output, dims=1))
    #prepare to update these values
    PmNC = zeros(Float64, 21, 21)
    PmPL = zeros(Float64, 21, 21)
    PmYC = zeros(Float64, 21, 21)
    PmD = zeros(Float64, 21, 21)

    PfNC = zeros(Float64, 21, 21)
    PfPL = zeros(Float64, 21, 21)
    PfYC = zeros(Float64, 21, 21)
    PfD = zeros(Float64, 21, 21)

    UmNC = zeros(Float64, 21)
    UmPL = zeros(Float64, 21)
    UmYC = zeros(Float64, 21)
    UmD = zeros(Float64, 21)
    UfNC = zeros(Float64, 21)
    UfPL = zeros(Float64, 21)
    UfYC = zeros(Float64, 21)
    UfD = zeros(Float64, 21)
    Pi = zeros(Float64, 21)

    get!(e, "hmNC", hmNC_0)
    get!(e, "hmPL", hmPL_0)
    get!(e, "hmYC", hmYC_0)
    get!(e, "hmD", hmD_0)
    get!(e, "hfNC", hfNC_0)
    get!(e, "hfPL", hfPL_0)
    get!(e, "hfYC", hfYC_0)
    get!(e, "hfD", hfD_0)

    get!(e, "umNC", umNC_0)
    get!(e, "umPL", umPL_0)
    get!(e, "umYC", umYC_0)
    get!(e, "umD", umD_0)
    get!(e, "ufNC", ufNC_0)
    get!(e, "ufPL", ufPL_0)
    get!(e, "ufYC", ufYC_0)
    get!(e, "ufD", ufD_0)

    get!(e, "v", v0)

    #output surpluses
    get!(e, "SmNC", SmNC)
    get!(e, "SmPL", SmPL)
    get!(e, "SmYC", SmYC)
    get!(e, "SmD", SmD)
    get!(e, "SfNC", SfNC)
    get!(e, "SfPL", SfPL)
    get!(e, "SfYC", SfYC)
    get!(e, "SfD", SfD)

    get!(e, "UmNC", UmNC)
    get!(e, "UmPL", UmPL)
    get!(e, "UmYC", UmYC)
    get!(e, "UmD", UmD)
    get!(e, "UfNC", UfNC)
    get!(e, "UfPL", UfPL)
    get!(e, "UfYC", UfYC)
    get!(e, "UfD", UfD)
    get!(e, "tau", tau)
    get!(e, "Pi", Pi)

    #Value function iterations
    maxit = 10000
    change_dist = zeros(maxit, 1)
    tol = 1e-9
    x0 = vcat(umNC_0, umPL_0, umYC_0, umD_0, ufNC_0, ufPL_0, ufYC_0, ufD_0, vec(hmNC_0), vec(hmYC_0), vec(hmPL_0), vec(hmD_0), vec(hfNC_0), vec(hfYC_0), vec(hfPL_0), vec(hfD_0), vec(v0))
    x = zeros(Float64, length(x0))
    #outer loop
    for itout = 1:maxit

        effU = sum(umNC_0 + ufNC_0 + fxp.s1 * (umYC_0 + umD_0 + ufYC_0 + ufD_0)) +
               fxp.s2 * sum(sum(hmNC_0 + hmYC_0 + hmD_0 + hfNC_0 + hfYC_0 + hfD_0, dims=1))
        V = sum(v0)
        fxp = merge(fxp, (lambdau_NC=fxp.theta / (effU * V)^0.5, lambdau_YC=fxp.theta / (effU * V)^0.5 * fxp.s1, lambdae=fxp.theta / (effU * V)^0.5 * fxp.s2))
        for it1 = 1:maxit
            #Value of Unemployment
            #Male NC
            UmNC = [x + fxp.beta * (fxp.gamma * y + fxp.chi * z +
                                    fxp.lambdau_NC * par.sigma * sum(v0 .* SmNCPos[i, :]) +
                                    (1 - fxp.gamma - fxp.chi) * m)
                    for (x, y, z, m, i) in zip(outputb, UmD_0, UmPL_0, UmNC_0, 1:21)]

            #Male PL
            UmPL = [x + fxp.beta * (fxp.gamma * y + fxp.etam * z +
                                    (1 - fxp.gamma - fxp.etam) * m)
                    for (x, y, z, m) in zip(outputb, UmD_0, UmYC_0, UmPL_0)]

            #Male YC
            UmYC = [x + fxp.beta * (fxp.gamma * y + fxp.chi * z +
                                    fxp.lambdau_YC * par.sigma * sum(v0 .* SmYCPos[i, :]) +
                                    (1 - fxp.gamma - fxp.chi) * m)
                    for (x, y, z, m, i) in zip(outputb, UmD_0, UmPL_0, UmYC_0, 1:21)]

            #Male D
            UmD = [x + fxp.beta * ((1 - fxp.phi) * y +
                                   fxp.lambdau_YC * par.sigma * sum(v0 .* SmDPos[i, :]))
                   for (x, y, i) in zip(outputb, UmD_0, 1:21)]

            #Female NC
            UfNC = [x + fxp.beta * (fxp.gamma * y + fxp.chi * z +
                                    fxp.lambdau_NC * par.sigma * sum(v0 .* SfNCPos[i, :]) +
                                    (1 - fxp.gamma - fxp.chi) * m)
                    for (x, y, z, m, i) in zip(outputb, UfD_0, UfPL_0, UfNC_0, 1:21)]

            #Female PL
            UfPL = [x + fxp.beta * (fxp.gamma * y + fxp.etaf * z +
                                    (1 - fxp.gamma - fxp.etaf) * m)
                    for (x, y, z, m) in zip(outputb, UfD_0, UfYC_0, UfPL_0)]
            #Female YC
            UfYC = [x + fxp.beta * (fxp.gamma * y + fxp.chi * z +
                                    fxp.lambdau_YC * par.sigma * sum(v0 .* SfYCPos[i, :]) +
                                    (1 - fxp.gamma - fxp.chi) * m)
                    for (x, y, z, m, i) in zip(outputb, UfD_0, UfPL_0, UfYC_0, 1:21)]
            #Female D
            UfD = [x + fxp.beta * ((1 - fxp.phi) * y +
                                   fxp.lambdau_YC * par.sigma * sum(v0 .* SfDPos[i, :]))
                   for (x, y, i) in zip(outputb, UfD_0, 1:21)]

            Pi = zeros(Float64, 21)
            for j = 1:kk
                dmNC = -1 .* SmNC .+ SmNC[:, j]
                dmYC = -1 .* SmYC .+ SmYC[:, j]
                dmD = -1 .* SmD .+ SmD[:, j]

                dfNC = -1 .* SfNC .+ SfNC[:, j]
                dfYC = -1 .* SfYC .+ SfYC[:, j]
                dfD = -1 .* SfD .+ SfD[:, j]

                #Vacancy Value Calculation
                Pi[j] = -fxp.c + fxp.beta * (
                    fxp.lambdau_NC * (1 - par.sigma) * sum(umNC_0 .* SmNCPos[:, j]) +
                    fxp.lambdau_YC * (1 - par.sigma) * (sum(umYC_0 .* SmYCPos[:, j]) + sum(umD_0 .* SmDPos[:, j])) +
                    fxp.lambdau_NC * (1 - par.sigma) * sum(ufNC_0 .* SfNCPos[:, j]) +
                    fxp.lambdau_YC * (1 - par.sigma) * (sum(ufYC_0 .* SfYCPos[:, j]) + sum(ufD_0 .* SfDPos[:, j])) +
                    fxp.lambdae * (1 - par.sigma) * (
                        sum(hmNC_0 .* dmNC .* logit(dmNC, fxp)) +
                        sum(hmYC_0 .* dmYC .* logit(dmYC, fxp)) +
                        sum(hmD_0 .* dmD .* logit(dmD, fxp)) +
                        sum(hfNC_0 .* dfNC .* logit(dfNC, fxp)) +
                        sum(hfYC_0 .* dfYC .* logit(dfYC, fxp)) +
                        sum(hfD_0 .* dfD .* logit(dfD, fxp))
                    )
                    + Pi0[j]
                )
            end

            # joint value P
            PmNC_hat = PmNC_0 .* (SmNC .> 0) + [x + y for x in UmNC, y in Pi] .* (SmNC .< 0)
            PmPL_hat = PmPL_0 .* (SmPL .> 0) + [x + y for x in UmPL, y in Pi] .* (SmPL .< 0)
            PmYC_hat = PmYC_0 .* (SmYC .> 0) + [x + y for x in UmYC, y in Pi] .* (SmYC .< 0)
            PmD_hat = PmD_0 .* (SmD .> 0) + [x + y for x in UmD, y in Pi] .* (SmD .< 0)

            PfNC_hat = PfNC_0 .* (SfNC .> 0) + [x + y for x in UfNC, y in Pi] .* (SfNC .< 0)
            PfPL_hat = PfPL_0 .* (SfPL .> 0) + [x + y for x in UfPL, y in Pi] .* (SfPL .< 0)
            PfYC_hat = PfYC_0 .* (SfYC .> 0) + [x + y for x in UfYC, y in Pi] .* (SfYC .< 0)
            PfD_hat = PfD_0 .* (SfD .> 0) + [x + y for x in UfD, y in Pi] .* (SfD .< 0)
            for i in 1:ii
                for j in 1:kk
                    PmNC[i, j] = getPmNC(i, j, v0, fxp, par, Pi, gem, SmNC, UmNC, PmD_hat, PmPL_hat, PmNC_hat, epsm, alpha, output, tau)
                    PmPL[i, j] = getPmPL(i, j, Pi, fxp, UmPL, PmD_hat, PmYC_hat, PmPL_hat, epsm, alpha, output)
                    PmYC[i, j] = getPmYC(i, j, v0, Pi, fxp, par, epsm, gem, tau, output, alpha, UmYC, SmYC, PmD_hat, PmPL_hat, PmYC_hat)
                    PmD[i, j] = getPmD(i, j, v0, Pi, fxp, par, output, gem, tau, epsm, alpha, SmD, PmD_hat, UmD)

                    PfNC[i, j] = getPfNC(i, j, v0, fxp, par, Pi, gef, SfNC, UfNC, PfD_hat, PfPL_hat, PfNC_hat, epsf, alpha, output, tau)
                    PfPL[i, j] = getPfPL(i, j, Pi, par, fxp, UfPL, PfD_hat, PfYC_hat, PfPL_hat, epsf, alpha, output)
                    PfYC[i, j] = getPfYC(i, j, v0, Pi, fxp, par, epsf, gef, tau, output, alpha, UfYC, SfYC, PfD_hat, PfPL_hat, PfYC_hat)
                    PfD[i, j] = getPfD(i, j, v0, Pi, fxp, par, output, gef, tau, epsf, alpha, SfD, PfD_hat, UfD)
                end
            end

            SmNC_0 = PmNC - [x + y for x in UmNC, y in Pi]
            SmPL_0 = PmPL - [x + y for x in UmPL, y in Pi]
            SmYC_0 = PmYC - [x + y for x in UmYC, y in Pi]
            SmD_0 = PmD - [x + y for x in UmD, y in Pi]

            SfNC_0 = PfNC - [x + y for x in UfNC, y in Pi]
            SfPL_0 = PfPL - [x + y for x in UfPL, y in Pi]
            SfYC_0 = PfYC - [x + y for x in UfYC, y in Pi]
            SfD_0 = PfD - [x + y for x in UfD, y in Pi]

            change1 = sum((SmNC - SmNC_0) .^ 2 + (SmYC - SmYC_0) .^ 2 + (SmPL - SmPL_0) .^ 2 + (SmD - SmD_0) .^ 2 + (SfNC - SfNC_0) .^ 2 + (SfYC - SfYC_0) .^ 2 + (SfPL - SfPL_0) .^ 2 + (SfD - SfD_0) .^ 2)

            if change1 < tol
                println("converged in " * string(it1) * " iterations, with change:" * string(change1))
                break
            end
            UmNC_0 = UmNC
            UmPL_0 = UmPL
            UmYC_0 = UmYC
            UmD_0 = UmD
            UfNC_0 = UfNC
            UfPL_0 = UfPL
            UfYC_0 = UfYC
            UfD_0 = UfD

            Pi0 = Pi

            PmNC_0 = PmNC
            PmYC_0 = PmYC
            PmPL_0 = PmPL
            PmD_0 = PmD
            PfNC_0 = PfNC
            PfYC_0 = PfYC
            PfPL_0 = PfPL
            PfD_0 = PfD

            SmNC = SmNC_0
            SmNCPos = copy(SmNC_0)
            SmNCPos[SmNCPos.<0] .= 0

            SmYC = SmYC_0
            SmYCPos = copy(SmYC)
            SmYCPos[SmYCPos.<0] .= 0

            SmPL = SmPL_0
            SmPLPos = copy(SmPL)
            SmPLPos[SmPLPos.<0] .= 0

            SmD = SmD_0
            SmDPos = copy(SmD)
            SmDPos[SmDPos.<0] .= 0

            SfNC = SfNC_0
            SfNCPos = copy(SfNC)
            SfNCPos[SfNCPos.<0] .= 0

            SfYC = SfYC_0
            SfYCPos = copy(SfYC)
            SfYCPos[SfYCPos.<0] .= 0

            SfPL = SfPL_0
            SmPLPos = copy(SmPL)
            SmPLPos[SmPLPos.<0] .= 0

            SfD = SfD_0
            SfDPos = copy(SfD)
            SfDPos[SfDPos.<0] .= 0
        end
        #For development purposes, will be deleted in the future
        if (cmp(purpose, "value_fn_dev") == 0)
            ValueDict = Dict("SmNC" => SmNC, "UfNC" => UfNC)
            get!(ValueDict, "SmPL", SmPL)
            get!(ValueDict, "SmYC", SmYC)
            get!(ValueDict, "SmD", SmD)

            get!(ValueDict, "SfNC", SfNC)
            get!(ValueDict, "SfPL", SfPL)
            get!(ValueDict, "SfYC", SfYC)
            get!(ValueDict, "SfD", SfD)

            get!(ValueDict, "UfPL", UfPL)
            get!(ValueDict, "UfYC", UfYC)
            get!(ValueDict, "UfD", UfD)

            get!(ValueDict, "UmNC", UmNC)
            get!(ValueDict, "UmPL", UmPL)
            get!(ValueDict, "UmYC", UmYC)
            get!(ValueDict, "UmD", UmD)

            get!(ValueDict, "Pi", Pi)
            return ValueDict
        end

        x0 = vcat(umNC_0, umPL_0, umYC_0, umD_0, ufNC_0, ufPL_0, ufYC_0, ufD_0, vec(hmNC_0), vec(hmYC_0), vec(hmPL_0), vec(hmD_0), vec(hfNC_0), vec(hfYC_0), vec(hfPL_0), vec(hfD_0), vec(v0))
        x_orig = x0

        #For development purposes, will be deleted in the future
        if (cmp(purpose, "solve_dist_dev") == 0)
            return fxp, SmNC, SmPL, SmYC, SmD, SfNC, SfPL, SfYC, SfD, e["stages"], x0
        end
        #For development purposes, will be deleted in the future
        if (cmp(purpose, "solve_dist_matrix") == 0)
            x = solveDist(x0, set, fxp, par, SmNC, SmPL, SmYC, SmD, SfNC, SfPL, SfYC, SfD, e["stages"])
            x[x.<1e-30] .= 0
            v0 = x[length(x)-set.Ny*set.Na+1:length(x)]
            #delete v0 from x
            x = deleteat!(x, length(x)-set.Ny*set.Na+1:length(x))
            x = reshape(x, set.Nx * set.Ne, :)

            umNC_0 = x[:, 1]
            umPL_0 = x[:, 2]
            umYC_0 = x[:, 3]
            umD_0 = x[:, 4]
            ufNC_0 = x[:, 5]
            ufPL_0 = x[:, 6]
            ufYC_0 = x[:, 7]
            ufD_0 = x[:, 8]

            hmNC_0 = x[:, 9:29]
            hmYC_0 = x[:, 30:50]
            hmPL_0 = x[:, 51:71]
            hmD_0 = x[:, 72:92]
            hfNC_0 = x[:, 93:113]
            hfYC_0 = x[:, 114:134]
            hfPL_0 = x[:, 135:155]
            hfD_0 = x[:, 156:176]

            solveDistMatrix = Dict("umNC_0" => umNC_0, "hmNC_0" => hmNC_0)
            get!(solveDistMatrix, "hmPL_0", hmPL_0)
            get!(solveDistMatrix, "hmYC_0", hmYC_0)
            get!(solveDistMatrix, "hmD_0", hmD_0)

            get!(solveDistMatrix, "hfNC_0", hfNC_0)
            get!(solveDistMatrix, "hfPL_0", hfPL_0)
            get!(solveDistMatrix, "hfYC_0", hfYC_0)
            get!(solveDistMatrix, "hfD_0", hfD_0)

            get!(solveDistMatrix, "umPL_0", umPL_0)
            get!(solveDistMatrix, "umYC_0", umYC_0)
            get!(solveDistMatrix, "umD_0", umD_0)

            get!(solveDistMatrix, "ufNC_0", ufNC_0)
            get!(solveDistMatrix, "ufPL_0", ufPL_0)
            get!(solveDistMatrix, "ufYC_0", ufYC_0)
            get!(solveDistMatrix, "ufD_0", ufD_0)

            return solveDistMatrix
        end

        for it2 = 1:maxit
            x = solveDist(x0, set, fxp, par, SmNC, SmPL, SmYC, SmD, SfNC, SfPL, SfYC, SfD, e["stages"])
            x[x.<1e-30] .= 0
            dev = (x - x0) ./ x0 #percentage deviation
            replace!(dev, Inf => NaN)
            replace!(dev, NaN => 0)
            change2 = sum(dev .^ 2)
            if (change2 < tol)
                println("converged in " * string(it2) * " iterations, with change:" * string(change2))
                break
            end
            x0 = copy(x)
        end
        #If this fixed point in distributions is different from the initial
        #distribution, then repeat the who process again
        dev = (x - x_orig) ./ x_orig        #percentage deviation
        replace!(dev, Inf => NaN)
        replace!(dev, NaN => 0)
        change_dist[itout] = sum(dev .^ 2)

        #Print Progress
        println("==================" * "Outer Iteration " * string(itout) * " with change " * string(change_dist[itout]) * "==================")
        #Update Solution for distributions
        v0 = x[length(x)-set.Ny*set.Na+1:length(x)]
        #delete v0 from x
        x = deleteat!(x, length(x)-set.Ny*set.Na+1:length(x))
        x = reshape(x, set.Nx * set.Ne, :)

        umNC_0 = x[:, 1]
        umPL_0 = x[:, 2]
        umYC_0 = x[:, 3]
        umD_0 = x[:, 4]
        ufNC_0 = x[:, 5]
        ufPL_0 = x[:, 6]
        ufYC_0 = x[:, 7]
        ufD_0 = x[:, 8]
        hmNC_0 = x[:, 9:29]
        hmYC_0 = x[:, 30:50]
        hmPL_0 = x[:, 51:71]
        hmD_0 = x[:, 72:92]
        hfNC_0 = x[:, 93:113]
        hfYC_0 = x[:, 114:134]
        hfPL_0 = x[:, 135:155]
        hfD_0 = x[:, 156:176]

        #output distribution
        e["hmNC"] = hmNC_0
        e["hmPL"] = hmPL_0
        e["hmYC"] = hmYC_0
        e["hmD"] = hmD_0

        e["hfNC"] = hfNC_0
        e["hfPL"] = hfPL_0
        e["hfYC"] = hfYC_0
        e["hfD"] = hfD_0

        e["umNC"] = umNC_0
        e["umPL"] = umPL_0
        e["umYC"] = umYC_0
        e["umD"] = umD_0

        e["ufNC"] = ufNC_0
        e["ufPL"] = ufPL_0
        e["ufYC"] = ufYC_0
        e["ufD"] = ufD_0

        e["v"] = v0

        #output surpluses
        e["SmNC"] = SmNC
        e["SmPL"] = SmPL
        e["SmYC"] = SmYC
        e["SmD"] = SmD

        e["SfNC"] = SfNC
        e["SfPL"] = SfPL
        e["SfYC"] = SfYC
        e["SfD"] = SfD

        e["UmNC"] = UmNC
        e["UmPL"] = UmPL
        e["UmYC"] = UmYC
        e["UmD"] = UmD

        e["UfNC"] = UfNC
        e["UfPL"] = UfPL
        e["UfYC"] = UfYC
        e["UfD"] = UfD

        e["Pi"] = Pi

        #Update Tax Rate
        tau = fxp.Transfer * sum(sum(hmPL_0 + hfPL_0, dims=1) .* output, dims=1) / (sum(sum(hmNC_0 + hfNC_0 + hmYC_0 + hfYC_0 + hmD_0 + hfD_0, dims=1) .* output, dims=1))
        e["tau"] = tau

        #for development purpose, will delete later
        if (cmp(purpose, "one_outer_loop") == 0)
            if (itout == 1)
                return fxp, e
            end
        end

        if (cmp(purpose, "two_outer_loop") == 0)
            if (itout == 2)
                return fxp, e
            end
        end

        #Convergence Check: bouncing between 2 iterations
        if (((itout > 3) && (abs(change_dist[itout-1] - change_dist[itout-3]) < tol / 100)) || itout > 30)
            get!(e, "Pineg", 0)
            get!(e, "par", par)
            get!(e, "fxp", fxp)
            findneg = (Pi .< 0)
            if (sum(findneg) > 0)
                get!(e, "Pineg", 1)
            end
            break
        end
        #converged
        if (change_dist[itout] < tol)
            println("After Iterations " * string(itout) * " converge with change " * string(change_dist[itout]))
            findneg = (Pi .< 0)
            if (sum(findneg) > 0)
                get!(e, "Pineg", 1)
            end
            break
        end
    end

    if (cmp(purpose, "solve_eqm_wo_wage") == 0)
        return fxp, e
    end

    #Solve for equilibrium wages
    # calculate wages even when didn't converge
    # if e["converge"] == 1 && e["Pineg"] == 0

    #quilibrium wage phi0 (update with surplus share)
    v0mat = repeat(vec(v0)', outer=[ii, 1])
    phi0m_NC = zeros(Float64, ii, kk)
    phi0m_PL = zeros(Float64, ii, kk)
    phi0m_YC = zeros(Float64, ii, kk)
    phi0m_D = zeros(Float64, ii, kk)
    phi0f_NC = zeros(Float64, ii, kk)
    phi0f_PL = zeros(Float64, ii, kk)
    phi0f_YC = zeros(Float64, ii, kk)
    phi0f_D = zeros(Float64, ii, kk)

    Wm_NC = [UmNC[x] + par.sigma * SmNC[x, y] for x in 1:21, y in 1:21] .* (SmNC .> 0)
    Wm_PL = [UmPL[x] + par.sigma * SmPL[x, y] for x in 1:21, y in 1:21] .* (SmPL .> 0)
    Wm_YC = [UmYC[x] + par.sigma * SmYC[x, y] for x in 1:21, y in 1:21] .* (SmYC .> 0)
    Wm_D = [UmD[x] + par.sigma * SmD[x, y] for x in 1:21, y in 1:21] .* (SmD .> 0)

    Wf_NC = [UfNC[x] + par.sigma * SfNC[x, y] for x in 1:21, y in 1:21] .* (SfNC .> 0)
    Wf_PL = [UfPL[x] + par.sigma * SfPL[x, y] for x in 1:21, y in 1:21] .* (SfPL .> 0)
    Wf_YC = [UfYC[x] + par.sigma * SfYC[x, y] for x in 1:21, y in 1:21] .* (SfYC .> 0)
    Wf_D = [UfD[x] + par.sigma * SfD[x, y] for x in 1:21, y in 1:21] .* (SfD .> 0)

    gem = zeros(set.Nx, set.Nx, set.Ny)
    gef = zeros(set.Nx, set.Nx, set.Ny)
    for y = 1:set.Ny
        for x0 = 1:set.Nx
            for x1 = 1:set.Nx
                if x1 == x0 + 1
                    gem[x0, x1, y] = par.d1m + par.d2m * fxp.p[y]
                    gef[x0, x1, y] = par.d1f + par.d2f * fxp.p[y]
                elseif x1 == x0 + 2
                    gem[x0, x1, y] = (par.d1m + par.d2m * fxp.p[y]) * fxp.d3
                    gef[x0, x1, y] = (par.d1f + par.d2f * fxp.p[y]) * fxp.d3
                end
            end
            gem[x0, x0, y] = 1 - sum(gem[x0, :, y])
            gef[x0, x0, y] = 1 - sum(gef[x0, :, y])

        end

        gem[set.Nx, set.Nx, y] = 1
        gef[set.Nx, set.Nx, y] = 1
    end
    A = ones(set.Nx, set.Nx)
    B = kron(Matrix(1.0I, set.Ne, set.Ne), A)
    gem = repeat(gem, outer=[set.Ne, set.Ne, set.Na])
    gem = gem .* repeat(B, outer=[1, 1, ii])
    gef = repeat(gef, outer=[set.Ne, set.Ne, set.Na])
    gef = gef .* repeat(B, outer=[1, 1, ii])
    gu = Matrix(1.0I, ii, ii)

    #phase 0
    for i = 1:ii
        for j = 1:kk
            gem_j = gem[i, :, j]     # for current HC i and job j, transition matrix from x to x'
            gef_j = gef[i, :, j]     # for current HC i and job j, transition matrix from x to x'
            gu_i = gu[i, :]

            #NC stage

            #Male
            vmI = repeat(SmNC[:, j] + UmNC, outer=[1, kk])   # max value incumbent willing to give to worker
            vmP = SmNC + repeat(UmNC, outer=[1, kk])       # poaching firm value
            dm = vmP - vmI
            Am = (vmP .> vmI) .* ((vmI + par.sigma * (vmP - vmI)) .* logit(dm, fxp) +
                                  vmI .* (1 .- logit(dm, fxp))) + (vmP .<= vmI) .* (vmP .> Wm_NC[i, j]) .* ((vmP + par.sigma * (vmI - vmP)) .* (1 .- logit(dm, fxp)) +
                                                                                vmP .* logit(dm, fxp)) + (vmP .<= Wm_NC[i, j]) .* (Wm_NC[i, j] .* (1 .- logit(dm, fxp)) +
                                                                                vmP .* logit(dm, fxp))

            #no wage update for HC growth
            Wmhat_NC = UmNC .* (SmNC[:, j] .<= 0) + Wm_NC[i, j] .* (SmNC[:, j] .> 0)
            Wmhat_PL = UmPL .* (SmPL[:, j] .<= 0) + Wm_NC[i, j] .* (SmPL[:, j] .> 0)
            Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm_NC[i, j] .* (SmD[:, j] .> 0)

            Em_jk = fxp.deltam_NC * UmNC + fxp.lambdae * sum(v0mat .* (Am - repeat(Wmhat_NC, outer=[1, kk])), dims=2) .* (SmNC[:, j] .> 0) +
                    fxp.gamma * Wmhat_D + fxp.chi * Wmhat_PL + (1 - fxp.deltam_NC - fxp.gamma - fxp.chi) * Wmhat_NC

            phi0m_NC[i, j] = Wm_NC[i, j] - epsm[i] * alpha[j] - first(fxp.beta * gem_j' * Em_jk)

            #female

            vfI = repeat(SfNC[:, j] + UfNC, outer=[1, kk])   # max value incumbent willing to give to worker
            vfP = SfNC + repeat(UfNC, outer=[1, kk])       # poaching firm value
            df = vfP - vfI

            Af = (vfP .> vfI) .* ((vfI + par.sigma * (vfP - vfI)) .* logit(df, fxp) +
                                  vfI .* (1 .- logit(df, fxp))) + (vfP .<= vfI) .* (vfP .> Wf_NC[i, j]) .* ((vfP + par.sigma * (vfI - vfP)) .* (1 .- logit(df, fxp)) +
                                                                                vfP .* logit(df, fxp)) + (vfP .<= Wf_NC[i, j]) .* (Wf_NC[i, j] .* (1 .- logit(df, fxp)) +
                                                                                vmP .* logit(df, fxp))

            #now wage update for HC growth
            Wfhat_NC = UfNC .* (SfNC[:, j] .<= 0) + Wf_NC[i, j] .* (SfNC[:, j] .> 0)
            Wfhat_PL = UfPL .* (SfPL[:, j] .<= 0) + Wf_NC[i, j] .* (SfPL[:, j] .> 0)
            Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf_NC[i, j] .* (SfD[:, j] .> 0)

            Ef_jk = fxp.deltaf_NC * UfNC + fxp.lambdae * sum(v0mat .* (Af - repeat(Wfhat_NC, outer=[1, kk])), dims=2) .* (SfNC[:, j] .> 0) +
                    fxp.gamma * Wfhat_D + fxp.chi * Wfhat_PL + (1 - fxp.deltaf_NC - fxp.gamma - fxp.chi) * Wfhat_NC

            phi0f_NC[i, j] = Wf_NC[i, j] - epsf[i] * alpha[j] - first(fxp.beta * gef_j' * Ef_jk)

            #PL stage wages never used

            #male
            Wmhat_PL = UmPL .* (SmPL[:, j] .<= 0) + Wm_PL[i, j] .* (SmPL[:, j] .> 0)
            Wmhat_YC = UmYC .* (SmYC[:, j] .<= 0) + Wm_PL[:, j] .* (SmYC[:, j] .> 0)
            Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm_PL[:, j] .* (SmD[:, j] .> 0)

            Em_jk = fxp.deltam_YC * UmPL + fxp.gamma * Wmhat_D + fxp.etam * Wmhat_YC + (1 - fxp.deltam_YC - fxp.gamma - fxp.etam) * Wmhat_PL
            phi0m_PL[i, j] = Wm_PL[i, j] - epsm[i] * alpha[j] - first(fxp.beta * gu_i' * Em_jk)

            #female
            Wfhat_PL = UfPL .* (SfPL[:, j] .<= 0) + Wf_PL[i, j] .* (SfPL[:, j] .> 0)
            Wfhat_YC = UfYC .* (SfYC[:, j] .<= 0) + Wf_PL[i, j] .* (SfYC[:, j] .> 0)
            Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf_PL[i, j] .* (SfD[:, j] .> 0)
            Ef_jk = fxp.deltaf_YC * UfPL + fxp.gamma * Wfhat_D + fxp.etaf * Wfhat_YC + (1 - fxp.deltaf_YC - fxp.gamma - fxp.etaf) * Wfhat_PL
            phi0f_PL[i, j] = Wf_PL[i, j] - (par.M + epsf[i]) * alpha[j] - first(fxp.beta * gu_i' * Ef_jk)

            #YC stage

            #male
            vmI = repeat(SmYC[:, j] + UmYC, outer=[1, kk])
            vmP = SmYC + repeat(UmYC, outer=[1, kk])
            dm = vmP - vmI
            Am = (vmP .> vmI) .* ((vmI + par.sigma * (vmP - vmI)) .* logit(dm, fxp) + vmI .* (1 .- logit(dm, fxp))) + (vmP .<= vmI) .* (vmP .> Wm_YC[i, j]) .* ((vmP + par.sigma * (vmI - vmP)) .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp)) + (vmP .<= Wm_YC[i, j]) .* (Wm_YC[i, j] .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp))

            #no wage update for HC growth
            Wmhat_YC = UmYC .* (SmYC[:, j] .<= 0) + Wm_YC[i, j] .* (SmYC[:, j] .> 0)
            Wmhat_PL = UmPL .* (SmPL[:, j] .<= 0) + Wm_YC[i, j] .* (SmPL[:, j] .> 0)
            Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm_YC[i, j] .* (SmD[:, j] .> 0)
            Em_jk = fxp.deltam_YC * UmYC + fxp.lambdae * sum(v0mat .* (Am - repeat(Wmhat_YC, outer=[1, kk])), dims=2) .* (SmYC[:, j] .> 0) + fxp.gamma * Wmhat_D + fxp.chi * Wmhat_PL + (1 - fxp.deltam_YC - fxp.gamma - fxp.chi) * Wmhat_YC
            phi0m_YC[i, j] = Wm_YC[i, j] - epsm[i] * alpha[j] - first(fxp.beta * gem_j' * Em_jk)

            #female
            vfI = repeat(SfYC[:, j] + UfYC, outer=[1, kk])
            vfP = SfYC + repeat(UfYC, outer=[1, kk])
            df = vfP - vfI
            Af = (vfP .> vfI) .* ((vfI + par.sigma * (vfP - vfI)) .* logit(df, fxp) + vfI .* (1 .- logit(df, fxp))) + (vfP .<= vfI) .* (vfP .> Wf_YC[i, j]) .* ((vfP + par.sigma * (vfI - vfP)) .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp)) + (vfP .<= Wf_YC[i, j]) .* (Wf_YC[i, j] .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp))

            Wfhat_YC = UfYC .* (SfYC[:, j] .<= 0) + Wf_YC[i, j] .* (SfYC[:, j] .> 0)
            Wfhat_PL = UfPL .* (SfPL[:, j] .<= 0) + Wf_YC[i, j] .* (SfPL[:, j] .> 0)
            Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf_YC[i, j] .* (SfD[:, j] .> 0)
            Ef_jk = fxp.deltaf_YC * UfYC + fxp.lambdae * sum(v0mat .* (Af - repeat(Wfhat_YC, outer=[1, kk])), dims=2) .* (SfYC[:, j] .> 0) + fxp.gamma * Wfhat_D + fxp.chi * Wfhat_PL + (1 - fxp.deltaf_YC - fxp.gamma - fxp.chi) * Wfhat_YC
            phi0f_YC[i, j] = Wf_YC[i, j] - (par.M + epsf[i]) * alpha[j] - first(fxp.beta * gef_j' * Ef_jk)

            #D stage
            #male
            vmI = repeat(SmD[:, j] + UmD, outer=[1, kk])
            vmP = SmD + repeat(UmD, outer=[1, kk])
            dm = vmP - vmI
            Am = (vmP .> vmI) .* ((vmI + par.sigma * (vmP - vmI)) .* logit(dm, fxp) + vmI .* (1 .- logit(dm, fxp))) + (vmP .<= vmI) .* (vmP .> Wm_D[i, j]) .* ((vmP + par.sigma * (vmI - vmP)) .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp)) .+ (vmP .<= Wm_D[i, j]) .* (Wm_D[i, j] .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp))

            Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm_D[i, j] .* (SmD[:, j] .> 0)
            Em_jk = fxp.delta * UmD + fxp.lambdae * sum(v0mat .* (Am - repeat(Wmhat_D, outer=[1, kk])), dims=2) .* (SmD[:, j] .> 0) + (1 - fxp.phi - fxp.delta) * Wmhat_D
            phi0m_D[i, j] = Wm_D[i, j] - epsm[i] * alpha[j] - first(fxp.beta * gem_j' * Em_jk)

            #Female
            vfI = repeat(SfD[:, j] + UfD, outer=[1, kk])
            vfP = SfD + repeat(UfD, outer=[1, kk])
            df = vfP - vfI
            Af = (vfP .> vfI) .* ((vfI + par.sigma * (vfP - vfI)) .* logit(df, fxp) + vfI .* (1 .- logit(df, fxp))) + (vfP .<= vfI) .* (vfP .> Wf_D[i, j]) .* ((vfP + par.sigma * (vfI - vfP)) .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp)) + (vfP .<= Wf_D[i, j]) .* (Wf_D[i, j] .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp))

            Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf_D[i, j] .* (SfD[:, j] .> 0)
            Ef_jk = fxp.delta * UfD + fxp.lambdae * sum(v0mat .* (Af - repeat(Wfhat_D, outer=[1, kk])), dims=2) .* (SfD[:, j] .> 0) + (1 - fxp.phi - fxp.delta) * Wfhat_D
            phi0f_D[i, j] = Wf_D[i, j] - epsf[i] * alpha[j] - first(fxp.beta * gef_j' * Ef_jk)
        end
    end
    phi0m_NC[SmNC.<=0] .= NaN
    phi0m_PL[SmPL.<=0] .= NaN
    phi0m_YC[SmYC.<=0] .= NaN
    phi0m_D[SmD.<=0] .= NaN
    phi0f_NC[SfNC.<=0] .= NaN
    phi0f_PL[SfPL.<=0] .= NaN
    phi0f_YC[SfYC.<=0] .= NaN
    phi0f_D[SfD.<=0] .= NaN

    #e qm wages phi1(y, y', x) -- wage at JTJ transitions from y to y'
    #make JTJ transition from y to y' when S(x, y')>S(x, y)
    Wm1_NC = zeros(Float64, kk, kk, ii)
    Wm1_YC = zeros(Float64, kk, kk, ii)
    Wm1_D = zeros(Float64, kk, kk, ii)
    Wf1_NC = zeros(Float64, kk, kk, ii)
    Wf1_YC = zeros(Float64, kk, kk, ii)
    Wf1_D = zeros(Float64, kk, kk, ii)
    phi1m_NC = zeros(Float64, kk, kk, ii)
    fill!(phi1m_NC, NaN)
    phi1m_YC = zeros(Float64, kk, kk, ii)
    fill!(phi1m_YC, NaN)
    phi1m_D = zeros(Float64, kk, kk, ii)
    fill!(phi1m_D, NaN)
    phi1f_NC = zeros(Float64, kk, kk, ii)
    fill!(phi1f_NC, NaN)
    phi1f_YC = zeros(Float64, kk, kk, ii)
    fill!(phi1f_YC, NaN)
    phi1f_D = zeros(Float64, kk, kk, ii)
    fill!(phi1f_D, NaN)
    #Wm_NC(k,j,i) where j is the poaching firm that succesfully poached from k (Incumbent)
    for i = 1:ii
        Sm_yP = repeat(SmNC[i, :]', outer=[kk, 1])
        Sm_yI = repeat(SmNC[i, :], outer=[1, kk])
        #JTJ transitions only happen when poaching>incumbent, and incumbent>0
        Wm1_NC[:, :, i] = (UmNC[i] * ones(kk, kk) + Sm_yI + par.sigma * (Sm_yP - Sm_yI)) .* (Sm_yP .> Sm_yI) .* (Sm_yI .> 0)

        Sm_yP = repeat(SmYC[i, :]', outer=[kk, 1])
        Sm_yI = repeat(SmYC[i, :], outer=[1, kk])
        #JTJ transitions only happen when poaching>incumbent, and incumbent>0
        Wm1_YC[:, :, i] = (UmYC[i] * ones(kk, kk) + Sm_yI + par.sigma * (Sm_yP - Sm_yI)) .* (Sm_yP .> Sm_yI) .* (Sm_yI .> 0)

        Sm_yP = repeat(SmD[i, :]', outer=[kk, 1])
        Sm_yI = repeat(SmD[i, :], outer=[1, kk])
        #JTJ transitions only happen when poaching>incumbent, and incumbent>0
        Wm1_D[:, :, i] = (UmD[i] * ones(kk, kk) + Sm_yI + par.sigma * (Sm_yP - Sm_yI)) .* (Sm_yP .> Sm_yI) .* (Sm_yI .> 0)

        Sf_yP = repeat(SfNC[i, :]', outer=[kk, 1])
        Sf_yI = repeat(SfNC[i, :], outer=[1, kk])
        #JTJ transitions only happen when poaching>incumbent, and incumbent>0
        Wf1_NC[:, :, i] = (UfNC[i] * ones(kk, kk) + Sf_yI + par.sigma * (Sf_yP - Sf_yI)) .* (Sf_yP .> Sf_yI) .* (Sf_yI .> 0)

        Sf_yP = repeat(SfYC[i, :]', outer=[kk, 1])
        Sf_yI = repeat(SfYC[i, :], outer=[1, kk])
        #JTJ transitions only happen when poaching>incumbent, and incumbent>0
        Wf1_YC[:, :, i] = (UfYC[i] * ones(kk, kk) + Sf_yI + par.sigma * (Sf_yP - Sf_yI)) .* (Sf_yP .> Sf_yI) .* (Sf_yI .> 0)

        Sf_yP = repeat(SfD[i, :]', outer=[kk, 1])
        Sf_yI = repeat(SfD[i, :], outer=[1, kk])
        #JTJ transitions only happen when poaching>incumbent, and incumbent>0
        Wf1_D[:, :, i] = (UfD[i] * ones(kk, kk) + Sf_yI + par.sigma * (Sf_yP - Sf_yI)) .* (Sf_yP .> Sf_yI) .* (Sf_yI .> 0)
    end
    #calculate expected future payoff for worker x at poaching firm y'
    #phi1(y,y',x) = phi1(k,j,i)
    for i = 1:ii
        for j = 1:kk
            for k = 1:kk
                gem_j = gem[i, :, j]
                gef_j = gef[i, :, j]

                #NC stage
                #male
                if Wm1_NC[k, j, i] != 0
                    vmI = repeat(SmNC[:, j] + UmNC, outer=[1, kk]) #max value incumbent willing to offer
                    vmP = SmNC + repeat(UmNC, outer=[1, kk])    #poaching firm value
                    dm = vmP - vmI
                    Am = (vmP .> vmI) .* ((vmI + par.sigma * (vmP - vmI)) .* logit(dm, fxp) + vmI .* (1 .- logit(dm, fxp))) + (vmP .<= vmI) .* (vmP .> Wm1_NC[k, j, i]) .* ((vmP + par.sigma * (vmI - vmP)) .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp)) + (vmP .<= Wm1_NC[k, j, i]) .* (Wm1_NC[k, j, i] .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp))

                    Wmhat_NC = UmNC .* (SmNC[:, j] .<= 0) + Wm1_NC[k, j, i] .* (SmNC[:, j] .> 0)
                    Wmhat_PL = UmPL .* (SmPL[:, j] .<= 0) + Wm1_NC[k, j, i] .* (SmPL[:, j] .> 0)
                    Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm1_NC[k, j, i] .* (SmD[:, j] .> 0)
                    Em_j = fxp.deltam_NC * UmNC + fxp.lambdae * sum(v0mat .* (Am - repeat(Wmhat_NC, outer=[1, kk])), dims=2) .* (SmNC[:, j] .> 0) + fxp.gamma * Wmhat_D + fxp.chi * Wmhat_PL + (1 - fxp.deltam_NC - fxp.gamma - fxp.chi) * Wmhat_NC
                    phi1m_NC[k, j, i] = Wm1_NC[k, j, i] - epsm[i] * alpha[j] - first(fxp.beta * gem_j' * Em_j)

                end
                if Wf1_NC[k, j, i] != 0
                    vfI = repeat(SfNC[:, j] + UfNC, outer=[1, kk]) #max value incumbent willing to offer
                    vfP = SfNC + repeat(UfNC, outer=[1, kk])    #poaching firm value
                    df = vfP - vfI

                    Af = (vfP .> vfI) .* ((vfI + par.sigma * (vfP - vfI)) .* logit(df, fxp) + vfI .* (1 .- logit(df, fxp))) + (vfP .<= vfI) .* (vfP .> Wf1_NC[k, j, i]) .* ((vfP + par.sigma * (vfI - vfP)) .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp)) + (vfP .<= Wf1_NC[k, j, i]) .* (Wf1_NC[k, j, i] .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp))

                    Wfhat_NC = UfNC .* (SfNC[:, j] .<= 0) + Wf1_NC[k, j, i] .* (SfNC[:, j] .> 0)
                    Wfhat_PL = UfPL .* (SfPL[:, j] .<= 0) + Wf1_NC[k, j, i] .* (SfPL[:, j] .> 0)
                    Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf1_NC[k, j, i] .* (SfD[:, j] .> 0)
                    Ef_j = fxp.deltaf_NC * UfNC + fxp.lambdae * sum(v0mat .* (Af - repeat(Wfhat_NC, outer=[1, kk])), dims=2) .* (SfNC[:, j] .> 0) + fxp.gamma * Wfhat_D + fxp.chi * Wfhat_PL + (1 - fxp.deltaf_NC - fxp.gamma - fxp.chi) * Wfhat_NC
                    phi1f_NC[k, j, i] = Wf1_NC[k, j, i] - epsf[i] * alpha[j] - first(fxp.beta * gef_j' * Ef_j)

                end

                #YC stage
                #Male
                if Wm1_YC[k, j, i] != 0
                    vmI = repeat(SmYC[:, j] + UmYC, outer=[1, kk])
                    vmP = SmYC + repeat(UmYC, outer=[1, kk])
                    dm = vmP - vmI

                    Am = (vmP .> vmI) .* ((vmI + par.sigma * (vmP - vmI)) .* logit(dm, fxp) + vmI .* (1 .- logit(dm, fxp))) + (vmP .<= vmI) .* (vmP .> Wm1_YC[k, j, i]) .* ((vmP + par.sigma * (vmI - vmP)) .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp)) + (vmP .<= Wm1_YC[k, j, i]) .* (Wm1_YC[k, j, i] .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp))

                    Wmhat_YC = UmYC .* (SmYC[:, j] .<= 0) + Wm1_YC[k, j, i] .* (SmYC[:, j] .> 0)
                    Wmhat_PL = UmPL .* (SmPL[:, j] .<= 0) + Wm1_YC[k, j, i] .* (SmPL[:, j] .> 0)
                    Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm1_YC[k, j, i] .* (SmD[:, j] .> 0)

                    Em_j = fxp.deltam_YC * UmYC + fxp.lambdae * sum(v0mat .* (Am - repeat(Wmhat_YC, outer=[1, kk])), dims=2) .* (SmYC[:, j] .> 0) + fxp.gamma * Wmhat_D + fxp.chi * Wmhat_PL + (1 - fxp.deltam_YC - fxp.gamma - fxp.chi) * Wmhat_YC
                    phi1m_YC[k, j, i] = Wm1_YC[k, j, i] - epsm[i] * alpha[j] - first(fxp.beta * gem_j' * Em_j)
                end

                #Female
                if Wf1_YC[k, j, i] != 0
                    vfI = repeat(SfYC[:, j] + UfYC, outer=[1, kk])
                    vfP = SfYC + repeat(UfYC, outer=[1, kk])
                    df = vfP - vfI

                    Af = (vfP .> vfI) .* ((vfI + par.sigma * (vfP - vfI)) .* logit(df, fxp) + vfI .* (1 .- logit(df, fxp))) + (vfP .<= vfI) .* (vfP .> Wf1_YC[k, j, i]) .* ((vfP + par.sigma * (vfI - vfP)) .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp)) + (vfP .<= Wf1_YC[k, j, i]) .* (Wf1_YC[k, j, i] .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp))

                    Wfhat_YC = UfYC .* (SfYC[:, j] .<= 0) + Wf1_YC[k, j, i] .* (SfYC[:, j] .> 0)
                    Wfhat_PL = UfPL .* (SfPL[:, j] .<= 0) + Wf1_YC[k, j, i] .* (SfPL[:, j] .> 0)
                    Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf1_YC[k, j, i] .* (SfD[:, j] .> 0)

                    Ef_j = fxp.deltaf_YC * UfYC + fxp.lambdae * sum(v0mat .* (Af - repeat(Wfhat_YC, outer=[1, kk])), dims=2) .* (SfYC[:, j] .> 0) + fxp.gamma * Wfhat_D + fxp.chi * Wfhat_PL + (1 - fxp.deltaf_YC - fxp.gamma - fxp.chi) * Wfhat_YC
                    phi1f_YC[k, j, i] = Wf1_YC[k, j, i] - (par.M + epsf[i]) * alpha[j] - first(fxp.beta * gef_j' * Ef_j)
                end

                #D stage
                #Male
                if Wm1_D[k, j, i] != 0
                    vmI = repeat(SmD[:, j] + UmD, outer=[1, kk])
                    vmP = SmD + repeat(UmD, outer=[1, kk])
                    dm = vmP - vmI

                    Am = (vmP .> vmI) .* ((vmI + par.sigma * (vmP - vmI)) .* logit(dm, fxp) + vmI .* (1 .- logit(dm, fxp))) + (vmP .<= vmI) .* (vmP .> Wm1_D[k, j, i]) .* ((vmP + par.sigma * (vmI - vmP)) .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp)) + (vmP .<= Wm1_D[k, j, i]) .* (Wm1_D[k, j, i] .* (1 .- logit(dm, fxp)) + vmP .* logit(dm, fxp))

                    Wmhat_D = UmD .* (SmD[:, j] .<= 0) + Wm1_D[k, j, i] .* (SmD[:, j] .> 0)

                    Em_j = fxp.delta * UmD + fxp.lambdae * sum(v0mat .* (Am - repeat(Wmhat_D, outer=[1, kk])), dims=2) .* (SmD[:, j] .> 0) + (1 - fxp.phi - fxp.delta) * Wmhat_D
                    phi1m_D[k, j, i] = Wm1_D[k, j, i] - epsm[i] * alpha[j] - first(fxp.beta * gem_j' * Em_j)
                end

                #Female
                if Wf1_D[k, j, i] != 0
                    vfI = repeat(SfD[:, j] + UfD, outer=[1, kk])
                    vfP = SfD + repeat(UfD, outer=[1, kk])
                    df = vfP - vfI

                    Af = (vfP .> vfI) .* ((vfI + par.sigma * (vfP - vfI)) .* logit(df, fxp) + vfI .* (1 .- logit(df, fxp))) + (vfP .<= vfI) .* (vfP .> Wf1_D[k, j, i]) .* ((vfP + par.sigma * (vfI - vfP)) .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp)) + (vfP .<= Wf1_D[k, j, i]) .* (Wf1_D[k, j, i] .* (1 .- logit(df, fxp)) + vfP .* logit(df, fxp))

                    Wfhat_D = UfD .* (SfD[:, j] .<= 0) + Wf1_D[k, j, i] .* (SfD[:, j] .> 0)

                    Ef_j = fxp.delta * UfD + fxp.lambdae * sum(v0mat .* (Af - repeat(Wfhat_D, outer=[1, kk])), dims=2) .* (SfD[:, j] .> 0) + (1 - fxp.phi - fxp.delta) * Wfhat_D
                    phi1f_D[k, j, i] = Wf1_D[k, j, i] - epsf[i] * alpha[j] - first(fxp.beta * gef_j' * Ef_j)
                end
            end
        end
    end
    #eqm wages phi2(y',y, x) -- when poaching firm triggers renegotiation
    #do not make JTJ transition: S(x, y')<S(x, y)

    Wm2_NC = zeros(Float64, kk, kk, ii)
    Wm2_YC = zeros(Float64, kk, kk, ii)
    Wm2_D = zeros(Float64, kk, kk, ii)
    Wf2_NC = zeros(Float64, kk, kk, ii)
    Wf2_YC = zeros(Float64, kk, kk, ii)
    Wf2_D = zeros(Float64, kk, kk, ii)
    phi2m_NC = zeros(Float64, kk, kk, ii)
    fill!(phi2m_NC, NaN)
    phi2m_YC = zeros(Float64, kk, kk, ii)
    fill!(phi2m_YC, NaN)
    phi2m_D = zeros(Float64, kk, kk, ii)
    fill!(phi2m_D, NaN)
    phi2f_NC = zeros(Float64, kk, kk, ii)
    fill!(phi2f_NC, NaN)
    phi2f_YC = zeros(Float64, kk, kk, ii)
    fill!(phi2f_YC, NaN)
    phi2f_D = zeros(Float64, kk, kk, ii)
    fill!(phi2f_D, NaN)

    for i = 1:ii
        #Male
        Sm_yP = repeat(SmNC[i, :], outer=[1, kk])
        Sm_yI = repeat(SmNC[i, :]', outer=[kk, 1])

        Wm2_NC[:, :, i] = (UmNC[i] * ones(kk, kk) + Sm_yP + par.sigma * (Sm_yI - Sm_yP)) .* (Sm_yP .<= Sm_yI) .* (Sm_yI .> 0) .* (Sm_yP .> 0)

        Sm_yP = repeat(SmYC[i, :], outer=[1, kk])
        Sm_yI = repeat(SmYC[i, :]', outer=[kk, 1])

        Wm2_YC[:, :, i] = (UmYC[i] * ones(kk, kk) + Sm_yP + par.sigma * (Sm_yI - Sm_yP)) .* (Sm_yP .<= Sm_yI) .* (Sm_yI .> 0) .* (Sm_yP .> 0)

        Sm_yP = repeat(SmD[i, :], outer=[1, kk])
        Sm_yI = repeat(SmD[i, :]', outer=[kk, 1])

        Wm2_D[:, :, i] = (UmD[i] * ones(kk, kk) + Sm_yP + par.sigma * (Sm_yI - Sm_yP)) .* (Sm_yP .<= Sm_yI) .* (Sm_yI .> 0) .* (Sm_yP .> 0)
        #Female
        Sf_yP = repeat(SfNC[i, :], outer=[1, kk])
        Sf_yI = repeat(SfNC[i, :]', outer=[kk, 1])

        Wf2_NC[:, :, i] = (UfNC[i] * ones(kk, kk) + Sf_yP + par.sigma * (Sf_yI - Sf_yP)) .* (Sf_yP .<= Sf_yI) .* (Sf_yI .> 0) .* (Sf_yP .> 0)

        Sf_yP = repeat(SfYC[i, :], outer=[1, kk])
        Sf_yI = repeat(SfYC[i, :]', outer=[kk, 1])

        Wf2_YC[:, :, i] = (UfYC[i] * ones(kk, kk) + Sf_yP + par.sigma * (Sf_yI - Sf_yP)) .* (Sf_yP .<= Sf_yI) .* (Sf_yI .> 0) .* (Sf_yP .> 0)

        Sf_yP = repeat(SfD[i, :], outer=[1, kk])
        Sf_yI = repeat(SfD[i, :]', outer=[kk, 1])

        Wf2_D[:, :, i] = (UfD[i] * ones(kk, kk) + Sf_yP + par.sigma * (Sf_yI - Sf_yP)) .* (Sf_yP .<= Sf_yI) .* (Sf_yI .> 0) .* (Sf_yP .> 0)
    end

    for i = 1:ii
        for j = 1:kk
            for k = 1:kk
                gem_j = gem[i,:,j];
                gef_j = gef[i,:,j];
                #NC
                #Male
                if Wm2_NC[k,j,i] != 0
                    vmI         = repeat(SmNC[:,j]+UmNC, outer=[1, kk]);
                    vmP         = SmNC + repeat(UmNC, outer=[1, kk]);
                    dm          = vmP-vmI;

                    Am          = (vmP .>vmI).*((vmI+ par.sigma*(vmP-vmI)).*logit(dm, fxp) + vmI.*(1 .-logit(dm, fxp))) + (vmP .<=vmI).*(vmP .>Wm2_NC[k,j,i]).*( (vmP+ par.sigma*(vmI-vmP)).*(1 .-logit(dm, fxp)) + vmP.*logit(dm, fxp)) + (vmP .<=Wm2_NC[k,j,i]).*( Wm2_NC[k,j,i].*(1 .-logit(dm, fxp)) + vmP.*logit(dm, fxp));

                    Wmhat_NC    = UmNC.*(SmNC[:,j].<=0)+Wm2_NC[k,j,i].*(SmNC[:,j].>0);
                    Wmhat_PL    = UmPL.*(SmPL[:,j].<=0)+Wm2_NC[k,j,i].*(SmPL[:,j].>0);
                    Wmhat_D     = UmD.*(SmD[:,j].<=0)+Wm2_NC[k,j,i].*(SmD[:,j].>0);
                    Em_j        = fxp.deltam_NC*UmNC + fxp.lambdae*sum(v0mat.*(Am-repeat(Wmhat_NC,outer=[1,kk])),dims=2).*(SmNC[:,j] .>0) + fxp.gamma*Wmhat_D + fxp.chi*Wmhat_PL + (1-fxp.deltam_NC-fxp.gamma-fxp.chi)*Wmhat_NC
                    phi2m_NC[k,j,i] = Wm2_NC[k,j,i] - epsm[i]*alpha[j]- first(fxp.beta*gem_j'*Em_j);
                end
                #Feale
                if Wf2_NC[k,j,i] != 0
                    vfI         = repeat(SfNC[:,j]+UfNC, outer=[1, kk]);
                    vfP         = SfNC + repeat(UfNC, outer=[1, kk]);
                    df          = vfP-vfI;

                    Af          = (vfP .>vfI).*((vfI+ par.sigma*(vfP-vfI)).*logit(df, fxp) + vfI.*(1 .-logit(df, fxp))) + (vfP .<=vfI).*(vfP .>Wf2_NC[k,j,i]).*( (vfP+ par.sigma*(vfI-vfP)).*(1 .-logit(df, fxp)) + vfP.*logit(df, fxp)) + (vfP .<=Wf2_NC[k,j,i]).*( Wf2_NC[k,j,i].*(1 .-logit(df, fxp)) + vfP.*logit(df, fxp));

                    Wfhat_NC    = UfNC.*(SfNC[:,j].<=0)+Wf2_NC[k,j,i].*(SfNC[:,j].>0);
                    Wfhat_PL    = UfPL.*(SfPL[:,j].<=0)+Wf2_NC[k,j,i].*(SfPL[:,j].>0);
                    Wfhat_D     = UfD.*(SfD[:,j].<=0)+Wf2_NC[k,j,i].*(SfD[:,j].>0);
                    Ef_j        = fxp.deltaf_NC*UfNC + fxp.lambdae*sum(v0mat.*(Af-repeat(Wfhat_NC,outer=[1,kk])),dims=2).*(SfNC[:,j] .>0) + fxp.gamma*Wfhat_D + fxp.chi*Wfhat_PL + (1-fxp.deltaf_NC-fxp.gamma-fxp.chi)*Wfhat_NC
                    phi2f_NC[k,j,i] = Wf2_NC[k,j,i] - epsf[i]*alpha[j]- first(fxp.beta*gef_j'*Ef_j)
                end

                #YC
                #Male
                if Wm2_YC[k,j,i] != 0
                    vmI         = repeat(SmYC[:,j]+UmYC, outer=[1, kk]);
                    vmP         = SmYC + repeat(UmYC, outer=[1, kk]);
                    dm          = vmP-vmI;

                    Am          = (vmP .>vmI).*((vmI+ par.sigma*(vmP-vmI)).*logit(dm, fxp) + vmI.*(1 .-logit(dm, fxp))) + (vmP .<=vmI).*(vmP .>Wm2_YC[k,j,i]).*( (vmP+ par.sigma*(vmI-vmP)).*(1 .-logit(dm, fxp)) + vmP.*logit(dm, fxp)) + (vmP .<=Wm2_YC[k,j,i]).*( Wm2_YC[k,j,i].*(1 .-logit(dm, fxp)) + vmP.*logit(dm, fxp));

                    Wmhat_YC    = UmYC.*(SmYC[:,j].<=0)+Wm2_YC[k,j,i].*(SmYC[:,j].>0);
                    Wmhat_PL    = UmPL.*(SmPL[:,j].<=0)+Wm2_YC[k,j,i].*(SmPL[:,j].>0);
                    Wmhat_D     = UmD.*(SmD[:,j].<=0)+Wm2_YC[k,j,i].*(SmD[:,j].>0);
                    Em_j        = fxp.deltam_YC*UmYC + fxp.lambdae*sum(v0mat.*(Am-repeat(Wmhat_YC,outer=[1,kk])),dims=2).*(SmNC[:,j] .>0) + fxp.gamma*Wmhat_D + fxp.chi*Wmhat_PL + (1-fxp.deltam_YC-fxp.gamma-fxp.chi)*Wmhat_YC
                    phi2m_YC[k,j,i] = Wm2_YC[k,j,i] - epsm[i]*alpha[j]- first(fxp.beta*gem_j'*Em_j)
                end

                #Female
                if Wf2_YC[k,j,i] != 0
                    vfI         = repeat(SfYC[:,j]+UfYC, outer=[1, kk]);
                    vfP         = SfYC + repeat(UfYC, outer=[1, kk]);
                    df          = vfP-vfI;

                    Af          = (vfP .>vfI).*((vfI+ par.sigma*(vfP-vfI)).*logit(df, fxp) + vfI.*(1 .-logit(df, fxp))) + (vfP .<=vfI).*(vfP .>Wf2_YC[k,j,i]).*( (vfP+ par.sigma*(vfI-vfP)).*(1 .-logit(df, fxp)) + vfP.*logit(df, fxp)) + (vfP .<=Wf2_YC[k,j,i]).*( Wf2_YC[k,j,i].*(1 .-logit(df, fxp)) + vfP.*logit(df, fxp));

                    Wfhat_YC    = UfYC.*(SfYC[:,j].<=0)+Wf2_YC[k,j,i].*(SfYC[:,j].>0);
                    Wfhat_PL    = UfPL.*(SfPL[:,j].<=0)+Wf2_YC[k,j,i].*(SfPL[:,j].>0);
                    Wfhat_D     = UfD.*(SfD[:,j].<=0)+Wf2_YC[k,j,i].*(SfD[:,j].>0);
                    Ef_j        = fxp.deltaf_YC*UfYC + fxp.lambdae*sum(v0mat.*(Af-repeat(Wfhat_YC,outer=[1,kk])),dims=2).*(SfNC[:,j] .>0) + fxp.gamma*Wfhat_D + fxp.chi*Wfhat_PL + (1-fxp.deltaf_YC-fxp.gamma-fxp.chi)*Wfhat_YC
                    phi2f_YC[k,j,i] = Wf2_YC[k,j,i] - (par.M+epsf[i])*alpha[j]- first(fxp.beta*gef_j'*Ef_j)
                end

                #D stage
                #Male
                if Wm2_D[k,j,i] != 0
                    vmI         = repeat(SmD[:,j]+UmD, outer=[1, kk]);
                    vmP         = SmD + repeat(UmD, outer=[1, kk]);
                    dm          = vmP-vmI;

                    Am          = (vmP .>vmI).*((vmI+ par.sigma*(vmP-vmI)).*logit(dm, fxp) + vmI.*(1 .-logit(dm, fxp))) + (vmP .<=vmI).*(vmP .>Wm2_D[k,j,i]).*( (vmP+ par.sigma*(vmI-vmP)).*(1 .-logit(dm, fxp)) + vmP.*logit(dm, fxp)) + (vmP .<=Wm2_D[k,j,i]).*( Wm2_D[k,j,i].*(1 .-logit(dm, fxp)) + vmP.*logit(dm, fxp));

                    Wmhat_D     = UmD.*(SmD[:,j].<=0)+Wm2_D[k,j,i].*(SmD[:,j].>0);
                    Em_j        = fxp.delta*UmD + fxp.lambdae*sum(v0mat.*(Am-repeat(Wmhat_D,outer=[1,kk])),dims=2).*(SmD[:,j] .>0) + (1-fxp.phi-fxp.delta)*Wmhat_D;
                    phi2m_D[k,j,i] = Wm2_D[k,j,i] - epsm[i]*alpha[j]- first(fxp.beta*gem_j'*Em_j)
                end

                #Female
                if Wf2_D[k,j,i] != 0
                    vfI         = repeat(SfD[:,j]+UfD, outer=[1, kk]);
                    vfP         = SfD + repeat(UfD, outer=[1, kk]);
                    df          = vfP-vfI;

                    Af          = (vfP .>vfI).*((vfI+ par.sigma*(vfP-vfI)).*logit(df, fxp) + vfI.*(1 .-logit(df, fxp))) + (vfP .<=vfI).*(vfP .>Wf2_D[k,j,i]).*( (vfP+ par.sigma*(vfI-vfP)).*(1 .-logit(df, fxp)) + vfP.*logit(df, fxp)) + (vfP .<=Wf2_D[k,j,i]).*( Wf2_D[k,j,i].*(1 .-logit(df, fxp)) + vfP.*logit(df, fxp));

                    Wfhat_D     = UfD.*(SfD[:,j].<=0)+Wf2_D[k,j,i].*(SfD[:,j].>0);
                    Ef_j        = fxp.delta*UfD + fxp.lambdae*sum(v0mat.*(Af-repeat(Wfhat_D,outer=[1,kk])),dims=2).*(SfD[:,j] .>0) + (1-fxp.phi-fxp.delta)*Wfhat_D;

                    phi2f_D[k,j,i] = Wf2_D[k,j,i] - epsf[i]*alpha[j]- first(fxp.beta*gef_j'*Ef_j)
                end
            end
        end
    end
    get!(e,"phi0m_NC",phi0m_NC);
    get!(e,"phi0m_PL",phi0m_PL);
    get!(e,"phi0m_YC",phi0m_YC);
    get!(e,"phi0m_D",phi0m_D);

    get!(e,"phi0f_NC",phi0f_NC);
    get!(e,"phi0f_PL",phi0f_PL);
    get!(e,"phi0f_YC",phi0f_YC);
    get!(e,"phi0f_D",phi0f_D);

    get!(e,"phi1m_NC",phi1m_NC);
    get!(e,"phi1m_YC",phi1m_YC);
    get!(e,"phi1m_D",phi1m_D);

    get!(e,"phi1f_NC",phi1f_NC);
    get!(e,"phi1f_YC",phi1f_YC);
    get!(e,"phi1f_D",phi1f_D);

    get!(e,"phi2m_NC",phi2m_NC);
    get!(e,"phi2m_YC",phi2m_YC);
    get!(e,"phi2m_D",phi2m_D);

    get!(e,"phi2f_NC",phi2f_NC);
    get!(e,"phi2f_YC",phi2f_YC);
    get!(e,"phi2f_D",phi2f_D);

    return e
end
export solve
end


text:
FRONT_VIEW = false

clear!(kps4)
kps4.stiffness_factor = 0.04

y0, yd0 = KiteModels.init(kps4)

find_steady_state!(kps4, prn=false)

if typeof(kps4) <: KPS4
    println("kite distance: $(norm(kps4.pos[end-2]))")
else
    println("kite distance: $(norm(kps4.pos[end]))")
end
println(KiteModels.spring_forces(kps4))
println("alpha_depower [deg]: $(rad2deg(kps4.alpha_depower))")
println("lift, drag    [N]  : $(KiteModels.lift_drag(kps4))")

include("plot2d.jl")
plot2d(kps4.pos; zoom=false, front=FRONT_VIEW)


text:
@testset "Topology" begin
    topology = Topology()
    @test size(topology) == 0

    # Creating some H2O2
    add_atom!(topology, Atom("H"))
    add_atom!(topology, Atom("O"))
    add_atom!(topology, Atom("O"))
    add_atom!(topology, Atom("H"))
    @test size(topology) == 4

    @test name(Atom(topology, 0)) == "H"
    @test name(Atom(topology, 1)) == "O"
    @test name(Atom(topology, 2)) == "O"
    @test name(Atom(topology, 3)) == "H"

    @test bonds_count(topology) == 0
    @test angles_count(topology) == 0
    @test dihedrals_count(topology) == 0

    add_bond!(topology, 0, 1)
    add_bond!(topology, 1, 2)
    add_bond!(topology, 2, 3, Chemfiles.AromaticBond)

    @test bonds_count(topology) == 3
    @test angles_count(topology) == 2
    @test dihedrals_count(topology) == 1

    @test bond_order(topology, 1, 2) == Chemfiles.UnknownBond
    @test bond_order(topology, 2, 3) == Chemfiles.AromaticBond

    @test bond_orders(topology) == [Chemfiles.UnknownBond, Chemfiles.UnknownBond, Chemfiles.AromaticBond]

    top_bonds = reshape(UInt64[0, 1,   1, 2,   2, 3], (2, 3))

    @test bonds(topology) == top_bonds
    @test angles(topology) == reshape(UInt64[0, 1, 2,   1, 2, 3,], (3, 2))
    @test dihedrals(topology) == reshape(UInt64[0, 1, 2, 3], (4,1))

    remove_bond!(topology, 2, 3)
    @test bonds_count(topology) == 2
    @test angles_count(topology) == 1
    @test dihedrals_count(topology) == 0

    clear_bonds!(topology)
    @test bonds_count(topology) == 0

    remove_atom!(topology, 3)
    @test size(topology) == 3

    resize!(topology, 42)
    @test size(topology) == 42

    copy = deepcopy(topology)
    @test size(copy) == 42

    resize!(copy, 25)
    @test size(copy) == 25
    @test size(topology) == 42

    @testset "Residues" begin
        topology = Topology()
        for i = 1:10
            add_atom!(topology, Atom("X"))
        end

        for atoms in [[2,3,6], [0,1,9], [4,5,8]]
            residue = Residue("X")
            for i in atoms
                add_atom!(residue, i)
            end
            add_residue!(topology, residue)
        end

        @test count_residues(topology) == 3

        first = residue_for_atom(topology, 2)
        second = residue_for_atom(topology, 0)

        @test residue_for_atom(topology, 7) === nothing

        @test first !== nothing
        @test second !== nothing
        remove_chemfiles_warning() do
            @test_throws ChemfilesError Residue(topology, 4)
        end

        @test are_linked(topology, first, second) == false

        add_bond!(topology, 6, 9)
        @test are_linked(topology, first, second) == true
    end

    @testset "Impropers" begin
        topology = Topology()
        add_atom!(topology, Atom("N"))
        add_atom!(topology, Atom("H"))
        add_atom!(topology, Atom("H"))
        add_atom!(topology, Atom("H"))

        add_bond!(topology, 0, 1)
        add_bond!(topology, 0, 2)
        add_bond!(topology, 0, 3)

        @test impropers_count(topology) == 1
        @test impropers(topology) == reshape(UInt64[1, 0, 2, 3], (4,1))
    end

    @testset "Topology indexing" begin
        topology = Topology()
        add_atom!(topology, Atom("N"))
        add_atom!(topology, Atom("H"))

        @test name(topology[0]) == "N"
        @test name(topology[1]) == "H"

        # @view is required to reference an atom without create a copy
        # Modifying an atom requires view, as normal indexing creates a copy
        set_type!(topology[0], "O")
        @test type(topology[0]) == "N"
        set_type!(@view(topology[0]), "O")
        @test type(topology[0]) == "O"
    end
end


text:
using Test

using CompScienceMeshes

e0 = point(0,0,0)
e1 = point(1,0,0)
e2 = point(0,1,0)
e3 = point(0,0,1)

p = simplex(e0, e1)
q = simplex((e0+e1)/2, (-e0+3*e1)/2)

pq = intersection(p,q)
@test length(pq) == 1
pq = pq[1]
@test volume(pq) == 1/2
@test pq.vertices[1] == q.vertices[1]
@test pq.vertices[2] == p.vertices[2]
@test 2*pq.tangents[1] == p.tangents[1] == q.tangents[1]

c = (e0+e1+e2)/3
p = simplex(e0, e1, e2)
q = simplex(e0, e1, c)

qq = intersection(q,q)
pq = intersection(p,q)
qp = intersection(q,p)

# For these case where points of the subject are on the clipper or
# the other way around, the precise computationof the intersection
# is subject to numerical rounding errors. For efficiency reasons
# this something that needs fixing. It should not affect the computation
# of the interaction elements, however.
@test sum([volume(r) for r in qq]) ≈ volume(q)
@test sum([volume(r) for r in pq]) ≈ volume(q)
@test sum([volume(r) for r in qp]) ≈ volume(q)

# This case caused NaNs in the RWG BC Gram matrix. Turned out to be caused
# by computing the intersection between two parallel lines. This needs a
# more robust solution in the long run.
p = point(-0.819550109062379,-0.394674531365384,-0.415415013001886)
q = point(-0.66680743687093,-0.618706884644307,-0.415415013001886)
r = point(-0.832690784854388,-0.535137873473586,-0.142314838273285)

a = point(-0.66680743687093,-0.618706884644307,-0.415415013001886)
b = point(-0.749749110862659,-0.5769223790589465,-0.27886492563758547)
c = point(-0.7730161102625658,-0.5161730964944257,-0.32438162142568566)

isct1 = CompScienceMeshes.sutherlandhodgman([p,q,r],[a,b,c])

@test isct1 == [
    point(-0.7730161102625658,-0.5161730964944257,-0.32438162142568566),
    point(-0.66680743687093,-0.618706884644307,-0.415415013001886),
    point(-0.66680743687093,-0.618706884644307,-0.415415013001886),
    point(-0.749749110862659,-0.5769223790589465,-0.27886492563758547)]

splx1 = simplex(
    point(0,0,0),
    point(1,0,0),
    point(0,1,0))
splx2 = simplex(
    point(0,-1,0),
    point(1,1,0),
    point(0.5,1,0))

# Verify that all parts of the intersection are oriented consistently
# and equal to the orientation of the inputs
isct2 = intersection(splx1, splx2)
@test length(isct2) == 2

for splx in isct2
    @test normal(splx) ≈ point(0,0,1)
end

# Repeat this test but now with the orientation of the second
# operand flipped
splx1 = simplex(
    point(0,0,0),
    point(1,0,0),
    point(0,1,0))
splx2 = simplex(
    point(1,1,0),
    point(0,-1,0),
    point(0.5,1,0))

isct2 = intersection(splx1, splx2)
@test length(isct2) == 2

for splx in isct2
    @test normal(splx) ≈ point(0,0,1)
end

# Repeat this test but now with the orientation of the first operand flipped
splx1 = simplex(
    point(1,0,0),
    point(0,0,0),
    point(0,1,0))
splx2 = simplex(
    point(0,-1,0),
    point(1,1,0),
    point(0.5,1,0))

isct2 = intersection(splx1, splx2)
@test length(isct2) == 2

for splx in isct2
    @test normal(splx) ≈ point(0,0,-1)
end