text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

const MortarElements3D = Union{Tri3,Tri6,Quad4}

function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector)
    return p - dot(p-x0, n0)*n0
end

function inv3(P::Matrix)
    n, m = size(P)
    @assert n == m == 3
    a, b, c, d, e, f, g, h, i = P
    A = e*i - f*h
    B = -d*i + f*g
    C = d*h - e*g
    D = -b*i + c*h
    E = a*i - c*g
    F = -a*h + b*g
    G = b*f - c*e
    H = -a*f + c*d
    I = a*e - b*d
    return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I]
end

function vertex_inside_polygon(q, P; atol=1.0e-3)
    N = length(P)
    angle = 0.0
    for i=1:N
        A = P[i] - q
        B = P[mod(i,N)+1] - q
        c = norm(A)*norm(B)
        isapprox(c, 0.0; atol=atol) && return true
        cosa = dot(A,B)/c
        isapprox(cosa, 1.0; atol=atol) && return false
        isapprox(cosa, -1.0; atol=atol) && return true
        #try
        angle += acos(cosa)
        #catch
        #    @info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
        #    @info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
        #    @info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
        #    @info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
        #    rethrow()
        #end
    end
    return isapprox(angle, 2*pi; atol=atol)
end

function calculate_centroid(P)
    N = length(P)
    P0 = P[1]
    areas = [norm(1/2*cross(P[i]-P0, P[mod(i,N)+1]-P0)) for i=2:N]
    centroids = [1/3*(P0+P[i]+P[mod(i,N)+1]) for i=2:N]
    C = 1/sum(areas)*sum(areas.*centroids)
    return C
end

function get_cells(P, C; allow_quads=false)
    N = length(P)
    cells = Vector[]
    # shared edge etc.
    N < 3 && return cells
    # trivial cases, polygon already triangle / quadrangle
    if N == 3
        return Vector[P]
    end
    if N == 4 && allow_quads
        return Vector[P]
    end
    cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N]
    return cells
end

""" Test does vector P contain approximately q. This function uses isapprox()
internally to make boolean test.

Examples
--------
julia> P = Vector[[1.0, 1.0], [2.0, 2.0]]
2-element Array{Array{T,1},1}:
 [1.0,1.0]
 [2.0,2.0]

julia> q = [1.0, 1.0] + eps(Float64)
2-element Array{Float64,1}:
 1.0
 1.0

julia> in(q, P)
false

julia> approx_in(q, P)
true

"""
function approx_in(q::T, P::Vector{T}; rtol=1.0e-4, atol=0.0) where T
    for p in P
        if isapprox(q, p; rtol=rtol, atol=atol)
            return true
        end
    end
    return false
end

function get_polygon_clip(xs::Vector{T}, xm::Vector{T}, n::T) where T
    # objective: search does line xm1 - xm2 clip xs
    nm = length(xm)
    ns = length(xs)
    P = T[]

    # 1. test is master point inside slave, if yes, add to clip
    for i=1:nm
        if vertex_inside_polygon(xm[i], xs)
            push!(P, xm[i])
        end
    end

    # 2. test is slave point inside master, if yes, add to clip
    for i=1:ns
        if vertex_inside_polygon(xs[i], xm)
            approx_in(xs[i], P) && continue
            push!(P, xs[i])
        end
    end

    for i=1:nm
        # 2. find possible intersection
        xm1 = xm[i]
        xm2 = xm[mod(i,nm)+1]
        # @info("intersecting line $xm1 -> $xm2")
        for j=1:ns
            xs1 = xs[j]
            xs2 = xs[mod(j,ns)+1]
            # @info("clipping polygon edge $xs1 -> $xs2")
            tnom = dot(cross(xm1-xs1, xm2-xm1), n)
            tdenom = dot(cross(xs2-xs1, xm2-xm1), n)
            isapprox(tdenom, 0) && continue
            t = tnom/tdenom
            (0 <= t <= 1) || continue
            q = xs1 + t*(xs2 - xs1)
            # @info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
            if vertex_inside_polygon(q, xm)
                approx_in(q, P) && continue
                push!(P, q)
            end
        end
    end

    return P
end

""" Project some vertex p to surface of element E using Newton's iterations. """
function project_vertex_to_surface(p, x0, n0,
                                   element::Element{E}, x, time;
                                   max_iterations=10, iter_tol=1.0e-6) where E
    basis(xi) = get_basis(element, xi, time)
    function dbasis(xi)
        return get_dbasis(element, xi, time)
    end
    nnodes = length(element)
    mul(a,b) = sum((a[:,i]*b[i]')' for i=1:length(b))

    function f(theta)
        b = [basis(theta[1:2])*collect(x)...;]
        b = b - theta[3]*n0 - p
        return b
    end

    L(theta) = inv3([mul(dbasis(theta[1:2]), x) -n0])
    theta = zeros(3)
    dtheta = zeros(3)
    for i=1:max_iterations
        invA = L(theta)
        b = f(theta)
        dtheta = invA * b
        theta -= dtheta
        if norm(dtheta) < iter_tol
            return theta[1:2], theta[3]
        end
    end
    #=
    @info("failed to project vertex from auxiliary plane back to surface")
    @info("element type: $E")
    @info("element connectivity: $(get_connectivity(element))")
    @info("auxiliary plane: x0 = $x0, n0 = $n0")
    @info("element geometry: $(x.data)")
    @info("vertex to project: $p")
    @info("parameter vector before giving up: $theta")
    @info("increment in parameter vector before giving up: $dtheta")
    @info("norm(dtheta) before giving up: $(norm(dtheta))")
    @info("f([0.0, 0.0, 0.0]) = $(f([0.0, 0.0, 0.0]))")
    @info("L([0.0, 0.0, 0.0]) = $(L([0.0, 0.0, 0.0]))")

    @info("iterations:")
    theta = zeros(3)
    dtheta = zeros(3)
    for i=1:max_iterations
        @info("iter $i, theta = $theta")
        @info("f = $(f(theta))")
        @info("L = $(L(theta))")
        dtheta = L(theta) * f(theta)
        @info("dtheta = $(dtheta)")
        theta -= dtheta
    end
    =#
    throw(error("project_point_to_surface: did not converge in $max_iterations iterations!"))
end

function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false)
    normals = Dict{Int64, Vector{Float64}}()
    for element in elements
        conn = get_connectivity(element)
        J = transpose(element([0.0, 0.0], time, Val{:Jacobian}))
        normal = cross(J[:,1], J[:,2])
        for nid in conn
            if haskey(normals, nid)
                normals[nid] += normal
            else
                normals[nid] = normal
            end
        end
    end
    # normalize to unit normal
    S = collect(keys(normals))
    for j in S
        normals[j] /= norm(normals[j])
    end
    if rotate_normals
        for j in S
            normals[j] = -normals[j]
        end
    end
    return normals
end

""" Given polygon P and normal direction n, check that polygon vertices are
ordered in counter clock wise direction with respect to surface normal and
sort if necessary. It is assumed that polygon is convex.

Examples
--------
Unit triangle, normal in z-direction:

julia> P = Vector[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]]
3-element Array{Array{T,1},1}:
 [0.0,0.0,0.0]
 [0.0,1.0,0.0]
 [1.0,0.0,0.0]

julia> n = [0.0, 0.0, 1.0]
3-element Array{Float64,1}:
 0.0
 0.0
 1.0

julia> check_orientation!(P, n)
3-element Array{Array{T,1},1}:
 [1.0,0.0,0.0]
 [0.0,0.0,0.0]
 [0.0,1.0,0.0]

"""
function check_orientation!(P, n)
    C = mean(P)
    np = length(P)
    s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np]
    all(s .< 0) && return
    # project points to new orthogonal basis Q and sort there
    t1 = (P[1]-C)/norm(P[1]-C)
    t2 = cross(n, t1)
    Q = [n t1 t2]
    sort!(P, lt=(A, B) -> begin
        A_proj = Q'*(A-C)
        B_proj = Q'*(B-C)
        a = atan(A_proj[3], A_proj[2])
        b = atan(B_proj[3], B_proj[2])
        return a > b
    end)
end

function convert_to_linear_element(element::Element{E}) where E
    return element
end

function convert_to_linear_element(element::Element{Tri6})
    new_element = Element(Tri3, element.connectivity[1:3])
    new_element.id = element.id
    new_element.fields = element.fields
    return new_element
end

function split_quadratic_element(element::Element{E}, time::Float64) where E
    return [element]
end

function split_quadratic_element(element::Element{Tri6}, time::Float64)
    element_maps = Vector{Int}[[1,4,6], [4,5,6], [4,2,5], [6,5,3]]
    new_elements = Element[]
    connectivity = get_connectivity(element)
    for elmap in element_maps
        new_element = Element(Tri3, connectivity[elmap])
        X = element("geometry", time)
        update!(new_element, "geometry", time => X[elmap])
        if haskey(element, "displacement")
            u = element("displacement", time)
            update!(new_element, "displacement", time => u[elmap])
        end
        if haskey(element, "normal")
            n = element("normal", time)
            update!(new_element, "normal", time => n[elmap])
        end
        push!(new_elements, new_element)
    end
    return new_elements
end

function split_quadratic_elements(elements::DVTI, time::Float64)
    return DVTI(split_quadratic_elements(elements.data, time))
end

""" Split quadratic surface elements to linear elements.  """
function split_quadratic_elements(elements::Vector, time::Float64)
    new_elements = Element[]
    for element in elements
        for splitted_element in split_quadratic_element(element, time)
            push!(new_elements, splitted_element)
        end
    end
    n1 = length(elements)
    n2 = length(new_elements)
    if n1 != n2
        @info("Splitted $n1 elements to $n2 (linear) sub-elements")
    end
    return new_elements
end

function get_mean_xi(element::Element)
    xi = zeros(2)
    coords = get_reference_coordinates(element)
    for (xi1,xi2) in coords
        xi[1] += xi1
        xi[2] += xi2
    end
    xi /= length(coords)
    return xi
end

""" Assemble linear surface element to problem.

Dual basis is constructed such that partially integrated slave segments are taken into account in a proper way.

Notes
-----
For full integrated slave element, coefficient matrix for Tri3 is
Ae = [3.0 -1.0 -1.0; -1.0 3.0 -1.0; -1.0 -1.0 3.0]

References
----------

[Popp2013] Popp, Alexander, et al. "Improved robustness and consistency of 3D contact algorithms based on a dual mortar approach." Computer Methods in Applied Mechanics and Engineering 264 (2013): 67-80.
"""
function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false) where E<:Union{Tri3, Quad4}

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    area = 0.0

    slave_element_nodes = get_connectivity(slave_element)
    nsl = length(slave_element)
    X1 = slave_element("geometry", time)
    n1 = slave_element("normal", time)

    # project slave nodes to auxiliary plane (x0, Q)
    xi = get_mean_xi(slave_element)
    N = vec(get_basis(slave_element, xi, time))
    x0 = interpolate(N, X1)
    n0 = interpolate(N, n1)
    S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]

    master_elements = slave_element("master elements", time)

    if props.dual_basis

        De = zeros(nsl, nsl)
        Me = zeros(nsl, nsl)

        for master_element in master_elements

            master_element_nodes = get_connectivity(master_element)
            nm = length(master_element)
            X2 = master_element("geometry", time)

            if norm(mean(X1) - mean(X2)) > problem.properties.distval
                continue
            end

            # 3.1 project master nodes to auxiliary plane and create polygon clipping
            M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]
            P = get_polygon_clip(S, M, n0)
            length(P) < 3 && continue # no clipping or shared edge (no volume)
            check_orientation!(P, n0)
            N_P = length(P)
            P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])

            if isapprox(P_area, 0.0)
                @info("Polygon P has zero area: $P_area")
                continue
            end

            # 4. loop integration cells
            C0 = calculate_centroid(P)
            all_cells = get_cells(P, C0)
            for cell in all_cells
                virtual_element = Element(Tri3, Int[])
                update!(virtual_element, "geometry", tuple(cell...))
                for ip in get_integration_points(virtual_element, 3)
                    detJ = virtual_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ
                    x_gauss = virtual_element("geometry", ip, time)
                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
                    N1 = slave_element(xi_s, time)
                    De += w*Matrix(Diagonal(vec(N1)))
                    Me += w*N1'*N1
                end
            end # integration cells done

        end # master elements done

        Ae = De*inv(Me)

        @info("Dual basis coefficient matrix: $Ae")

    else
        Ae = Matrix(1.0I, nsl, nsl)
    end

    for master_element in master_elements

        master_element_nodes = get_connectivity(master_element)
        nm = length(master_element)
        X2 = master_element("geometry", time)

        if norm(mean(X1) - mean(X2)) > problem.properties.distval
            continue
        end

        # 3.1 project master nodes to auxiliary plane and create polygon clipping
        M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]
        P = get_polygon_clip(S, M, n0)
        length(P) < 3 && continue # no clipping or shared edge (no volume)
        check_orientation!(P, n0)
        N_P = length(P)
        P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])

        if isapprox(P_area, 0.0)
            @info("Polygon P has zero area: $P_area")
            continue
        end

        C0 = calculate_centroid(P)

        De = zeros(nsl, nsl)
        Me = zeros(nsl, nm)
        ge = zeros(field_dim*nsl)

        # 4. loop integration cells
        all_cells = get_cells(P, C0)
        for cell in all_cells
            virtual_element = Element(Tri3, Int[])
            update!(virtual_element, "geometry", tuple(cell...))

            # 5. loop integration point of integration cell
            for ip in get_integration_points(virtual_element, 3)
                detJ = virtual_element(ip, time, Val{:detJ})
                w = ip.weight*detJ

                # project gauss point from auxiliary plane to master and slave element
                x_gauss = virtual_element("geometry", ip, time)

                xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
                xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)

                # add contributions
                N1 = vec(get_basis(slave_element, xi_s, time))
                N2 = vec(get_basis(master_element, xi_m, time))
                Phi = Ae*N1
                # Phi = [3.0-4.0*xi_s[1]-4.0*xi_s[2], 4.0*xi_s[1]-1.0, 4.0*xi_s[2]-1.0]
                De += w*Phi*N1'
                Me += w*Phi*N2'
                if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement")
                    u1 = slave_element("displacement", time)
                    u2 = master_element("displacement", time)
                    x_s = interpolate(N1, map(+,X1,u1))
                    x_m = interpolate(N2, map(+,X2,u2))
                    ge += w*vec((x_m-x_s)*Phi')
                end
                area += w
            end # integration points done

        end # integration cells done

        # 6. add contribution to contact virtual work
        sdofs = get_gdofs(problem, slave_element)
        mdofs = get_gdofs(problem, master_element)

        for i=1:field_dim
            lsdofs = sdofs[i:field_dim:end]
            lmdofs = mdofs[i:field_dim:end]
            add!(problem.assembly.C1, lsdofs, lsdofs, De)
            add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
            add!(problem.assembly.C2, lsdofs, lsdofs, De)
            add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
        end
        add!(problem.assembly.g, sdofs, ge)

    end # master elements done

    return area
end

""" Assemble quadratic surface element to problem.

In polygon clipping element is divided to linear sub-elements proposed in [Puso2008].

References
----------

[Puso2008] Puso, Michael A., T. A. Laursen, and Jerome Solberg. "A segment-to-segment mortar contact method for quadratic elements and large deformations." Computer Methods in Applied Mechanics and Engineering 197.6 (2008): 555-566.

[Popp1012] Popp, Alexander, et al. "Dual quadratic mortar finite element methods for 3D finite deformation contact." SIAM Journal on Scientific Computing 34.4 (2012): B421-B446.

"""
function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false) where E<:Union{Tri6}

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    area = 0.0

    Xs = slave_element("geometry", time)

    alp = props.alpha

    if alp != 0.0
        T = [
                1.0 0.0 0.0 0.0 0.0 0.0
                0.0 1.0 0.0 0.0 0.0 0.0
                0.0 0.0 1.0 0.0 0.0 0.0
                alp alp 0.0 1.0-2*alp 0.0 0.0
                0.0 alp alp 0.0 1.0-2*alp 0.0
                alp 0.0 alp 0.0 0.0 1.0-2*alp
            ]
    else
        T = Matrix(1.0I, 6, 6)
    end

    #=
    invT = [
        1.0 0.0 0.0 0.0 0.0 0.0
        0.0 1.0 0.0 0.0 0.0 0.0
        0.0 0.0 1.0 0.0 0.0 0.0
        -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
        0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
        -alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
    ]
    =#

    if props.dual_basis
        # @info("Creating dual basis for element $(slave_element.id)")
        nsl = length(slave_element)
        De = zeros(nsl, nsl)
        Me = zeros(nsl, nsl)

        # split slave element to linear sub-elements and loop
        for sub_slave_element in split_quadratic_element(slave_element, time)

            slave_element_nodes = get_connectivity(sub_slave_element)
            nsl = length(sub_slave_element)
            X1 = sub_slave_element("geometry", time)
            n1 = sub_slave_element("normal", time)

            # create auxiliary plane
            xi = get_mean_xi(sub_slave_element)
            N = vec(get_basis(sub_slave_element, xi, time))
            x0 = interpolate(N, X1)
            n0 = interpolate(N, n1)

            # project slave nodes to auxiliary plane
            S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]

            # 3. loop all master elements
            master_elements = slave_element("master elements", time)

            for master_element in master_elements

                Xm = master_element("geometry", time)

                if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
                    continue
                end

                # split master element to linear sub-elements and loop
                for sub_master_element in split_quadratic_element(master_element, time)

                    master_element_nodes = get_connectivity(sub_master_element)
                    nm = length(sub_master_element)
                    X2 = sub_master_element("geometry", time)

                    # 3.1 project master nodes to auxiliary plane
                    M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]

                    # create polygon clipping P
                    P = get_polygon_clip(S, M, n0)
                    length(P) < 3 && continue # no clipping or shared edge (no volume)
                    check_orientation!(P, n0)
                    N_P = length(P)
                    P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])

                    C0 = calculate_centroid(P)

                    # 4. loop integration cells
                    all_cells = get_cells(P, C0)
                    for cell in all_cells
                        virtual_element = Element(Tri3, Int[])
                        update!(virtual_element, "geometry", tuple(cell...))
                        for ip in get_integration_points(virtual_element, 3)
                            x_gauss = virtual_element("geometry", ip, time)
                            xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
                            detJ = virtual_element(ip, time, Val{:detJ})
                            w = ip.weight*detJ
                            N1 = vec(slave_element(xi_s, time)*T)
                            De += w*Matrix(Diagonal(N1))
                            Me += w*N1*N1'
                        end

                    end # integration cells done

                end # sub master elements done

            end # master elements done

        end # sub slave elements done

        Ae = De*inv(Me)
        # @info("Dual basis construction finished.")
        # @info("Slave element geometry = $Xs")
        # @info("De = $De")
        # @info("Me = $Me")
        # @info("Dual basis coefficient matrix: $Ae")

    else
        nsl = length(slave_element)
        Ae = Matrix(1.0I, nsl, nsl)
    end

    # split slave element to linear sub-elements and loop
    for sub_slave_element in split_quadratic_element(slave_element, time)

        slave_element_nodes = get_connectivity(sub_slave_element)
        nsl = length(sub_slave_element)
        X1 = sub_slave_element("geometry", time)
        n1 = sub_slave_element("normal", time)

        # create auxiliary plane
        xi = get_mean_xi(sub_slave_element)
        N = vec(get_basis(sub_slave_element, xi, time))
        x0 = interpolate(N, X1)
        n0 = interpolate(N, n1)

        # project slave nodes to auxiliary plane
        S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]

        # 3. loop all master elements
        master_elements = slave_element("master elements", time)

        for master_element in master_elements

            Xm = master_element("geometry", time)

            if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
                continue
            end

            # split master element to linear sub-elements and loop
            for sub_master_element in split_quadratic_element(master_element, time)

                master_element_nodes = get_connectivity(sub_master_element)
                nm = length(sub_master_element)
                X2 = sub_master_element("geometry", time)

                # 3.1 project master nodes to auxiliary plane
                M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]

                # create polygon clipping P
                P = get_polygon_clip(S, M, n0)
                length(P) < 3 && continue # no clipping or shared edge (no volume)
                check_orientation!(P, n0)
                N_P = length(P)
                P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])

                if isapprox(P_area, 0.0)
                    @warn("Polygon P has zero area: $P_area")
                    continue
                end

                C0 = calculate_centroid(P)

                # while our polygon clipping algorithm is working in linear sub elements
                # contributions is calculated using quadratic shape functions
                De = zeros(length(slave_element), length(slave_element))
                Me = zeros(length(slave_element), length(master_element))
                ge = zeros(field_dim*length(slave_element))

                # 4. loop integration cells
                all_cells = get_cells(P, C0)
                for cell in all_cells
                    virtual_element = Element(Tri3, Int[])
                    update!(virtual_element, "geometry", tuple(cell...))

                    # 5. loop integration point of integration cell
                    for ip in get_integration_points(virtual_element, 3)

                        x_gauss = virtual_element("geometry", ip, time)
                        xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
                        xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time)

                        # add contributions
                        N1 = vec(slave_element(xi_s, time)*T)
                        N2 = vec(master_element(xi_m, time))
                        Phi = Ae*N1

                        detJ = virtual_element(ip, time, Val{:detJ})
                        w = ip.weight*detJ

                        De += w*Phi*N1'
                        Me += w*Phi*N2'
                        if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement")
                            u1 = slave_element("displacement", time)
                            u2 = master_element("displacement", time)
                            xs = interpolate(N1, map(+,Xs,u1))
                            xm = interpolate(N2, map(+,Xm,u2))
                            ge += w*vec((xm-xs)*Phi')
                        end
                        area += w
                    end # integration points done

                end # integration cells done

                # 6. add contribution to contact virtual work
                sdofs = get_gdofs(problem, slave_element)
                mdofs = get_gdofs(problem, master_element)

                for i=1:field_dim
                    lsdofs = sdofs[i:field_dim:end]
                    lmdofs = mdofs[i:field_dim:end]
                    add!(problem.assembly.C1, lsdofs, lsdofs, De)
                    add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
                    add!(problem.assembly.C2, lsdofs, lsdofs, De)
                    add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
                end
                add!(problem.assembly.g, sdofs, ge)

            end # sub aster elements done

        end # master elements done

    end # sub slave elements done

    return area
end

function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}})

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)
    area = 0.0

    #=
    if props.split_quadratic_slave_elements
        if !props.linear_surface_elements
            @warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
        end
        slave_elements = split_quadratic_elements(slave_elements, time)
    end
    =#

    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
    normals = calculate_normals(slave_elements, time, Val{2};
                                rotate_normals=props.rotate_normals)

    update!(slave_elements, "normal", time => normals)

    # 2. loop all slave elements
    first_slave_element = true

    for slave_element in slave_elements

        area += assemble!(problem, slave_element, time; first_slave_element=first_slave_element)
        first_slave_element = false

    end # slave elements done, contact virtual work ready

    C1 = sparse(problem.assembly.C1)
    C2 = sparse(problem.assembly.C2)

    maxdim = maximum(size(C1))
    if problem.properties.alpha != 0.0
        alp = problem.properties.alpha
        Te = [
                1.0 0.0 0.0 0.0 0.0 0.0
                0.0 1.0 0.0 0.0 0.0 0.0
                0.0 0.0 1.0 0.0 0.0 0.0
                alp alp 0.0 1.0-2*alp 0.0 0.0
                0.0 alp alp 0.0 1.0-2*alp 0.0
                alp 0.0 alp 0.0 0.0 1.0-2*alp
            ]
        invTe = [
            1.0 0.0 0.0 0.0 0.0 0.0
            0.0 1.0 0.0 0.0 0.0 0.0
            0.0 0.0 1.0 0.0 0.0 0.0
            -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
            0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
            -alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
        ]
        # construct global transformation matrices T and invT
        T = SparseMatrixCOO()
        invT = SparseMatrixCOO()
        for element in slave_elements
            dofs = get_gdofs(problem, element)
            for i=1:field_dim
                ldofs = dofs[i:field_dim:end]
                add!(T, ldofs, ldofs, Te)
                add!(invT, ldofs, ldofs, invTe)
            end
        end
        T = sparse(T, maxdim, maxdim, (a, b) -> b)
        invT = sparse(invT, maxdim, maxdim, (a, b) -> b)
        # fill diagonal
        d = ones(size(T, 1))
        d[get_nonzero_rows(T)] .= 0.0
        T += sparse(Diagonal(d))
        invT += sparse(Diagonal(d))
        #invT2 = sparse(inv(full(T)))
        #@info("invT == invT2? ", invT == invT2)
        #maxabsdiff = maximum(abs(invT - invT2))
        #@info("max diff = $maxabsdiff")
        C1 = C1*invT
        C2 = C2*invT
    end

    tol = problem.properties.drop_tolerance
    SparseArrays.droptol!(C1, tol)
    SparseArrays.droptol!(C2, tol)

    problem.assembly.C1 = C1
    problem.assembly.C2 = C2

end


text:
"""
**Coefficient of variation**

Corrected for the sample size.

"""
function coefficient_of_variation(x)
    cv = std(x) / mean(x)
    norm = 1 + 1 / (4 * length(x))
    return norm * cv
end

"""
**Number of surviving species**

Number of species with a biomass larger than the `threshold`. The threshold is
by default set at `eps()`, which should be close to 10^-16.
"""
function species_richness(parameters; threshold::Float64=eps(), last::Int64=1000)
    @assert last <= size(parameters[:B], 1)
    measure_on = parameters[:B][end-(last-1):end,:]
    if sum(measure_on) == 0
        return NaN
    end
    richness = vec(sum(measure_on .> threshold, dims = 2))
    return mean(richness)
end

"""
**Proportion of surviving species**

Proportion of species with a biomass larger than the `threshold`. The threshold is
by default set at `eps()`, which should be close to 10^-16.
"""
function species_persistence(parameters; threshold::Float64=eps(), last::Int64=1000)
    r = species_richness(parameters, threshold=threshold, last=last)
    m = size(parameters[:B], 2) # Number of species is the number of columns in the biomass matrix
    return r/m
end

"""
**Population stability**

Population stability is measured as the mean of the negative coefficient
of variations of all species with an abundance higher than `threshold`. By
default, the stability is measured over the last `last=1000` timesteps.
"""
function population_stability(parameters; threshold::Float64=eps(), last=1000)
    @assert last <= size(parameters[:B], 1)
    non_extinct = parameters[:B][end,:] .> threshold
    measure_on = parameters[:B][end-(last-1):end,non_extinct]
    if sum(measure_on) == 0
        return NaN
    end
    stability = -mapslices(coefficient_of_variation, measure_on, dims = 1)
    return mean(stability)
end

"""
**Total biomass**

Returns the sum of biomass, averaged over the last `last` timesteps.

"""
function total_biomass(parameters; last=1000)
    @assert last <= size(parameters[:B], 1)
    measure_on = parameters[:B][end-(last-1):end,:]
    if sum(measure_on) == 0
        return NaN
    end
    biomass = vec(sum(measure_on, dims = 2))
    return mean(biomass)
end

"""
**Per species biomass**

Returns the average biomass of all species, over the last `last` timesteps.
"""
function population_biomass(parameters; last=1000)
    @assert last <= size(parameters[:B], 1)
    measure_on = parameters[:B][end-(last-1):end,:]
    if sum(measure_on) == 0
        return NaN
    end
    biomass = vec(mean(measure_on, dims = 1))
    return biomass
end

"""
**Shannon's entropy**

Corrected for the number of species, removes negative and null values, return
`NaN` in case of problem.
"""
function shannon(n)
    x = copy(n)
    x = filter((k) -> k > 0.0, x)
    try
        if length(x) > 1
            p = x ./ sum(x)
            corr = log.(length(x))
            p_ln_p = p .* log.(p)
            return -(sum(p_ln_p)/corr)
        else
            return NaN
        end
    catch
        return NaN
    end
end

"""
**Food web diversity**

Based on the average of Shannon's entropy (corrected for the number of
species) over the last `last` timesteps. Values close to 1 indicate that
all populations have equal biomasses.

"""
function foodweb_evenness(parameters; last=1000)
    @assert last <= size(parameters[:B], 1)
    measure_on = parameters[:B][end-(last-1):end,:]
    if sum(measure_on) == 0
        return NaN
    end
    shan = [shannon(vec(measure_on[i,:])) for i in 1:size(measure_on, 1)]
    return mean(shan)
end

# """
# **Save the output of a simulation**
#
# Takes a simulation output as a mandatory argument. The two keyword arguments
# are `as` (can be `:json` or `:jld`), defining the file format, and `filename`
# (without an extension, defaults to `NaN`). If `:jld` is used, the variable
# is named `befwm_simul` unless a `varname` is given.
#
# Called with the defaults, this function will write `befwm_xxxxxxxx.json`
# with the current simulation output, where `xxxxxxxx` is a hash of the `p`
# output (ensuring that all output files are unique).
#
# This function is *not* exported, so it must be called with `BioEnergeticFoodWebs.save`.
#
# """
# function save(parameters::Dict{Symbol,Any}; as::Symbol=:json, filename=nothing, varname=nothing)
#     if as == :JSON
#         as = :json
#     end
#     if as == :JLD
#         as = :jld
#     end
#     @assert as ∈ vec([:json :jld])
#     if filename == nothing
#         filename = "befwm_" * string(hash(parameters))
#     end
#     if varname == nothing
#         varname = "befwm_simul"
#     end
#     if as == :json
#         filename = filename * ".json"
#         f = open(filename, "w")
#         JSON.print(f, parameters)
#         close(f)
#     end
#     if as == :jld
#         filename = filename * ".jld"
#         JLD.save(filename, varname, parameters)
#     end
# end

"""
**Producers growth rate**

This function takes the simulation outputs from `simulate` and returns the producers
growth rates. Depending on the value given to the keyword `out_type`, it can return
more specifically:
- growth rates for each producer at each time step form end-last to last (`out_type = :all`)
- the mean growth rate for each producer over the last `last` time steps (`out_type = :mean`)
- the standard deviation of the growth rate for each producer over the last `last` time steps (`out_type = :std`)
"""
function producer_growth(out::Dict{Symbol,Any}; last::Int64 = 1000, out_type::Symbol = :all)
    parameters = out[:p] #extract parameters
    @assert last <= size(out[:B], 1)
    measure_on = out[:B][end-(last-1):end,:] #extract the biomasses that will be used
    measure_on_mat = [measure_on[i,:] for i = 1:last] #make it an array of array so we can use the map function
    if parameters[:productivity] == :nutrients #if the producers do NOT rely on nutrients for their growth
        c = out[:C][end-(last-1):end,:] #extract the timesteps of interest for the nutrients concentration
        c_mat = [c[i,:] for i = 1:last] #make it an array of array
        gr = map((x,y) -> get_growth(parameters,x,c=y), measure_on_mat, c_mat)
        growth = hcat(map(x -> x[1], gr)...)'
    else
        gr = map(x -> get_growth(parameters,x), measure_on_mat)
        growth = hcat(map(x -> x[1], gr)...)'
    end
    growth[:,.!parameters[:is_producer]] .= 0.0
    if out_type == :all #return all growth rates (each producer at each time step)
        return growth
    elseif out_type == :mean #return the producers mean growth rate over the last `last` time steps
        return mean(growth, dims = 1)
    elseif out_type == :std #return the growth rate standard deviation over the last `last` time steps (for each producer)
        return std(growth, dims = 1)
    else #if the keyword used is not one of :mean, :all or :std, print an error
        error("out_type should be one of :all, :mean or :std")
    end
end

"""
**Nutrients intake**

This function takes the simulation outputs from `simulate` and returns the producers
nutrient intake. Depending on the value given to the keyword `out_type`, it can return
more specifically:

- nutrient intake for each producer at each time step form end-last to last (`out_type = :all`)
- the mean nutrient intake for each producer over the last `last` time steps (`out_type = :mean`)
- the standard deviation of the nutrient intake for each producer over the last `last` time steps (`out_type = :std`)
"""
function nutrient_intake(out::Dict{Symbol,Any}; last::Int64 = 1000, out_type::Symbol = :all)
    parameters = out[:p] #extract parameters
    @assert last <= size(out[:B], 1)
    @assert parameters[:productivity] == :nutrients
    measure_on = out[:B][end-(last-1):end,:] #extract the biomasses that will be used
    measure_on_mat = [measure_on[i,:] for i = 1:last] #make it an array of array so we can use the map function
    c = out[:C][end-(last-1):end,:] #extract the timesteps of interest for the nutrients concentration
    c_mat = [c[i,:] for i = 1:last] #make it an array of array
    gr = map((x,y) -> get_growth(parameters,x,c=y), measure_on_mat, c_mat)
    intake = hcat(map(x -> x[2], gr)...)'
    if out_type == :all #return all growth rates (each producer at each time step)
        return intake
    elseif out_type == :mean #return the producers mean growth rate over the last `last` time steps
        return mean(intake, dims = 1)
    elseif out_type == :std #return the growth rate standard deviation over the last `last` time steps (for each producer)
        return std(intake, dims = 1)
    else #if the keyword used is not one of :mean, :all or :std, print an error
        error("out_type should be one of :all, :mean or :std")
    end
end

"""
**Consumers' biomass intake**

This function takes the simulation outputs from `simulate` and returns the consumers
biomass intake. Depending on the value given to the keyword `out_type`, it can return
more specifically:

- biomass intake for each species at each time step form end-last to last (`out_type = :all`)
- the mean biomass intake for each species over the last `last` time steps (`out_type = :mean`)
- the standard deviation of the biomass intake for each species over the last `last` time steps (`out_type = :std`)

"""
function consumer_intake(out::Dict{Symbol,Any}; last::Int64 = 1000, out_type::Symbol = :all)
    parameters = out[:p] #extract parameters
    @assert last <= size(out[:B], 1)
    measure_on = out[:B][end-(last-1):end,:] #extract the biomasses that will be used
    measure_on_mat = [measure_on[i,:] for i = 1:last] #make it an array of array so we can use the map function
    cons = map(x -> consumption(parameters, x), measure_on_mat)
    gains = hcat(map(x-> x[1], cons)...)'
    #losses = hcat(map(x-> x[:loss], cons)...)'
    if out_type == :all #return all growth rates (each producer at each time step)
        return gains
    elseif out_type == :mean #return the producers mean growth rate over the last `last` time steps
        return mean(gains, dims = 1)
    elseif out_type == :std #return the growth rate standard deviation over the last `last` time steps (for each producer)
        return std(gains, dims = 1)
    else #if the keyword used is not one of :mean, :all or :std, print an error
        error("out_type should be one of :all, :mean or :std")
    end
end

"""
**Metabolic loss**

This function takes the simulation outputs from `simulate` and returns the species
metabolic losses. Depending on the value given to the keyword `out_type`, it can return
more specifically:

- metabolic losses for each species at each time step form end-last to last (`out_type = :all`)
- the mean metabolic loss for each species over the last `last` time steps (`out_type = :mean`)
- the standard deviation of the metabolic losses for each species over the last `last` time steps (`out_type = :std`)
"""
function metabolism(out::Dict{Symbol,Any}; last::Int64 = 1000, out_type::Symbol = :all)
    parameters = out[:p] #extract parameters
    @assert last <= size(out[:B], 1)
    measure_on = out[:B][end-(last-1):end,:] #extract the biomasses that will be used
    measure_on_mat = [measure_on[i,:] for i = 1:last] #make it an array of array so we can use the map function
    function metab(b,p)
        if parameters[:productivity] == :nutrients
            m = parameters[:x] .* b
        else
            m = (parameters[:x] .* b) .* .!parameters[:is_producer]
        end
        return m
    end
    metabolic_losses = hcat(map(x -> metab(x, parameters), measure_on_mat)...)'
    if out_type == :all #return all growth rates (each producer at each time step)
        return metabolic_losses
    elseif out_type == :mean #return the producers mean growth rate over the last `last` time steps
        return mean(metabolic_losses, dims = 1)
    elseif out_type == :std #return the growth rate standard deviation over the last `last` time steps (for each producer)
        return std(metabolic_losses, dims = 1)
    else #if the keyword used is not one of :mean, :all or :std, print an error
        error("out_type should be one of :all, :mean or :std")
    end
end

ext: jl

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
# Command: Find Elements from Element
"""
    Elements(element::Element, location_strategy::AbstractString, value::AbstractString)::Vector{Element}

Find Elements

See also: [`Element`](@ref)
"""
function Elements(
    element::Element,
    location_strategy::AbstractString,
    value::AbstractString,
)::Vector{Element}
    location_strategy ∈
    ["css selector", "link text", "partial link text", "tag name", "xpath"] ||
    throw(ArgumentError("""location_strategy must be "css selector", "link text", "partial link text", "tag name" or "xpath"."""))
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.post(
        "$addr/session/$id/element/$element_id/elements",
        [("Content-Type" => "application/json")],
        JSON3.write(Dict("using" => location_strategy, "value" => value)),
    )
    @assert response.status == 200
    Element.(element, first(values(node)) for node in JSON3.read(response.body).value)
end


text:
module NeuralNetwork

using LinearAlgebra
using NNlib
using Zygote
using Zygote: @adjoint

export NetWork, Dense, Conv, Flatten, Dropout, MaxPool, add_layer!, train!, @epochs, params, KeepSize

include("utils.jl")
include("activations.jl")
include("Layers/Network.jl")
include("Optimize/train.jl")

end


text:
using Graphs
using SimpleGraphs

function save_json(arg,fname::ASCIIString,gname::ASCIIString)
    d = JSON.json(arg)
    f = open(@sprintf("%s",joinpath(fname,gname)),"w")
    JSON.print(f,d)
    close(f)
end

function plot_to_file(g::GenericGraph,fname::ASCIIString,gname::ASCIIString)
    write_to_file(g,fname,gname)
    pname = joinpath(fname,gname)
    run(`neato -Tpng $pname".dot" -o $pname".png"`)
end

function write_to_file(g::GenericGraph,fname::ASCIIString,gname::ASCIIString)
    if !ispath(fname)
        mkdir(fname)
    end
    f = open(@sprintf("%s.dot",joinpath(fname,gname)),"w")
    write(f,to_dot(g))
    close(f)
end

function plot_to_file(g::SimpleGraphs.SimpleGraph,fname::ASCIIString,gname::ASCIIString)
    (G,d,dinv) = convert_simple(g)
    plot_to_file(G,fname,gname)
end

function write_to_file(g::SimpleGraphs.SimpleGraph,fname::ASCIIString,gname::ASCIIString)
    (G,d,dinv) = convert_simple(g)
    write_to_file(G,fname,gname)
end

function simple_NE_sgraph(n::Int,e::Int)
    g = SimpleGraphs.SimpleGraph()
    for i = 1:e
        pts = rand(collect(1:n),2)
        add!(g,pts[1],pts[2])
    end
    if NE(g) != e
        # add in missing edges
        for i = 1:(e-NE(g))
            ns = nt = 1
            while ns == nt || has(g,ns,nt)
                ns = rand(1:NV(g))
                nt = rand(1:NV(g))
            end
            add!(g,ns,nt)
        end
    end
    return g
end

function simple_NE_graph(N,E;parallel_edges=false,self_loops=false)
    g = simple_graph(N,is_directed=false)
    for i = 1:E
        if !parallel_edges
            edge_placed = false
            while !edge_placed
                pts = rand(collect(1:N),2)
                if !(pts[2] in collect(Graphs.out_neighbors(pts[1],g)))
                    if !self_loops && (pts[2] != pts[1])
                        add_edge!(g,pts[1],pts[2])
                        edge_placed = true
                    elseif self_loops
                        add_edge!(g,pts[1],pts[2])
                        edge_placed = true
                    end
                else
                    continue
                end
            end
        else
            pts = rand(collect(1:N),2)
            add_edge!(g,pts[1],pts[2])
        end
    end
    return g
end


text:
export PrioritizedSweepingSamplingModel

using DataStructures: PriorityQueue, dequeue!

import StatsBase: sample

"""
    PrioritizedSweepingSamplingModel(θ::Float64=1e-4)
See more details at Section (8.4) on Page 168 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
mutable struct PrioritizedSweepingSamplingModel <: AbstractEnvironmentModel
    experiences::Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}
    PQueue::PriorityQueue{Tuple{Any,Any},Float64}
    predecessors::Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}
    θ::Float64
    sample_count::Int
    PrioritizedSweepingSamplingModel(θ::Float64 = 1e-4) = new(
        Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}(),
        PriorityQueue{Tuple{Any,Any},Float64}(Base.Order.Reverse),
        Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}(),
        θ,
        0,
    )
end

function RLBase.update!(
    m::PrioritizedSweepingSamplingModel,
    t::AbstractTrajectory,
    p::AbstractPolicy,
    ::AbstractEnv,
    ::Union{PreActStage,PostEpisodeStage},
)
    if length(t[:terminal]) > 0
        transition = (
            t[:state][end-1],
            t[:action][end-1],
            t[:reward][end],
            t[:terminal][end],
            t[:state][end],
        )
        pri = RLBase.priority(p, transition)
        update!(m, (transition..., pri))
    end
end

function RLBase.update!(m::PrioritizedSweepingSamplingModel, transition::Tuple)
    s, a, r, d, s′, P = transition
    m.experiences[(s, a)] = (r, d, s′)
    if P >= m.θ
        m.PQueue[(s, a)] = P
    end
    if !haskey(m.predecessors, s′)
        m.predecessors[s′] = Set{Tuple{Any,Any,Float64,Bool}}()
    end
    push!(m.predecessors[s′], (s, a, r, d))
end

function sample(m::PrioritizedSweepingSamplingModel)
    if length(m.PQueue) > 0
        s, a = dequeue!(m.PQueue)
        r, d, s′ = m.experiences[(s, a)]
        m.sample_count += 1
        s, a, r, d, s′
    else
        nothing
    end
end


text:
# Deterministic equivalent generation #
# ========================== #
function generate!(stochasticprogram::StochasticProgram{N}, structure::DeterministicEquivalent{N}) where N
    # Set the optimizer
    structure.model.moi_backend = optimizer(stochasticprogram)
    # Prepare decisions
    structure.model.ext[:decisions] = structure.decisions
    add_decision_bridges!(structure.model)
    # Generate all stages
    for stage in 1:N
        generate!(stochasticprogram, structure, stage)
        # Register stage decisions
        if stage < N
            append!(structure.decision_variables[stage], all_decision_variables(structure.model))
            for previous_stage in structure.decision_variables[1:stage-1]
                filter!(d -> d ∈ previous_stage, structure.decision_variables[stage])
            end
        end
    end
    return nothing
end

function generate!(stochasticprogram::StochasticProgram{N}, structure::DeterministicEquivalent{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    dep_model = structure.model
    if stage == 1
        # Define first-stage problem
        generator(stochasticprogram, :stage_1)(dep_model, stage_parameters(stochasticprogram, 1))
    else
        # Sanity check on scenario probabilities
        if num_scenarios(stochasticprogram, stage) > 0
            p = stage_probability(stochasticprogram, stage)
            abs(p - 1.0) <= 1e-6 || @warn "Scenario probabilities do not add up to one. The probability sum is given by $p"
        end
        # Cache current objective
        dep_obj = objective_function(dep_model)
        obj_sense = objective_sense(dep_model)
        obj_sense = obj_sense == MOI.FEASIBILITY_SENSE ? MOI.MIN_SENSE : obj_sense
        # Define second-stage problems, renaming variables according to scenario.
        stage_two_params = stage_parameters(stochasticprogram, 2)
        visited_objs = collect(keys(object_dictionary(dep_model)))
        for (i, scenario) in enumerate(scenarios(stochasticprogram))
            stage_key = Symbol(:stage_, stage)
            generator(stochasticprogram, stage_key)(dep_model, stage_parameters(stochasticprogram, stage), scenario)
            sub_sense = objective_sense(dep_model)
            sub_obj = objective_function(dep_model)
            if obj_sense == sub_sense
                dep_obj += probability(scenario) * objective_function(dep_model)
                push!(structure.sub_objectives[stage - 1], moi_function(sub_obj))
            else
                dep_obj -= probability(scenario) * objective_function(dep_model)
                push!(structure.sub_objectives[stage - 1], -moi_function(sub_obj))
            end
            for (objkey,obj) ∈ filter(kv->kv.first ∉ visited_objs, object_dictionary(dep_model))
                newkey = if isa(obj, VariableRef) || isa(obj, DecisionRef)
                    varname = if N > 2
                        varname = add_subscript(add_subscript(JuMP.name(obj), stage), i)
                    else
                        varname = add_subscript(JuMP.name(obj), i)
                    end
                    set_name(obj, varname)
                    newkey = Symbol(varname)
                elseif isa(obj, AbstractArray{<:VariableRef}) || isa(obj, AbstractArray{<:DecisionRef})
                    arrayname = if N > 2
                        arrayname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        arrayname = add_subscript(objkey, i)
                    end
                    for var in obj
                        splitname = split(name(var), "[")
                        varname = if N > 2
                            varname = add_subscript(splitname[1], stage)
                        else
                            varname = splitname[1]
                        end
                        set_name(var, @sprintf("%s[%s", add_subscript(varname,i), splitname[2]))
                    end
                    newkey = Symbol(arrayname)
                elseif isa(obj, JuMP.ConstraintRef) ||
                       isa(obj, AbstractArray{<:ConstraintRef}) ||
                       isa(obj, AbstractArray{<:GenericAffExpr}) ||
                       isa(obj, AbstractArray{<:DecisionAffExpr})
                    arrayname = if N > 2
                        arrayname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        arrayname = add_subscript(objkey, i)
                    end
                    newkey = Symbol(arrayname)
                elseif isa(obj, AbstractJuMPScalar)
                    exprname = if N > 2
                        exprname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        exprname = add_subscript(objkey, i)
                    end
                    newkey = Symbol(exprname)
                else
                    continue
                end
                dep_model.obj_dict[newkey] = obj
                delete!(dep_model.obj_dict, objkey)
                push!(visited_objs, newkey)
            end
        end
        set_objective_function(dep_model, dep_obj)
        set_objective_sense(dep_model, obj_sense)
    end
    return nothing
end

function clear!(dep::DeterministicEquivalent)
    # Clear deterministic equivalent model
    empty!(dep.model)
    return nothing
end


text:
using Test
using NURBS

@testset "NURBS.jl" begin
    @testset "BSplineBasis.jl" begin
        include("BSplineBasis.jl")
    end
    @testset "NURBSBasis.jl" begin
        include("NURBSBasis.jl")
    end
end


text:
# We will be wrapping MLJ machines/resampling machines
# (cf. MLJBase/src/resampling.jl). It is for the wrapped type that we
# will be overloading the methods of `IterationControl`.

const ERR_EVALUATION =
    ArgumentError("There are no evaluation objects if `resampling=nothing`. ")

mlj_model(mach::Machine) = mach.model
mlj_model(mach::Machine{<:Resampler}) = mach.model.model

# training machine wrappper (constructor resets mlj model iterations
# to zero):
struct ICModel{M}
    mach::M
    iteration_parameter::Union{Symbol,Expr}
    verbosity::Integer
    n::Base.RefValue{Int64} # num calls to train!
    i::Base.RefValue{Int64} # num iterations
    Δi::Base.RefValue{Int64}
    function ICModel(mach::M, iteration_parameter, verbosity) where M
        model = mlj_model(mach)
        rset!(model, iteration_parameter, 0)
        return new{M}(mach,
                      iteration_parameter,
                      verbosity,
                      Ref(0),
                      Ref(0),
                      Ref(0))
    end
end

mlj_model(ic_model::ICModel) = mlj_model(ic_model.mach)

# overloading `train!`:
function IterationControl.train!(m::ICModel, n::Int)

    # increment training call counter:
    m.n[] = m.n[] + 1

    model = mlj_model(m)

    # update iteration parameter value stored on ic_model:
    m.i[] = rget(model, m.iteration_parameter) + n

    # update the iteration parameter in the mlj model:
    rset!(model, m.iteration_parameter, m.i[])

    fit!(m.mach, verbosity=m.verbosity - 2)
    m.Δi[] = n
end

# overloading `expose`- for `resampling === nothing`:
IterationControl.expose(ic_model::ICModel) = ic_model.mach

# overloading `expose`- for `resampling isa Holdout` or
# other resampling strategy:
IterationControl.expose(ic_model::ICModel{<:Machine{<:Resampler}}) =
    MLJBase.fitted_params(ic_model.mach).machine

# overloading `loss` - for `resampling === nothing`:
function IterationControl.loss(m::ICModel)
    losses = training_losses(IterationControl.expose(m))
    losses isa Nothing && return nothing
    return last(losses)
end

# overloading `loss` - for `resampling isa Holdout`:
IterationControl.loss(m::ICModel{<:Machine{<:Resampler}}) =
    evaluate(m.mach).measurement |> first

# overloading `training_loss`:
function IterationControl.training_losses(m::ICModel)
    mach = IterationControl.expose(m)
    losses = training_losses(mach)
    losses isa Nothing && return nothing
    s = length(losses)
    return view(losses, (s - m.Δi[] + 1):s)
end

_evaluation(ic_model) = throw(ERR_EVALUATION)
_evaluation(ic_model::ICModel{<:Machine{<:Resampler}}) =
    MLJBase.evaluate(ic_model.mach)

# interface for user-defined controls:
function Base.getproperty(m::ICModel, property::Symbol)
    property === :machine && return IterationControl.expose(m)
    property === :model && return mlj_model(m)
    property === :n_cycles && return m.n[]
    property === :n_iterations && return m.i[]
    property === :Δiterations && return m.Δi[]
    property === :loss && return IterationControl.loss(m)
    property === :training_losses && return IterationControl.training_losses(m)
    property === :evaluation && return _evaluation(m)
    return getfield(m, property)
end
Base.propertynames(::ICModel) = (:machine,
                                 :model,
                                 :n_cycles,
                                 :n_iterations,
                                 :Δiterations,
                                 :loss,
                                 :training_losses,
                                 :evaluation)

max_stars_repo_name: ablaom/MLJIteration

max_issues_repo_name: ablaom/MLJIteration

max_forks_repo_name: ablaom/MLJIteration


text:
using BlockArrays, DataFrames, LinearAlgebra, MixedModels, NamedArrays
using Random, RData, SparseArrays, Statistics, Tables, Test

if !@isdefined(dat) || !isa(dat, Dict{Symbol, DataFrame})
    const dat = Dict(Symbol(k) => v for (k, v) in
        load(joinpath(dirname(pathof(MixedModels)), "..", "test", "dat.rda")))
end

const LMM = LinearMixedModel

@testset "Dyestuff" begin
    fm1 = LMM(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff])

    @test BlockArrays.nblocks(fm1.A) == (3, 3)
    @test size(fm1.reterms) == (1, )
    @test BlockArrays.nblocks(fm1.L) == (3, 3)
    @test lowerbd(fm1) == zeros(1)
    @test fm1.lowerbd == zeros(1)
    @test fm1.θ == ones(1)
    fm1.θ = ones(1)
    @test fm1.θ == ones(1)
#    @test_warn "Model has not been fit" show(fm1)

    @test objective(updateL!(setθ!(fm1, [0.713]))) ≈ 327.34216280955366
    MixedModels.describeblocks(IOBuffer(), fm1)

    fit!(fm1);
    @test :θ in propertynames(fm1)
    @test objective(fm1) ≈ 327.3270598811428 atol=0.001
    @test fm1.θ ≈ [0.752580] atol=1.e-5
    @test fm1.λ ≈ [LowerTriangular(reshape(fm1.θ, (1,1)))] atol=1.e-5
    @test deviance(fm1) ≈ 327.32705988 atol=0.001
    @test aic(fm1) ≈ 333.3270598811394 atol=0.001
    @test bic(fm1) ≈ 337.5306520261259 atol=0.001
    @test fixef(fm1) ≈ [1527.5]
    @test first(first(fm1.σs)) ≈ 37.26034462135931 atol=0.0001
    @test fm1.β ≈ [1527.5]
    @test dof(fm1) == 3
    @test nobs(fm1) == 30
    @test MixedModels.fixef!(zeros(1),fm1) ≈ [1527.5]
    @test coef(fm1) ≈ [1527.5]
    fm1β = fm1.βs
    @test fm1β isa NamedTuple
    @test isone(length(fm1β))
    @test first(values(fm1β)) ≈ 1527.5
    fm1σρ = fm1.σρs
    @test fm1σρ isa NamedTuple
    @test isone(length(fm1σρ))
    @test isone(length(getproperty(first(fm1σρ), :σ)))
    @test isempty(getproperty(first(fm1σρ), :ρ))
    @test fm1.σ == sdest(fm1)
    @test fm1.b == ranef(fm1)
    @test fm1.u == ranef(fm1, uscale=true)
    @test fm1.stderror == stderror(fm1)
    @test isone(length(fm1.pvalues))
    @test fm1.objective == objective(fm1)
    @test fm1.σ ≈ 49.510099986291145 atol=1.e-5
    @test fm1.X == ones(30,1)
    @test fm1.y == dat[:Dyestuff][!, :Y]
    @test cond(fm1) == ones(1)
    cm = coeftable(fm1)
    @test length(cm.rownms) == 1
    @test length(cm.colnms) == 4
    @test fnames(fm1) == (:G,)
    @test response(fm1) == dat[:Dyestuff][!, :Y]
    rfu = ranef(fm1, uscale = true)
    rfb = ranef(fm1)
    @test abs(sum(rfu[1])) < 1.e-5
    cv = condVar(fm1)
    @test length(cv) == 1
    @test size(first(cv)) == (1, 1, 6)
    show(IOBuffer(), fm1.optsum)

    @test logdet(fm1) ≈ 8.06014522999825 atol=0.001
    @test varest(fm1) ≈ 2451.2501089607676 atol=0.001
    @test pwrss(fm1) ≈ 73537.49947885796 atol=0.001
    @test stderror(fm1) ≈ [17.69455188898009] atol=0.0001

    vc = VarCorr(fm1)
    show(IOBuffer(), vc)
    @test vc.s == sdest(fm1)

    fit!(fm1, REML=true)
    @test objective(fm1) ≈ 319.65427684225216 atol=0.0001
    @test_throws ArgumentError loglikelihood(fm1)
    @test dof_residual(fm1) ≥ 0
    print(IOBuffer(), fm1)
end

@testset "Dyestuff2" begin
    fm = fit(MixedModel, @formula(Y ~ 1 + (1 | G)), dat[:Dyestuff2])
    @test lowerbd(fm) == zeros(1)
    show(IOBuffer(), fm)
    @test fm.θ ≈ zeros(1)
    @test objective(fm) ≈ 162.87303665382575
    @test abs(std(fm)[1][1]) < 1.0e-9
    @test std(fm)[2] ≈ [3.653231351374652]
    @test stderror(fm) ≈ [0.6669857396443261]
    @test coef(fm) ≈ [5.6656]
    @test logdet(fm) ≈ 0.0
    refit!(fm, dat[:Dyestuff][!, :Y])
    @test objective(fm) ≈ 327.3270598811428 atol=0.001
end

@testset "penicillin" begin
    fm = LMM(@formula(Y ~ 1 + (1 | G) + (1 | H)), dat[:Penicillin]);
    @test size(fm) == (144, 1, 30, 2)
    @test fm.θ == ones(2)
    @test lowerbd(fm) == zeros(2)

    fit!(fm)

    @test objective(fm) ≈ 332.18834867227616 atol=0.001
    @test coef(fm) ≈ [22.97222222222222] atol=0.001
    @test fixef(fm) ≈ [22.97222222222222] atol=0.001
    @test coef(fm)[1] ≈ mean(dat[:Penicillin][!, :Y])
    @test stderror(fm) ≈ [0.7445960346851368] atol=0.0001
    @test fm.θ ≈ [1.5375772376554968, 3.219751321180035] atol=0.001
    @test first(std(fm)) ≈ [0.8455645948223015] atol=0.0001
    @test std(fm)[2] ≈ [1.770647779277388] atol=0.0001
    @test varest(fm) ≈ 0.3024263987592062 atol=0.0001
    @test logdet(fm) ≈ 95.74614821367786 atol=0.001
    rfu = ranef(fm, uscale=true)
    rfb = ranef(fm)
    @test length(rfb) == 2
end

@testset "pastes" begin
    fm = LMM(@formula(Y ~ (1 | G) + (1 | H)), dat[:Pastes])
    @test size(fm) == (60, 1, 40, 2)
    @test fm.θ == ones(2)
    @test lowerbd(fm) == zeros(2)

    fit!(fm);

    @test objective(fm) ≈ 247.99446586289676 atol=0.001
    @test coef(fm) ≈ [60.05333333333329] atol=0.001
    @test fixef(fm) ≈ [60.05333333333329] atol=0.001
    @test stderror(fm) ≈ [0.6421359883527029] atol=0.0001
    @test fm.θ ≈ [3.5268858714382905, 1.3299230213750168] atol=0.001
    @test first(std(fm)) ≈ [2.904069002535747] atol=0.001
    @test std(fm)[2] ≈ [1.095070371687089] atol=0.0001
    @test std(fm)[3] ≈ [0.8234088395243269] atol=0.0001
    @test varest(fm) ≈ 0.6780020742644107 atol=0.0001
    @test logdet(fm) ≈ 101.0381339953986 atol=0.001
end

@testset "InstEval" begin
    fm1 = LMM(@formula(Y ~ 1 + A + (1 | G) + (1 | H) + (1 | I)), dat[:InstEval])
    @test size(fm1) == (73421, 2, 4114, 3)
    @test fm1.θ == ones(3)
    @test lowerbd(fm1) == zeros(3)

    fit!(fm1);

    @test objective(fm1) ≈ 237721.7687745563 atol=0.001
    ftd1 = fitted(fm1);
    @test size(ftd1) == (73421, )
    @test ftd1 == predict(fm1)
    @test first(ftd1) ≈ 3.17876 atol=0.0001
    resid1 = residuals(fm1);
    @test size(resid1) == (73421, )
    @test first(resid1) ≈ 1.82124 atol=0.00001
    @test length(fm1.rePCA) == 3

    fm2 = fit!(LinearMixedModel(@formula(Y ~ 1 + A*I + (1 | G) + (1 | H)), dat[:InstEval]))
    @test objective(fm2) ≈ 237585.5534151694 atol=0.001
    @test size(fm2) == (73421, 28, 4100, 2)
end

@testset "sleep" begin
    fm = LinearMixedModel(@formula(Y ~ 1 + U + (1 + U | G)), dat[:sleepstudy]);
    @test lowerbd(fm) == [0.0, -Inf, 0.0]
    A11 = fm.A[Block(1,1)]
    @test isa(A11, UniformBlockDiagonal{Float64})
    @test isa(fm.L[Block(1, 1)], UniformBlockDiagonal{Float64})
    @test size(A11) == (36, 36)
    @test A11.facevec[1] == [10. 45.; 45. 285.]
    @test length(A11.facevec) == 18
    updateL!(fm);
    b11 = LowerTriangular(fm.L[Block(1, 1)].facevec[1])
    @test b11 * b11' == fm.A[Block(1, 1)].facevec[1] + I
    @test count(!iszero, Matrix(fm.L[Block(1, 1)])) == 18 * 4
    @test rank(fm) == 2

    fit!(fm)

    @test objective(fm) ≈ 1751.9393444647046
    @test fm.θ ≈ [0.929221307, 0.01816838, 0.22264487096] atol=1.e-6
    @test pwrss(fm) ≈ 117889.46144025437
    @test logdet(fm) ≈ 73.90322021999222 atol=0.001
    @test stderror(fm) ≈ [6.632257721914501, 1.5022354739749826] atol=0.0001
    @test coef(fm) ≈ [251.40510484848477,10.4672859595959]
    @test fixef(fm) ≈  [251.40510484848477,10.4672859595959]
    @test std(fm)[1] ≈ [23.780468100188497, 5.716827903196682] atol=0.01
    @test logdet(fm) ≈ 73.90337187545992 atol=0.001
    @test cond(fm) ≈ [4.175251] atol=0.0001
    @test loglikelihood(fm) ≈ -875.9696722323523
    σs = fm.σs
    @test length(σs) == 1
    @test keys(σs) == (:G,)
    @test length(σs.G) == 2
    @test first(values(σs.G)) ≈ 23.780468626896745 atol=0.0001
    @test last(values(first(σs))) ≈ 5.716827808126002 atol=0.0001
    show(IOBuffer(), fm)

    u3 = ranef(fm, uscale=true)
    @test length(u3) == 1
    @test size(first(u3)) == (2, 18)
    @test first(u3)[1, 1] ≈ 3.030300122575336 atol=0.001
    u3n = first(ranef(fm, uscale=true, named=true))
    @test u3n isa NamedArrays.NamedArray

    b3 = ranef(fm)
    @test length(b3) == 1
    @test size(first(b3)) == (2, 18)
    @test first(first(b3)) ≈ 2.815819441982976 atol=0.001

    simulate!(fm)  # to test one of the unscaledre methods

    fmnc = zerocorr!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]));
    @test size(fmnc) == (180,2,36,1)
    @test fmnc.θ == ones(2)
    @test lowerbd(fmnc) == zeros(2)

    fit!(fmnc)

    @test deviance(fmnc) ≈ 1752.0032551398835 atol=0.001
    @test objective(fmnc) ≈ 1752.0032551398835 atol=0.001
    @test coef(fmnc) ≈ [251.40510484848585, 10.467285959595715]
    @test fixef(fmnc) ≈ [251.40510484848477, 10.467285959595715]
    @test stderror(fmnc) ≈ [6.707710260366577, 1.5193083237479683] atol=0.001
    @test fmnc.θ ≈ [0.9458106880922268, 0.22692826607677266] atol=0.0001
    @test first(std(fmnc)) ≈ [24.171449463289047, 5.799379721123582]
    @test last(std(fmnc)) ≈ [25.556130034081047]
    @test logdet(fmnc) ≈ 74.46952585564611 atol=0.001
    ρ = first(fmnc.σρs.G.ρ)
    @test ρ === -0.0   # test that systematic zero correlations are returned as -0.0

    fmnc2 = LinearMixedModel(@formula(Y ~ 1 + U + zerocorr(1+U|G)),
                             dat[:sleepstudy])
    @test size(fmnc2) == (180,2,36,1)
    @test fmnc2.θ == ones(2)
    @test lowerbd(fmnc2) == zeros(2)

    fit!(fmnc2)

    @test deviance(fmnc2) ≈ 1752.0032551398835 atol=0.001
    @test objective(fmnc2) ≈ 1752.0032551398835 atol=0.001
    @test coef(fmnc2) ≈ [251.40510484848585, 10.467285959595715]
    @test fixef(fmnc2) ≈ [251.40510484848477, 10.467285959595715]
    @test stderror(fmnc2) ≈ [6.707710260366577, 1.5193083237479683] atol=0.001
    @test fmnc2.θ ≈ [0.9458106880922268, 0.22692826607677266] atol=0.0001
    @test first(std(fmnc2)) ≈ [24.171449463289047, 5.799379721123582]
    @test last(std(fmnc2)) ≈ [25.556130034081047]
    @test logdet(fmnc2) ≈ 74.46952585564611 atol=0.001

#    MixedModels.lrt(fm, fmnc)

    MixedModels.likelihoodratiotest(fm, fmnc)

    fmrs = LinearMixedModel(@formula(Y ~ 1 + U + (0 + U|G)), dat[:sleepstudy]);
    fit!(fmrs)
    @test objective(fmrs) ≈ 1774.080315280528 rtol=0.00001
    @test fmrs.θ ≈ [0.24353985679033105] rtol=0.00001

    fm_ind = LinearMixedModel(@formula(Y ~ 1 + U + (1|G) + (0 + U|G)), dat[:sleepstudy]);
    fit!(fm_ind)
    @test objective(fm_ind) ≈ objective(fmnc)
    @test coef(fm_ind) ≈ coef(fmnc)
    @test fixef(fm_ind) ≈ fixef(fmnc)
    @test stderror(fm_ind) ≈ stderror(fmnc)
    @test fm_ind.θ ≈ fmnc.θ
    @test std(fm_ind) ≈ std(fmnc)
    @test logdet(fm_ind) ≈ logdet(fmnc)

    # combining [ReMat{T,S1}, ReMat{T,S2}] for S1 ≠ S2
    slpcat = categorical!(deepcopy(dat[:sleepstudy]), [:U])
    fm_cat = LMM(@formula(Y ~ 1 + U + (1|G) + (0 + U|G)), slpcat)
    @test fm_cat isa LMM
    σρ = fit!(fm_cat).σρs
    @test σρ isa NamedTuple
    @test isone(length(σρ))
    @test first(keys(σρ)) == :G
    @test keys(σρ.G) == (:σ, :ρ)
    @test length(σρ.G) == 2
    @test length(first(σρ.G)) == 11
    @test length(σρ.G.ρ) == 55
    @test iszero(σρ.G.ρ[46])
    @test σρ.G.ρ[46] === -0.0
end

@testset "d3" begin
    fm = updateL!(LMM(@formula(Y ~ 1 + U + (1+U|G) + (1+U|H) + (1+U|I)), dat[:d3]));
    @test pwrss(fm) ≈ 5.1261847180180885e6 rtol = 1e-6
    @test objective(fm) ≈ 901641.2930413672 rtol = 1e-6
    fit!(fm)
    @test objective(fm) ≈ 884957.5540213 rtol = 1e-6
    @test coef(fm) ≈ [0.4991229873, 0.31130780953] atol = 1.e-4
    bv = ranef(fm)
    @test length(bv) == 3
end

@testset "simulate!" begin
    fm = fit(MixedModel, @formula(Y ~ 1 + (1 | G)), dat[:Dyestuff])
    refit!(simulate!(Random.MersenneTwister(1234321), fm))
    @test deviance(fm) ≈ 339.0218639362958 atol=0.001
    refit!(fm, dat[:Dyestuff][!, :Y])
    Random.seed!(1234321)
    refit!(simulate!(fm))
    @test deviance(fm) ≈ 339.0218639362958 atol=0.001
    simulate!(fm, θ = fm.θ)
    @test_throws DimensionMismatch refit!(fm, zeros(29))
    bsamp = parametricbootstrap(MersenneTwister(1234321), 10, fm)
    @test length(bsamp.objective) == 10
    @test keys(bsamp.bstr) == (:objective, :σ, :β₁, :θ)
end

@testset "Rank deficient" begin
    Random.seed!(0)
    x = rand(100)
    data = columntable((x = x, x2 = 1.5 .* x, y = rand(100), z = repeat('A':'T', 5)))
    model = fit!(LinearMixedModel(@formula(y ~ x + x2 + (1|z)), data))
    @test length(fixef(model)) == 2
    @test rank(model) == 2
    @test length(coef(model)) == 3
end


text:
using SafeTestsets

@safetestset "solving 1 FDE" begin include("1FDE_test.jl") end

@safetestset "solving 2 parametric FDEs" begin include("2FDE_par_test.jl") end

@safetestset "solving 3 parametric FDEs with inner function" begin include("3FDE_par_test.jl") end

@safetestset "solving  Simple Harmonic Motion (order=2) " begin include("SimpleHarmonicMotion_test.jl") end

@safetestset "sample test" begin include("sample_test.jl") end


text:
using AutoDiff
#=
x=Variable(2.0)
y=Variable(4.0)

z=(x-1)^2+(y-3)^2
println(z)
minimize(z)
println("x: $(val(x))")
println("y: $(val(y))")

#same computation as above, vectorized
z=sum((Lift([Variable(2.0), Variable(5.0)])-[1.0,3.0]).^2)
println(z)
minimize(z)
println("x: $(val(x))")
println("y: $(val(y))")
=#

using MNIST
using Images,ImageView

N=60000
data=traindata()[1][:,1:N]/256
raw_labels=traindata()[2][1:N]
labels=zeros(Float64,(10,N))
for i in 1:length(raw_labels)
	labels[raw_labels[i]+1,i]=1.0
end

currdigit=0
function getdigit()
	global currdigit=(currdigit+1)%N
	data[:,currdigit+1]
end
function getlabel(digit)
	labels[:,currdigit+1]
end

digit=Call(getdigit,())
label=Call(getlabel,(digit,))

weights1=Variable(0.01*randn((10,784)))
biases1=Variable(0.01*randn((10,)))
function prediction(digit)
	layer1=weights1*digit.+biases1
	tmp=exp(layer1)
	tmp./sum(tmp)
end

begin #viewing
	function f(x)
		a=x-minimum(x)
		a/=maximum(a)
		grayim(transpose(reshape(a,(28,28))))
	end
	c = canvasgrid(4,5)
	function updateview()
		for i in 1:10
			view(c[i],f(val(weights1)[i,:]))
		end
	end
end

error=sum((prediction(digit)-label).^2)
running_mean=0
using ProfileView
minimize(error,delta=0.0001,its=10)
Profile.clear_malloc_data()
minimize(error,delta=0.0001,its=10000000,f=(it->begin
	global running_mean
	running_mean=0.9999*running_mean+0.0001*val(error)
	if it%1000==0
		println("Mean error: $(running_mean)")
		updateview()
	end
end
))
ProfileView.view()


text:
# # Unicycle Car Model
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Sherlock-Benchmark-10-Unicycle.ipynb)

module Unicycle  #jl

using NeuralNetworkAnalysis, MAT
using NeuralNetworkAnalysis: UniformAdditivePostprocessing

# This benchmark is that of a unicycle model of a car [^DCS19] taken from Benchmark
# 10 of the Sherlock tool. It models the dynamics of a car involving 4
# variables, specifically the $x$ and $y$ coordinates on a 2 dimensional
# plane, as well as velocity magnitude (speed) and steering angle.

# ## Model
#
# The model is a four-dimensional system given by the following equations:
#
# ```math
# \left\{ \begin{array}{lcl}
# \dot{x}_1 &=& x_4 \cos(x_3) \\
# \dot{x}_2 &=& x_4 \sin(x_3) \\
# \dot{x}_3 &=& u_2 \\
# \dot{x}_4 &=& u_1 + w
# \end{array} \right.
# ```
# where $w$ is a bounded error in the range $[−1e−4, 1e−4]$. A neural
# network controller was trained for this system. The trained network has 1
# hidden layer with 500 neurons. Note that the output of the neural network
# $f(x)$ needs to be normalized in order to obtain $(u_1, u_2)$, namely
# $u_i = f(x)_i − 20$. The sampling time for this controller is 0.2s.

# We model the error $w$ as a nondeterministically assigned constant.
@taylorize function unicycle!(dx, x, p, t)
    x₁, x₂, x₃, x₄, w, u₁, u₂ = x

    dx[1] = x₄ * cos(x₃)
    dx[2] = x₄ * sin(x₃)
    dx[3] = u₂
    dx[4] = u₁ + w
    dx[5] = zero(x[5])
    dx[6] = zero(x[6])
    dx[7] = zero(x[7])
    return dx
end

controller = read_nnet_mat(@modelpath("Sherlock-Benchmark-10-Unicycle",
                                      "controllerB_nnv.mat");
                           act_key="act_fcns");

# ## Specification
#
# The verification problem here is that of reachability. For an initial set of,
# $x_1 ∈ [9.5,9.55], x_2 ∈ [−4.5,−4.45], x_3 ∈ [2.1,2.11], x_4 ∈ [1.5,1.51]$,
# its is required to prove that the system reaches the set
# $x_1 ∈ [−0.6,0.6], x_2 ∈ [−0.2,0.2], x_3 ∈ [−0.06,0.06], x_4 ∈ [−0.3,0.3]$
# within a time window of 10s.

X₀ = Hyperrectangle(low=[9.5, -4.5, 2.1, 1.5, -1e-4],
                    high=[9.55, -4.45, 2.11, 1.51, 1e-4])
U₀ = ZeroSet(2)
vars_idx = Dict(:state_vars=>1:4, :input_vars=>[5], :control_vars=>6:7)
ivp = @ivp(x' = unicycle!(x), dim: 7, x(0) ∈ X₀ × U₀)

period = 0.2  # control period

control_postprocessing = UniformAdditivePostprocessing(-20.0)

prob = ControlledPlant(ivp, controller, vars_idx, period;
                       postprocessing=control_postprocessing)

## Safety specification: [x[1], x[2], x[3], x[4]] ∈ ±[0.6, 0.2, 0.06, 0.3] for some t ≤ T
T = 10.0
T_warmup = 2 * period  # shorter time horizon for dry run

target_states = cartesian_product(Hyperrectangle(zeros(4), [0.6, 0.2, 0.06, 0.3]),
                                  Universe(3))
predicate = X -> X ⊆ target_states
predicate_sol = sol -> any(predicate(R) for F in sol for R in F);

## sufficient check: only look at the final time point
predicate_R_tend = R -> overapproximate(R, Zonotope, tend(R)) ⊆ target_states
predicate_R_all = R -> R ⊆ target_states
predicate_sol_suff = sol -> predicate_R_all(sol[end]);

# ## Results

alg = TMJets(abstol=1e-15, orderT=10, orderQ=1)
alg_nn = Ai2()
splitter = BoxSplitter([3, 1, 8, 1])

function benchmark(; T=T, silent::Bool=false)
    ## We solve the controlled system:
    silent || println("flowpipe construction")
    res_sol = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg,
                                 splitter=splitter)
    sol = res_sol.value
    silent || print_timed(res_sol)

    ## Next we check the property for an overapproximated flowpipe:
    silent || println("property checking")
    res_pred = @timed predicate_sol_suff(sol)
    silent || print_timed(res_pred)
    if res_pred.value
        silent || println("The property is satisfied.")
    else
        silent || println("The property may be violated.")
    end
    return sol
end;

benchmark(T=T_warmup, silent=true)  # warm-up
res = @timed benchmark()  # benchmark
sol = res.value
println("total analysis time")
print_timed(res);

# We also compute some simulations:

import DifferentialEquations

println("simulation")
res = @timed simulate(prob, T=T; trajectories=10, include_vertices=false)
sim = res.value
print_timed(res);

# Finally we plot the results:

using Plots
import DisplayAs

solz = overapproximate(sol, Zonotope)
Tint = try convert(Int, T) catch; T end

function plot_helper(fig, vars; show_simulation::Bool=true)
    if vars[1] == 0
        target_states_projected = project(target_states, [vars[2]])
        time = Interval(0, T)
        target_states_projected = cartesian_product(time, target_states_projected)
    else
        target_states_projected = project(target_states, vars)
    end
    plot!(fig, solz, vars=vars, color=:yellow, lab="")
    plot!(fig, target_states_projected, color=:cyan, alpha=0.5, lab="target states")
    if show_simulation
        plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
    end
end

vars = (1, 2)
fig = plot(xlab="x₁", ylab="x₂")
plot_helper(fig, vars)
## savefig("Unicycle-x1-x2.png")
fig = DisplayAs.Text(DisplayAs.PNG(fig))

#-

fig = plot(xlab="x₁", ylab="x₂")
plot_helper(fig, vars; show_simulation=false)
xlims!(0, 1)
ylims!(-0.5, 0.5)
plot!(fig, overapproximate(sol[end], Zonotope, tend(solz[end])), vars=vars,
      color=:orange, lab="reach set at t = $Tint")
## savefig("Unicycle-close-x1-x2.png")
fig = DisplayAs.Text(DisplayAs.PNG(fig))

#-

vars=(3, 4)
fig = plot(xlab="x₃", ylab="x₄", leg=:bottomright)
plot_helper(fig, vars)
## savefig("Unicycle-x3-x4.png")
fig = DisplayAs.Text(DisplayAs.PNG(fig))

#-

fig = plot(xlab="x₃", ylab="x₄", leg=:bottomright)
plot_helper(fig, vars; show_simulation=false)
xlims!(-0.1, 0.1)
ylims!(-0.4, 0)
plot!(fig, overapproximate(sol[end], Zonotope, tend(solz[end])), vars=vars,
      color=:orange, lab="reach set at t = $Tint")
## savefig("Unicycle-close-x3-x4.png")
fig = DisplayAs.Text(DisplayAs.PNG(fig))

# ## References

# [^DCS19]: Souradeep Dutta, Xin Chen, and Sriram Sankaranarayanan. *Reachability
#           analysis for neural feedback systems using regressive polynomial rule
#           inference.* In [Proceedings of the 22nd ACMInternational Conference on Hybrid
#           Systems: Computation and Control, HSCC 2019, Montreal,QC, Canada,
#           April 16-18, 2019., pages 157–168, 2019](https://dl.acm.org/doi/abs/10.1145/3302504.3311807).

end  #jl
nothing  #jl


text:
"""
    OneHotArray

A type to hold a one-hot encoded array efficiently.
Implements efficient matrix multiplication.
"""
struct OneHotArray{N,A} <: AbstractArray{Bool,N}
    c::A
    q::Int
    function OneHotArray(c::Union{Int, AbstractArray{Int}}, q::Int = maximum(c))
        if !all(1 ≤ x ≤ q for x in c)
            throw(ArgumentError("c must be ≥ 1 and ≤ $q; got $x"))
        end
        return new{ndims(c) + 1, typeof(c)}(c, q)
    end
end

OneHotArray(X::AbstractArray{Bool}) = OneHotArray(decode(X))
OneHotArray(X::OneHotArray) = X

Base.size(a::OneHotArray) = (a.q, size(a.c)...)

@inline function Base.getindex(a::OneHotArray{N}, i::Vararg{Int,N}) where {N}
    @boundscheck checkbounds(a, i...)
    return @inbounds first(i) == a.c[tail(i)...]
end

@inline function Base.getindex(a::OneHotArray, ::Colon, i...)
    @boundscheck checkbounds(a.c, i...)
    return @inbounds OneHotArray(a.c[i...], a.q)
end

const OneHotVector = OneHotArray{1,Int}
const OneHotMatrix = OneHotArray{2,Vector{Int}}
const OneHotVecOrMat = Union{OneHotVector, OneHotMatrix}

function Base.:*(A::AbstractMatrix, B::OneHotVecOrMat)
    if size(A, 2) ≠ size(B, 1)
        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
    end
    return A[:, B.c]
end

"""
    hotmul(A, B::OneHotArray, Val(dims))

`A*B` contracting first `dims` dimensions of `A` with first `dims` dimensions of `B`
(which must match). Exploits the fact that `B` is a `OneHotArray`.
"""
function hotmul(A::AbstractArray, B::OneHotArray, ::Val{dims}) where {dims}
	dims::Int
	@assert ndims(A) ≥ dims && ndims(B) ≥ dims

	A_contracted = ntuple(d -> size(A, d), Val(dims))
	B_contracted = ntuple(d -> size(B, d), Val(dims))
	@assert A_contracted == B_contracted
	A_non_contracted = ntuple(d -> size(A, dims + d), Val(ndims(A) - dims))
	B_non_contracted = ntuple(d -> size(B, dims + d), Val(ndims(B) - dims))
	@assert size(A) == (A_contracted..., A_non_contracted...)
	@assert size(B) == (B_contracted..., B_non_contracted...)

	Amat = reshape(A, :, prod(A_non_contracted))
	Bmat = reshape(B, :, prod(B_non_contracted))
	Cmat = Amat' * Bmat
	C_size = (A_non_contracted..., B_non_contracted...)
	C = reshape(Cmat, C_size)
	return C
end

"""
    tensormul_ff(A, B, Val(dims))

`A*B` contracting first `dims` dimensions of `A` with first `dims` dimensions of `B`
(which must match).
"""
function tensormul_ff(A::AbstractArray, B::AbstractArray, ::Val{dims}) where {dims}
	dims::Int
	@assert ndims(A) ≥ dims && ndims(B) ≥ dims
	A_contracted = ntuple(d -> size(A, d), Val(dims))
	B_contracted = ntuple(d -> size(B, d), Val(dims))
	A_non_contracted = ntuple(d -> size(A, dims + d), Val(ndims(A) - dims))
	B_non_contracted = ntuple(d -> size(B, dims + d), Val(ndims(B) - dims))
	@assert A_contracted == B_contracted
	@assert size(A) == (A_contracted..., A_non_contracted...)
	@assert size(B) == (B_contracted..., B_non_contracted...)
	Amat = reshape(A, :, prod(A_non_contracted))
	Bmat = reshape(B, :, prod(B_non_contracted))
	Cmat = Amat' * Bmat
	C_size = (A_non_contracted..., B_non_contracted...)
	C = reshape(Cmat, C_size)
	return C
end


text:
#######################################
### helper functions to check input ###
#######################################

#' Check Input for Functions in the bigleaf Package
#'
#' Checks length and type of the provided input variables.
#'
#' - data   Input DataFrame or matrix
#' - ...    Input variables. Either a list or individual vectors
#'
#' # Note
#' This function can be run for named variables (in which case the return
#'       value will be named according to the given name), or for placeholder
#'       variables that are assigned and named according to e.g. entries of a character
#'       vector. In the latter case, the input variable has to be named as `"var"` or
#'       `"var_qc"`.
#'
#' @keywords internal
function check_input(data,...)

  vars = check_length(list(...))

  if (missing(data))
    data = nothing
end

  varlist  = match_call()[-c(1:2)]
  varnames = c(unlist(sapply(varlist,as_character)))
  varnames = varnames[!varnames in c("c","list")]

  for (v in seq_along(vars))

    var     = vars[[v]]
    varname = ifelse(varnames[v] in c("var","var_qc"),gsub("\"","",deparse(substitute(var))),varnames[v])

    if (is_character(var))
      if (!missing(data) & !isnothing(data))
        if (length(var) == 1)
          if (var in colnames(data))
            var = data[,var]
            if (is_numeric(var))
              assign(varname,var,pos=sys_frame(-1))
else
              stop("column representing '",varname,"' in the input matrix/DataFrame must be numeric",call.=false)
end
else
            stop ("there is no column named '",var,"' in the input matrix/DataFrame. Indicate the name of the column representing variable '",varname,"', or alternatively, provide a numeric vector of the same length as the input matrix/DataFrame or of length 1.",call.=false)
end
else
          stop("name of variable '",varname,"' must have length 1",call.=false)
end
else
        if ("data" in names(formals(sys_function(which=-1))))
          if (var in as_character(unlist(match_call(definition=sys_function(-1),call=sys_call(-1))[-1])))
            stop("variable '",var,"' is of type character and interpreted as a column name, but no input matrix/DataFrame is provided. Provide '",var,"' as a numeric vector, or an input matrix/DataFrame with a column named '",var,"'",call.=false)
else
            stop("variable '",var,"' is not provided",call.=false)
end
else
          stop("variable '",var,"' must be numeric",call.=false)
end
end
else
      if (length(var) < 2)
        if (isnothing(var))
          assign(varname,var,pos=sys_frame(-1))
          next
elseif (ismissing(var))
          assign(varname,var,pos=sys_frame(-1))
          next
end
end
      if (!missing(data) & !isnothing(data))
        if (is_numeric(var) & length(var) == nrow(data))
          assign(varname,var,envir=sys_frame(-1))
elseif (is_numeric(var) & length(var) != nrow(data))
          if (length(var) == 1)
            var = rep(var,length=nrow(data))
            assign(varname,var,envir=sys_frame(-1))
else
            stop("variable '",varname,"' must have the same length as the input matrix/DataFrame or length 1. Do NOT provide an input matrix/DataFrame if none of its variables are used!",call.=false)
end
elseif (!is_numeric(var))
          stop("variable '",varname,"' must be numeric",call.=false)
end
else
        if (is_numeric(var))
          assign(varname,var,envir=sys_frame(-1))
else
          stop("variable '",varname,"' must be numeric",call.=false)
end
end
end
end
end

#' Test Variables for Equal Length
#'
#' - varlist List of variables for which the length has to be compared
#'
#' # Note
#' This function only plays a role if no input DataFrame or matrix are
#'       provided. In this case it ensures that provided vectors have the same
#'       length to avoid trouble later up the function call.
#'
#' @keywords internal
function check_length(varlist)

  if (is_list(unlist(varlist,recursive=false)))
    varlist = unlist(varlist,recursive=false)
end

  length_vars = sapply(varlist,length)
  length_vars = length_vars[length_vars > 0]

  if (length(unique(length_vars)) >= 2)
    if (sort(unique(length_vars))[1] != 1 | length(unique(length_vars)) > 2)
      stop("All input variables must have the same length or a length of 1!",call.=false)
end
end
  return(varlist)
end


text:
using MeshCore
using MeshCore: ShapeColl, IncRel, ir_skeleton
using MeshCore: shapedesc, nshapes, ir_code, attribute
import Base.show
using StaticArrays

"""
    Mesh

The type of the mesh.

It stores the incidence relations keyed by the code of the relation.

The incidence relation code is `(d1, d2)`, where `d1` is the manifold dimension
of the shape collection on the left, and `d2` is the manifold dimension of the
shape collection on the right.

The incidence relations are stored with a key consisting of the code and a
string tag. If the string tag is unspecified, it is assumed to be an empty
string.
"""
struct Mesh
    name::String # name of the mesh
    _increls::Dict{Tuple{Tuple{Int64, Int64}, String}, IncRel} # dictionary of incidence relations
end

"""
    Mesh()

Define the mesh with default name and empty dictionary of incidence relations.
"""
function Mesh()
    Mesh("mesh", Dict{Tuple{Tuple{Int64, Int64}, String}, IncRel}())
end

"""
    Mesh(s::String)

Define the mesh named `s` with an empty dictionary of incidence relations.
"""
function Mesh(s::String)
    Mesh(s, Dict{Tuple{Tuple{Int64, Int64}, String}, IncRel}())
end

"""
    load(m::Mesh, filename::String)

Load a mesh (incidence relation) from a MESH file.

!!! note

    No check is performed that the loaded incidence relation is compatible
    with the existing incidence relations in the mesh.
"""
function load(m::Mesh, filename::String)
    conns = import_MESH(filename)
    return attach!(m, conns[1])
end

"""
    save(m::Mesh, filename::String)

Save a mesh base incidence relation to a MESH file.
"""
function save(m::Mesh, filename::String)
    ir = increl(m, basecode(m))
    return export_MESH(filename, ir)
end

"""
    increl(m::Mesh, irc::Tuple{Int64, Int64})

Retrieve the named incidence relation based on the code.

Any tag is matched.
"""
function increl(m::Mesh, irc::Tuple{Int64, Int64})
    for (k, v) in zip(keys(m._increls), values(m._increls))
        if k[1] == irc
            return v
        end
    end
    return nothing
end

"""
    increl(m::Mesh, fullirc::Tuple{Tuple{Int64, Int64}, String})

Retrieve the named incidence relation based on the full key (code + tag).
"""
increl(m::Mesh, fullirc::Tuple{Tuple{Int64, Int64}, String}) = m._increls[fullirc]

"""
    attach!(m::Mesh, increl::IncRel)

Attach the incidence relation under its code and empty tag.

The code of the incidence relation combined with an empty tag (`""`) is the key
under which this relation is stored in the mesh.
"""
function attach!(m::Mesh, ir::IncRel)
    return attach!(m, ir, "")
end

"""
    attach!(m::Mesh, increl::IncRel, tag::String)

Attach the incidence relation under its code and given tag.

The code of the incidence relation combined with the tag is the key
under which this relation is stored in the mesh.
"""
function attach!(m::Mesh, ir::IncRel, tag::String)
    m._increls[(ir_code(ir), tag)] = ir
    return m
end

"""
    basecode(m::Mesh)

Compute the code of the base relation.

The base incidence relation is `(d, 0)` that represents the elements of the
interior of the domain.
"""
function basecode(m::Mesh)
    maxd = 0
    for irc in keys(m._increls)
        maxd = max(maxd, irc[1][1])
    end
    return (maxd, 0)
end

"""
    nspacedims(m::Mesh)

Furnish the dimension of the space in which the mesh lives.
"""
function nspacedims(m::Mesh)
    ir = increl(m, basecode(m))
    a = attribute(ir.right, "geom")
    return length(a[1])
end

"""
    baseincrel(m::Mesh)

Retrieve the base incidence relation for the mesh.
"""
baseincrel(m::Mesh) = increl(m, basecode(m))

"""
    baseincrel(m::Mesh, tag::String)

Retrieve the base incidence relation for the mesh distinguished by its tag.
"""
baseincrel(m::Mesh, tag::String) = increl(m, (basecode(m), tag))

"""
    geometry(m::Mesh)

Retrieve the geometry attribute from the vertices.
"""
function geometry(m::Mesh)
    ir = increl(m, basecode(m))
    return attribute(ir.right, "geom")
end

"""
    Base.summary(m::Mesh)

Form a brief summary of the mesh.
"""
function Base.summary(m::Mesh)
    s = "Mesh $(m.name):"
    for ir in m._increls
        s = s * " $(ir[1]) = " * summary(ir[2]) * "; "
    end
    return s
end

"""
    Base.summary(io::IO, m::Mesh)

Form a brief summary of the mesh.
"""
function Base.summary(io::IO, m::Mesh)
    print(io, summary(m), "\n")
end

"""
    vselect(m::Mesh; kwargs...)

Select vertices. Return as an incidence relation.

Refer to `vselect` that works with the geometry attribute.
"""
function vselect(m::Mesh; kwargs...)
	ir = increl(m, basecode(m))
	geom = attribute(ir.right, "geom")
	list =  vselect(geom; (kwargs...))
	return IncRel(ShapeColl(MeshCore.P1, length(list)), ir.right, [[idx] for idx in list])
end

"""
    boundary(m::Mesh)

Compute the boundary of the mesh.

The incidents relation is stored in the mesh with the tag "boundary".
"""
function boundary(m::Mesh)
	ir = increl(m, basecode(m))
    sir = ir_skeleton(ir, "skeleton") # compute the skeleton of the base incidence relation
    attach!(m, sir) # insert the skeleton into the mesh
    # Now construct the boundary incidence relation
    isboundary = sir.left.attributes["isboundary"]
    ind = [i for i in 1:length(isboundary) if isboundary[i]]
    lft = ShapeColl(shapedesc(sir.left), length(ind), "facets")
    bir = IncRel(lft, sir.right, deepcopy(sir._v[ind]))
    attach!(m, bir, "boundary")
    return bir
end

"""
    vertices(m::Mesh)

Compute the `(0, 0)` incidence relation for the vertices of the base incidence
relation.
"""
function vertices(m::Mesh)
    ir = increl(m, basecode(m))
    return IncRel(ir.right, ir.right, [SVector{1, Int64}([idx]) for idx in 1:nshapes(ir.right)], "vertices")
end

"""
    submesh(m::Mesh, list)

Extract a submesh constructed of a subset of the base relation.
"""
function submesh(m::Mesh, list)
    ir = increl(m, basecode(m))
    lft = ShapeColl(shapedesc(ir.left), length(list))
    v = [ir[idx] for idx in list]
    nir = IncRel(lft, ir.right, v)
    nm = Mesh()
    return attach!(nm, nir)
end

function _label(sc, list, lab)
    if !("label" in keys(sc.attributes))
        sc.attributes["label"] = VecAttrib([zero(typeof(lab)) for idx in 1:nshapes(sc)])
    end
    a = sc.attributes["label"]
    for i in 1:length(list)
        a[list[i]] = lab
    end
end

"""
    label(m::Mesh, irc, list, lab)

Label shapes in `list` with the label `lab`.

Label the shapes on the `shapecoll` of the incidence relation.
`shapecoll` must be either `:left` or `:right`.
"""
function label(m::Mesh, irc, shapecoll, list, lab)
    ir = increl(m, irc)
    if shapecoll == :left
        _label(ir.left, list, lab)
    else
        shapecoll == :right
        _label(ir.right, list, lab)
    end
end


text:
# dataset: https://archive.ics.uci.edu/ml/datasets/Dota2+Games+Results
using CSV
using DataFrames
using Flux: onehotbatch
using Flux.Data: DataLoader
using Parameters: @with_kw
using StatsBase
using Statistics: var
include("../../src/preprocessing/ohe_cat_features.jl")

@with_kw struct TabularDataloaders
    train_loader
    test_loader
    cat_dict
end

function getdata(; batchsize=4096)
    train = CSV.File("input/dota2Train.csv", delim=",", header=0) |> DataFrame
    test = CSV.File("input/dota2Test.csv", delim=",", header=0) |> DataFrame

    X_train_cat = select(train, Not([:Column1, :Column2, :Column3])) |> Array{Float32} |> transpose |> copy
    X_test_cat = select(test, Not([:Column1, :Column2, :Column3])) |> Array{Float32} |> transpose |> copy

    X_train_cont = train[!, [:Column2, :Column3]]|> Array{Float32} |>  transpose
    X_test_cont = test[!, [:Column2, :Column3]] |> Array{Float32} |>  transpose

    # drop cat variables with 1 category
    cat_var = var.(eachrow(X_train_cat))
    X_train_cat = X_train_cat[cat_var.!=0.0, :]
    X_test_cat = X_test_cat[cat_var.!=0.0, :]

    # get categorical dictionary
    cat_dict = get_category_dict(X_train_cat)

    # standardize continious data
    dt = fit(ZScoreTransform, X_train_cont, dims=2)
    X_train_cont = StatsBase.transform(dt, X_train_cont)
    X_test_cont = StatsBase.transform(dt, X_test_cont)

    # One-hot-encode the labels
    y_train = onehotbatch(train[!, :Column1], [-1.0, 1.0])
    y_test = onehotbatch(test[!, :Column1], [-1.0, 1.0])

    # ohe categorical columns
    # results in vector of ohe columns
    X_train_cat = ohe_cat_features(X_train_cat, cat_dict)
    X_test_cat = ohe_cat_features(X_test_cat, cat_dict)

    # Create DataLoaders (mini-batch iterators)
    train_loader = DataLoader((y_train, X_train_cont, X_train_cat...), batchsize=batchsize,  shuffle=true)
    test_loader = DataLoader((y_test, X_test_cont, X_test_cat...), batchsize=batchsize, shuffle=false)

    tabular_dataloaders = TabularDataloaders(train_loader, test_loader, cat_dict)
    return tabular_dataloaders
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Basics (Modia3D/Basics/_module.jl)
#

"""
    getAndCheckFullLibraryPath(dir,libname)

Return joinpath(dir,libname). The returned full path must be a DLL.
It is checked whether this DLL can be opened
"""
function getAndCheckFullLibraryPath(dir,libname)
   full_libname = joinpath(dir,libname)

   # Check that library can be opened (in order to get an early error message)
   dll = Base.Libdl.dlopen_e(full_libname)
   if  dll == C_NULL
      error("SimVis: Unable to load DLL ", full_libname)
   else
      # close dll
      Base.Libdl.dlclose(dll)
   end

   return full_libname
end

"""
    value = getEnvironmentVariable(name, description)

Returns the value of the environment variable `name`.
If `name` does not exist, an error message is triggered:

    Environment variable <name> not defined
    (= <description>)
"""
function getEnvironmentVariable(name::String, description::String)::String
   if haskey(ENV,name)
      return ENV[name]
   else
      error("\n\nEnvironment variable $name not defined\n(= $description)\n")
   end
end


text:
"""
Wrapper for minimization function
"""

## Minimization function
"""
Method to resize bounds if float is passed
"""
function _resize_bounds(lx, nx::Int)
    if length(lx) == 1 && nx > 1
        lx = lx*ones(nx)
    end
    return lx
end

"""
    minimize(func!::Function, x0::Vector, ng::Int; kwargs...)

Minimize function, common interface to IPOPT and SNOPT

Objective function structure:

```julia
function func!(g, x)
    <modify g>
    return <objective value to be minimized>
end
```

# Arguments
    - `func!::Function`: fitness function with signature `rosenbrock!(g, x)` and returning objective
    - `x0::Vector: initial guess
    - `ng::Int: number of constraints
    - `lx::Vector: lower bound on x
    - `ux::Vector: upper bound on x
    - `lg::Vector: lower bound on constraints
    - `ug::Vector: upper bound on constraints
    - `solver::String`: solver, "ipopt" or "snopt"
    - `options::Dict`: options used by solver
    - `sparsity::AbstractSparsityPattern`: sparsity pattern, defualt is DensePattern()
    - `derivatives::AbstractDiffMethod`: derivative method, ForwardFD() or CentralFD() or ForwardAD() or ReverseAD()
    - `outputfile::Boolean`: whether to create output file
    - `verbosity::Int`: verbosirty level, >= 0

# Returns
    - `Tuple`: xstar, fstar, info
"""
function minimize(func!::Function, x0::Vector, ng::Int; kwargs...)
    # unpack values
    lx = _assign_from_kwargs(Dict(kwargs), :lx, -Inf)
    ux = _assign_from_kwargs(Dict(kwargs), :ux, Inf)
    lg = _assign_from_kwargs(Dict(kwargs), :lg, -Inf)
    ug = _assign_from_kwargs(Dict(kwargs), :ug, 0.0)
    solver      = _assign_from_kwargs(Dict(kwargs), :solver, "ipopt")
    options     = _assign_from_kwargs(Dict(kwargs), :options,  Dict())
    sparsity    = _assign_from_kwargs(Dict(kwargs), :sparsity, DensePattern())
    derivatives = _assign_from_kwargs(Dict(kwargs), :derivatives, ForwardFD())
    outputfile  = _assign_from_kwargs(Dict(kwargs), :outputfile, false)
    verbosity   = _assign_from_kwargs(Dict(kwargs), :verbosity, 0)

    # snopt-specific settings
    lencw       = _assign_from_kwargs(Dict(kwargs), :lencw, 500)
    iSumm       = _assign_from_kwargs(Dict(kwargs), :iSumm, 6)   # print to screen if 6

    # initialize number of decision variables
    nx = length(x0)

    # resize bounds if input length is 1
    lx = _resize_bounds(lx, nx)
    ux = _resize_bounds(ux, nx)
    lg = _resize_bounds(lg, ng)
    ug = _resize_bounds(ug, ng)

    # create cache
    # if verbosity>0
    #     println("derivatives: $derivatives")
    # end
    cache = _create_cache(sparsity, derivatives, func!, nx, ng)

    # determine sparsity pattern
    rows, cols = _get_sparsity(sparsity, nx, ng)

    if cmp(solver, "ipopt") == 0
        xstar, fstar, info = minimize_ipopt(options, cache, x0, lx, ux, lg, ug, rows, cols, outputfile)
    elseif cmp(solver, "snopt") == 0
        xstar, fstar, info = minimize_snopt(options, cache, x0, lx, ux, lg, ug, rows, cols, outputfile, lencw, iSumm)
    end
    #print("typeof-cache: ")
    #println(typeof(cache))
    return xstar, fstar, info
end

"""
    _assign_from_kwargs(keyargs::Dict, keyname, default)

Utility to unpack kwargs
"""
function _assign_from_kwargs(keyargs::Dict, keyname, default)
    if haskey(keyargs, keyname)==true
        var = keyargs[keyname]
    else
        var = default
    end
    return var
end

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise