text:
# Some general utility functions. Many of them were part of Flux.jl

"""
    unsqueeze(x; dims)

Return `x` reshaped into an array one dimensionality higher than `x`,
where `dims` indicates in which dimension `x` is extended.

See also [`flatten`](@ref), [`stack`](@ref).

# Examples

```jldoctest
julia> unsqueeze([1 2; 3 4], dims=2)
2×1×2 Array{Int64, 3}:
[:, :, 1] =
 1
 3

[:, :, 2] =
 2
 4

julia> xs = [[1, 2], [3, 4], [5, 6]]
3-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]
 [5, 6]

julia> unsqueeze(xs, dims=1)
1×3 Matrix{Vector{Int64}}:
 [1, 2]  [3, 4]  [5, 6]
```
"""
function unsqueeze(x::AbstractArray; dims::Int)
    sz = ntuple(i -> i < dims ? size(x, i) : i == dims ? 1 : size(x, i - 1), ndims(x) + 1)
    return reshape(x, sz)
end

"""
    unsqueeze(; dims)

Returns a function which, acting on an array, inserts a dimension of size 1 at `dims`.

# Examples

```jldoctest
julia> rand(21, 22, 23) |> unsqueeze(dims=2) |> size
(21, 1, 22, 23)
```
"""
unsqueeze(; dims::Int) = Base.Fix2(_unsqueeze, dims)
_unsqueeze(x, dims) = unsqueeze(x; dims)

Base.show_function(io::IO, u::Base.Fix2{typeof(_unsqueeze)}, ::Bool) = print(io, "unsqueeze(dims=", u.x, ")")

"""
    stack(xs; dims)

Concatenate the given array of arrays `xs` into a single array along the
given dimension `dims`.

See also [`stack`](@ref) and [`batch`](@ref).

# Examples

```jldoctest
julia> xs = [[1, 2], [3, 4], [5, 6]]
3-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]
 [5, 6]

julia> stack(xs, dims=1)
3×2 Matrix{Int64}:
 1  2
 3  4
 5  6

julia> stack(xs, dims=2)
2×3 Matrix{Int64}:
 1  3  5
 2  4  6

julia> stack(xs, dims=3)
2×1×3 Array{Int64, 3}:
[:, :, 1] =
 1
 2

[:, :, 2] =
 3
 4

[:, :, 3] =
 5
 6
```
"""
stack(xs; dims::Int) = cat(unsqueeze.(xs; dims)...; dims)

"""
    unstack(xs; dims)

Unroll the given `xs` into an array of arrays along the given dimension `dims`.

See also [`stack`](@ref) and [`unbatch`](@ref).

# Examples

```jldoctest
julia> unstack([1 3 5 7; 2 4 6 8], dims=2)
4-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]
 [5, 6]
 [7, 8]
```
"""
unstack(xs; dims::Int) = [copy(selectdim(xs, dims, i)) for i in 1:size(xs, dims)]

"""
    chunk(x, n; [dims])

Split `x` into `n` parts. The parts contain the same number of elements
except possibly for the last one that can be smaller.

If `x` is an array, `dims` can be used to specify along which dimension to
split (defaults to the last dimension).

# Examples

```jldoctest
julia> chunk(1:10, 3)
3-element Vector{UnitRange{Int64}}:
 1:4
 5:8
 9:10

julia> x = reshape(collect(1:20), (5, 4))
5×4 Matrix{Int64}:
 1   6  11  16
 2   7  12  17
 3   8  13  18
 4   9  14  19
 5  10  15  20

julia> xs = chunk(x, 2, dims=1)
2-element Vector{SubArray{Int64, 2, Matrix{Int64}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}}:
 [1 6 11 16; 2 7 12 17; 3 8 13 18]
 [4 9 14 19; 5 10 15 20]

julia> xs[1]
3×4 view(::Matrix{Int64}, 1:3, :) with eltype Int64:
 1  6  11  16
 2  7  12  17
 3  8  13  18
```
"""
chunk(x, n::Int) = collect(Iterators.partition(x, ceil(Int, length(x) / n)))

function chunk(x::AbstractArray, n::Int; dims::Int=ndims(x))
    idxs = _partition_idxs(x, n, dims)
    [selectdim(x, dims, i) for i in idxs]
end

function _partition_idxs(x, n, dims)
    bs = ceil(Int, size(x, dims) / n)
    Iterators.partition(axes(x, dims), bs)
end

function rrule(::typeof(chunk), x::AbstractArray, n::Int; dims::Int=ndims(x))
    # this is the implementation of chunk
    idxs = _partition_idxs(x, n, dims)
    y = [selectdim(x, dims, i) for i in idxs]
    valdims = Val(dims)
    chunk_pullback(dy) = (NoTangent(), ∇chunk(unthunk(dy), x, idxs, valdims), NoTangent())

    return y, chunk_pullback
end

# Similar to ∇eachslice  https://github.com/JuliaDiff/ChainRules.jl/blob/8108a77a96af5d4b0c460aac393e44f8943f3c5e/src/rulesets/Base/indexing.jl#L77
function ∇chunk(dys, x::AbstractArray, idxs, vd::Val{dim}) where {dim}
    i1 = findfirst(dy -> !(dy isa AbstractZero), dys)
    if i1 === nothing  # all slices are Zero!
        return _zero_fill!(similar(x, float(eltype(x))))
    end
    T = promote_type(eltype(dys[i1]), eltype(x))
    # The whole point of this gradient is that we can allocate one `dx` array:
    dx = similar(x, T)
    for (k, i) in enumerate(idxs)
        slice = selectdim(dx, dim, i)
        if dys[k] isa AbstractZero
            _zero_fill!(slice)  # Avoids this: copyto!([1,2,3], ZeroTangent()) == [0,2,3]
        else
            copyto!(slice, dys[k])
        end
    end
    return ProjectTo(x)(dx)
end

_zero_fill!(dx::AbstractArray{<:Number}) = fill!(dx, zero(eltype(dx)))
_zero_fill!(dx::AbstractArray) = map!(zero, dx, dx)

function rrule(::typeof(∇chunk), dys, x, idxs, vd::Val{dim}) where dim
    n = length(dys)
    function ∇∇chunk(dz_raw)
        dz = chunk(unthunk(dz_raw), n; dims=dim)
        return (NoTangent(), dz, NoTangent(), NoTangent(), NoTangent())
    end
    return ∇chunk(dys, x, idxs, vd), ∇∇chunk
end

"""
    group_counts(x)

Count the number of times that each element of `x` appears.

See also [`group_indices`](@ref)
# Examples

```jldoctest
julia> group_counts(['a', 'b', 'b'])
Dict{Char, Int64} with 2 entries:
  'a' => 1
  'b' => 2
```
"""
function group_counts(x)
    fs = Dict{eltype(x),Int}()
    for a in x
        fs[a] = get(fs, a, 0) + 1
    end
    return fs
end

"""
    group_indices(x) -> Dict

Computes the indices of elements in the vector `x` for each distinct value contained.
This information is useful for resampling strategies, such as stratified sampling.

See also [`group_counts`](@ref).

# Examples

```jldoctest
julia> x = [:yes, :no, :maybe, :yes];

julia> group_indices(x)
Dict{Symbol, Vector{Int64}} with 3 entries:
  :yes   => [1, 4]
  :maybe => [3]
  :no    => [2]
```
"""
function group_indices(classes::T) where T<:AbstractVector
    dict = Dict{eltype(T), Vector{Int}}()
    for (idx, elem) in enumerate(classes)
        if !haskey(dict, elem)
            push!(dict, elem => [idx])
        else
            push!(dict[elem], idx)
        end
    end
    return dict
end

"""
    batch(xs)

Batch the arrays in `xs` into a single array with
an extra dimension.

If the elements of `xs` are tuples, named tuples, or dicts,
the output will be of the same type.

See also [`unbatch`](@ref).

# Examples

```jldoctest
julia> batch([[1,2,3],
              [4,5,6]])
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

julia> batch([(a=[1,2], b=[3,4])
               (a=[5,6], b=[7,8])])
(a = [1 5; 2 6], b = [3 7; 4 8])
```
"""
function batch(xs)
# Fallback for generric iterables
    @assert length(xs) > 0 "Input should be non-empty"
    data = first(xs) isa AbstractArray ?
        similar(first(xs), size(first(xs))..., length(xs)) :
        Vector{eltype(xs)}(undef, length(xs))
    for (i, x) in enumerate(xs)
        data[batchindex(data, i)...] = x
    end
    return data
end

batchindex(xs, i) = (reverse(Base.tail(reverse(axes(xs))))..., i)

function batch(xs::AbstractArray{<:AbstractArray})
    # Don't use stack(xs, dims=N+1), it is much slower.
    # Here we do reduce(vcat, xs) along with some reshapes.
    szxs = size(xs)
    @assert length(xs) > 0 "Minimum batch size is 1."
    szx = size(xs[1])
    @assert all(x -> size(x) == szx, xs) "All arrays must be of the same size."
    vxs = vec(vec.(xs))
    y = reduce(vcat, vxs)
    return reshape(y, szx..., szxs...)
end

function batch(xs::Vector{<:Tuple})
    @assert length(xs) > 0 "Input should be non-empty"
    n = length(first(xs))
    @assert all(length.(xs) .== n) "Cannot batch tuples with different lengths"
    return ntuple(i -> batch([x[i] for x in xs]), n)
end

function batch(xs::Vector{<:NamedTuple})
    @assert length(xs) > 0 "Input should be non-empty"
    all_keys = [sort(collect(keys(x))) for x in xs]
    ks = all_keys[1]
    @assert all(==(ks), all_keys) "Cannot batch named tuples with different keys"
    NamedTuple(k => batch([x[k] for x in xs]) for k in ks)
end

function batch(xs::Vector{<:Dict})
    @assert length(xs) > 0 "Input should be non-empty"
    all_keys = [sort(collect(keys(x))) for x in xs]
    ks = all_keys[1]
    @assert all(==(ks), all_keys) "cannot batch dicts with different keys"
    Dict(k => batch([x[k] for x in xs]) for k in ks)
end

"""
    unbatch(x)

Reverse of the [`batch`](@ref) operation,
unstacking the last dimension of the array `x`.

See also [`unstack`](@ref).

# Examples

```jldoctest
julia> unbatch([1 3 5 7;
                     2 4 6 8])
4-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]
 [5, 6]
 [7, 8]
```
"""
unbatch(x::AbstractArray) = [getobs(x, i) for i in 1:numobs(x)]
unbatch(x::AbstractVector) = x

"""
    rpad(v::AbstractVector, n::Integer, p)

Return the given sequence padded with `p` up to a maximum length of `n`.

# Examples

```jldoctest
julia> rpad([1, 2], 4, 0)
4-element Vector{Int64}:
 1
 2
 0
 0

julia> rpad([1, 2, 3], 2, 0)
3-element Vector{Int64}:
 1
 2
 3
```
"""
Base.rpad(v::AbstractVector, n::Integer, p) = [v; fill(p, max(n - length(v), 0))]
# TODO Piracy

"""
    batchseq(seqs, pad)

Take a list of `N` sequences, and turn them into a single sequence where each
item is a batch of `N`. Short sequences will be padded by `pad`.

# Examples

```jldoctest
julia> batchseq([[1, 2, 3], [4, 5]], 0)
3-element Vector{Vector{Int64}}:
 [1, 4]
 [2, 5]
 [3, 0]
```
"""
function batchseq(xs, pad = nothing, n = maximum(length(x) for x in xs))
    xs_ = [rpad(x, n, pad) for x in xs]
    [batch([xs_[j][i] for j = 1:length(xs_)]) for i = 1:n]
end

"""
    flatten(x::AbstractArray)

Reshape arbitrarly-shaped input into a matrix-shaped output,
preserving the size of the last dimension.

See also [`unsqueeze`](@ref).

# Examples

```jldoctest
julia> rand(3,4,5) |> flatten |> size
(12, 5)
```
"""
function flatten(x::AbstractArray)
    return reshape(x, :, size(x)[end])
end

"""
    normalise(x; dims=ndims(x), ϵ=1e-5)

Normalise the array `x` to mean 0 and standard deviation 1 across the dimension(s) given by `dims`.
Per default, `dims` is the last dimension.

`ϵ` is a small additive factor added to the denominator for numerical stability.
"""
function normalise(x::AbstractArray; dims=ndims(x), ϵ=ofeltype(x, 1e-5))
    μ = mean(x, dims=dims)
    #   σ = std(x, dims=dims, mean=μ, corrected=false) # use this when Zygote#478 gets merged
    σ = std(x, dims=dims, corrected=false)
    return (x .- μ) ./ (σ .+ ϵ)
end

ofeltype(x, y) = convert(float(eltype(x)), y)
epseltype(x) = eps(float(eltype(x)))

"""
    ones_like(x, [element_type=eltype(x)], [dims=size(x)]))

Create an array with the given element type and size, based upon the given source array `x`.
All element of the new array will be set to 1.
The second and third arguments are both optional, defaulting to the given array's eltype and
size. The dimensions may be specified as an integer or as a tuple argument.

See also [`zeros_like`](@ref) and [`fill_like`](@ref).

# Examples

```julia-repl
julia> x = rand(Float32, 2)
2-element Vector{Float32}:
 0.8621633
 0.5158395

julia> ones_like(x, (3, 3))
3×3 Matrix{Float32}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> using CUDA

julia> x = CUDA.rand(2, 2)
2×2 CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}:
 0.82297   0.656143
 0.701828  0.391335

julia> ones_like(x, Float64)
2×2 CuArray{Float64, 2, CUDA.Mem.DeviceBuffer}:
 1.0  1.0
 1.0  1.0
```
"""
ones_like(x::AbstractArray, T::Type, sz=size(x)) = fill!(similar(x, T, sz), 1)
ones_like(x::AbstractArray, sz=size(x)) = ones_like(x, eltype(x), sz)

"""
    zeros_like(x, [element_type=eltype(x)], [dims=size(x)]))

Create an array with the given element type and size, based upon the given source array `x`.
All element of the new array will be set to 0.
The second and third arguments are both optional, defaulting to the given array's eltype and
size. The dimensions may be specified as an integer or as a tuple argument.

See also [`ones_like`](@ref) and [`fill_like`](@ref).

# Examples

```julia-repl
julia> x = rand(Float32, 2)
2-element Vector{Float32}:
 0.4005432
 0.36934233

julia> zeros_like(x, (3, 3))
3×3 Matrix{Float32}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> using CUDA

julia> x = CUDA.rand(2, 2)
2×2 CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}:
 0.0695155  0.667979
 0.558468   0.59903

julia> zeros_like(x, Float64)
2×2 CuArray{Float64, 2, CUDA.Mem.DeviceBuffer}:
 0.0  0.0
 0.0  0.0
```
"""
zeros_like(x::AbstractArray, T::Type, sz=size(x)) = fill!(similar(x, T, sz), 0)
zeros_like(x::AbstractArray, sz=size(x)) = zeros_like(x, eltype(x), sz)

"""
    rand_like([rng=default_rng()], x, [element_type=eltype(x)], [dims=size(x)])

Create an array with the given element type and size, based upon the given source array `x`.
All element of the new array will be set to a random value.
The last two arguments are both optional, defaulting to the given array's eltype and
size. The dimensions may be specified as an integer or as a tuple argument.

The default random number generator is used, unless a custom one is passed in explicitly
as the first argument.

See also `Base.rand` and [`randn_like`](@ref).

# Examples

```julia-repl
julia> x = ones(Float32, 2)
2-element Vector{Float32}:
 1.0
 1.0

julia> rand_like(x, (3, 3))
3×3 Matrix{Float32}:
 0.780032  0.920552  0.53689
 0.121451  0.741334  0.5449
 0.55348   0.138136  0.556404

julia> using CUDA

julia> CUDA.ones(2, 2)
2×2 CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}:
 1.0  1.0
 1.0  1.0

julia> rand_like(x, Float64)
2×2 CuArray{Float64, 2, CUDA.Mem.DeviceBuffer}:
 0.429274  0.135379
 0.718895  0.0098756
```
"""
rand_like(x::AbstractArray, T::Type, sz=size(x)) = rand!(similar(x, T, sz))
rand_like(x::AbstractArray, sz=size(x)) = rand_like(x, eltype(x), sz)
rand_like(rng::AbstractRNG, x::AbstractArray, T::Type, sz=size(x)) = rand!(rng, similar(x, T, sz))
rand_like(rng::AbstractRNG, x::AbstractArray, sz=size(x)) = rand_like(rng, x, eltype(x), sz)

"""
    randn_like([rng=default_rng()], x, [element_type=eltype(x)], [dims=size(x)])

Create an array with the given element type and size, based upon the given source array `x`.
All element of the new array will be set to a random value drawn from a normal distribution.
The last two arguments are both optional, defaulting to the given array's eltype and
size. The dimensions may be specified as an integer or as a tuple argument.

The default random number generator is used, unless a custom one is passed in explicitly
as the first argument.

See also `Base.randn` and [`rand_like`](@ref).

# Examples
```julia-repl
julia> x = ones(Float32, 2)
2-element Vector{Float32}:
 1.0
 1.0

julia> randn_like(x, (3, 3))
3×3 Matrix{Float32}:
 -0.385331    0.956231   0.0745102
  1.43756    -0.967328   2.06311
  0.0482372   1.78728   -0.902547

julia> using CUDA

julia> CUDA.ones(2, 2)
2×2 CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}:
 1.0  1.0
 1.0  1.0

julia> randn_like(x, Float64)
2×2 CuArray{Float64, 2, CUDA.Mem.DeviceBuffer}:
 -0.578527   0.823445
 -1.01338   -0.612053
```
"""
randn_like(x::AbstractArray, T::Type, sz=size(x)) = randn!(similar(x, T, sz))
randn_like(x::AbstractArray, sz=size(x)) = randn_like(x, eltype(x), sz)
randn_like(rng::AbstractRNG, x::AbstractArray, T::Type, sz=size(x)) = randn!(rng, similar(x, T, sz))
randn_like(rng::AbstractRNG, x::AbstractArray, sz=size(x)) = randn_like(rng, x, eltype(x), sz)

"""
    fill_like(x, val, [element_type=eltype(x)], [dims=size(x)]))

Create an array with the given element type and size, based upon the given source array `x`.
All element of the new array will be set to `val`.
The third and fourth arguments are both optional, defaulting to the given array's eltype and
size. The dimensions may be specified as an integer or as a tuple argument.

See also [`zeros_like`](@ref) and [`ones_like`](@ref).

# Examples

```julia-repl
julia> x = rand(Float32, 2)
2-element Vector{Float32}:
 0.16087806
 0.89916044

julia> fill_like(x, 1.7, (3, 3))
3×3 Matrix{Float32}:
 1.7  1.7  1.7
 1.7  1.7  1.7
 1.7  1.7  1.7

julia> using CUDA

julia> x = CUDA.rand(2, 2)
2×2 CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}:
 0.803167  0.476101
 0.303041  0.317581

julia> fill_like(x, 1.7, Float64)
2×2 CuArray{Float64, 2, CUDA.Mem.DeviceBuffer}:
 1.7  1.7
 1.7  1.7
```
"""
fill_like(x::AbstractArray, val, T::Type, sz=size(x)) = fill!(similar(x, T, sz), val)
fill_like(x::AbstractArray, val, sz=size(x)) = fill_like(x, val, eltype(x), sz)

@non_differentiable zeros_like(::Any...)
@non_differentiable ones_like(::Any...)
@non_differentiable rand_like(::Any...)
@non_differentiable randn_like(::Any...)

function rrule(::typeof(fill_like), x::AbstractArray, val, T::Type, sz)
    function fill_like_pullback(Δ)
        return (NoTangent(), ZeroTangent(), sum(Δ), NoTangent(), NoTangent())
    end
    return fill_like(x, val, T, sz), fill_like_pullback
end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
module PurlinDesigner

#Import local modules.
using ..AISIS10016
using ..AISIS10024
using ..ThinWalledBeam
using ..InternalForces
using ..Mesh
using ..BeamColumn
using ..CrossSection

# export line_strength, free_flange_define
export define_purlin_cross_section

function define_purlin_cross_section(cross_section_type, t, d1, b1, h, b2, d2, α1, α2, α3, α4, α5, r1, r2, r3, r4, n, n_radius)

    if cross_section_type == "Z"

        #Define the top Z purlin surface.   For the top flange, this means the out-to-out dimensions.  For the bottom flange, the interior outside dimensions should be used.

        #First calculate the correction on out-to-out length to go from the outside surface to the inside bottom flange surface.
        delta_lip_bottom = t / tan((π - deg2rad(abs(α2 - α1))) / 2)
        delta_web_bottom = t / tan((π - deg2rad(abs(α3 - α2))) / 2)

        #Note here that the bottom flange and lip dimensions are smaller here.
        ΔL = [d1 - delta_lip_bottom, b1 - delta_lip_bottom - delta_web_bottom, h - delta_web_bottom, b2, d2]
        θ = [α1, α2, α3, α4, α5]

        #Note that the outside radius is used at the top flange, and the inside radius is used for the bottom flange.
        radius = [r1 - t, r2 - t, r3, r4]

        closed_or_open = 1

        purlin = CrossSection.Feature(ΔL, θ, n, radius, n_radius, closed_or_open)

        #Calculate the out-to-out purlin surface coordinates.
        xcoords_out, ycoords_out = CrossSection.get_xy_coordinates(purlin)

        #Calculate centerline purlin coordinates.
        unitnormals = CrossSection.surface_normals(xcoords_out, ycoords_out, closed_or_open)
        nodenormals = CrossSection.avg_node_normals(unitnormals, closed_or_open)
        xcoords_center, ycoords_center = CrossSection.xycoords_along_normal(xcoords_out, ycoords_out, nodenormals, -t/2)

        #Shift y coordinates so that the bottom purlin face is at y = 0.
        ycoords_center = ycoords_center .- minimum(ycoords_center) .+ t/2

        #Shift x coordinates so that the purlin web centerline is at x = 0.
        index = floor(Int, length(xcoords_center)/2)
        xcoords_center = xcoords_center .- xcoords_center[index]

    end

    #Add C section here at some point.

    #Package nodal geometry.
    node_geometry = [xcoords_center ycoords_center]

    #Define cross-section element connectivity and thicknesses.
    num_cross_section_nodes = length(xcoords_center)
    element_info = [1:(num_cross_section_nodes - 1) 2:num_cross_section_nodes ones(num_cross_section_nodes - 1) * t]

    return node_geometry, element_info

end

# struct PurlinLineProperties

#     Fy::Array{Float64, 1}
#     E::Array{Float64, 1}
#     μ::Array{Float64, 1}
#     G::Array{Float64, 1}
#     ho::Array{Float64, 1}
#     h::Array{Float64, 1}
#     t::Array{Float64, 1}
#     b_top::Array{Float64, 1}
#     b_bottom::Array{Float64, 1}
#     d::Array{Float64, 1}
#     θc::Array{Float64, 1}
#     Ixx::Array{Float64, 1}
#     Iyy::Array{Float64, 1}
#     c_y_top::Array{Float64, 1}
#     c_y_bottom::Array{Float64, 1}
#     c_x_left::Array{Float64, 1}
#     c_x_right::Array{Float64, 1}
#     Cw::Array{Float64, 1}
#     Wn::Array{Float64, 1}
#     Iyy_flange::Array{Float64, 1}
#     c_x_flange_left::Array{Float64, 1}
#     c_x_flange_right::Array{Float64, 1}
#     Sxx_top::Array{Float64, 1}
#     Sxx_bottom::Array{Float64, 1}
#     Syy_left::Array{Float64, 1}
#     Syy_right::Array{Float64, 1}
#     Syy_flange_left::Array{Float64, 1}
#     Syy_flange_right::Array{Float64, 1}
#     My_xx_top::Array{Float64, 1}
#     My_xx_bottom::Array{Float64, 1}
#     My_xx::Array{Float64, 1}
#     My_yy_left::Array{Float64, 1}
#     My_yy_right::Array{Float64, 1}
#     My_yy::Array{Float64, 1}
#     My_yy_flange_left::Array{Float64, 1}
#     My_yy_flange_right::Array{Float64, 1}
#     My_yy_flange::Array{Float64, 1}
#     Mcrℓ_xx::Array{Float64, 1}
#     Mcrℓ_yy::Array{Float64, 1}
#     Bcrℓ::Array{Float64, 1}
#     Mcrℓ_yy_flange::Array{Float64, 1}
#     CorZ::Array{Int, 1}
#     f1::Array{Float64, 1}
#     f2::Array{Float64, 1}
#     Lm::Array{Float64, 1}
#     M1::Array{Float64, 1}
#     M2::Array{Float64, 1}
#     kϕ_dist::Array{Float64, 1}
#     Sf_dist::Array{Float64, 1}
#     Mcrd::Array{Float64, 1}
#     a::Array{Float64, 1}
#     kv::Array{Float64, 1}
#     Fcrv::Array{Float64, 1}
#     Mne_xx::Array{Float64, 1}
#     Mne_yy::Array{Float64, 1}
#     Mne_flange_yy::Array{Float64, 1}

# end

# function define(member_definitions, dm, dz, material_properties, section_properties, cross_section_dimensions)

#     #MATERIAL PROPERTIES

#     #Define purlin yield stress.
#     Fy = Mesh.create_line_element_property_array(member_definitions, dm, dz, material_properties, 4, 3)

#     #Define purlin steel elastic modulus.
#     E = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, material_properties, 4, 1)

#     #Define purlin steel Poisson's ratio.
#     μ = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, material_properties, 4, 2)

#     #Define purlin steel shear modulus.
#     G= E./(2 .*(1 .+ μ))

#     #CROSS-SECTION DIMENSIONS

#     #Define out to out purlin web height.
#     ho = Mesh.create_line_element_property_array(member_definitions, dm, dz, cross_section_dimensions, 7, 2)

#     #Define flat purlin web width.
#     h = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 7, 7)

#     #Define purlin base metal thickness.
#     t = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, cross_section_dimensions, 7, 1)

#     #Define purlin out-to-out top flange width.
#     b_top = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, cross_section_dimensions, 3, 3)

#     #Define purlin out-to-out bottom flange width.
#     b_bottom = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, cross_section_dimensions, 3, 3)

#     #Define purlin out-to-out flange lip length.
#     d = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, cross_section_dimensions, 3, 4)

#     #Define purlin lip angle from the horizontal, in degrees.
#     θc = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, cross_section_dimensions, 3, 5)

#     #SECTION PROPERTIES

#     #Define purlin centroidal x-x moment of inertia.
#     Ixx = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)

#     #Define purlin centroidal y-y moment of inertia.
#     Iyy = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 2)

#     #Define y distance from purlin centroid to the top fiber.
#     c_y_top =  Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)  #check this...

#     #Define y distance from purlin centroid to the bottom fiber.
#     c_y_bottom =  Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)  #check this...

#     #Define x distance from purlin centroid to the leftmost fiber.
#     c_x_left =  Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)  #check this...

#     #Define x distance from purlin centroid to the rightmost fiber.
#     c_x_right =  Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)  #check this...

#     #Define purlin warping torsion constant.
#     Cw = Mesh.create_line_element_property_array(member_definitions, dm, dz, sectionProperties, 3, 5)

#     #Define normalized maximum warping stress in the purlin cross-section.
#     Wn = Mesh.create_line_element_property_array(member_definitions, dm, dz, sectionProperties, 3, 6)

#     #Define moment of inertia of the bottom purlin flange.
#     Iyy_flange = Mesh.create_line_element_property_array(member_definitions, dm, dz, FlangeProperties, 3, 2)

#     #Define x distance from purlin bottom flange centroid to the leftmost fiber.
#     c_x_flange_left = Mesh.create_line_element_property_array(member_definitions, dm, dz, FlangeProperties, 3, 2) #check this...

#     #Define x distance from purlin bottom flange centroid to the rightmost fiber.
#     c_x_flange_right = Mesh.create_line_element_property_array(member_definitions, dm, dz, FlangeProperties, 3, 2) #check this...

#     #Calculate the purlin section modulus for the top fiber about the x-x ais.
#     Sxx_top = Ixx ./ c_y_top

#     #Calculate the purlin section modulus for the bottom fiber about the x-x axis.
#     Sxx_bottom = Ixx ./ c_y_bottom

#     #Calculate the purlin section modulus for the leftmost fiber about the y-y ais.
#     Syy_left = Iyy ./ c_y_left

#     #Calculate the purlin section modulus for the rightmost fiber about the y-y axis.
#     Syy_right = Iyy ./ c_y_right

#     #Calculate the bottom purlin flange section modulus for the leftmost fiber about the y-y ais.
#     Syy_flange_left = Iyy_flange ./ c_x_flange_left

#     #Calculate the bottom purlin flange section modulus for the rightmost fiber about the y-y axis.
#     Syy_flange_right = Iyy_flange ./ c_x_flange_right

#     #FIRST YIELD MOMENTS

#     #Define purlin first yield moment about x-x centroidal axis.
#     My_xx_top = Fy .* Sxx_top
#     My_xx_bottom = Fy .* Sxx_bottom
#     My_xx = min.([My_xx_top  My_xx_bottom])

#     #Define purlin first yield moment about y-y centroidal axis.
#     My_yy_left = Fy .* Syy_left
#     My_yy_right = Fy .* Syy_right
#     My_yy = min.([My_yy_left  My_yy_right])

#     #Define purlin bottom flange first yield moment about y-y centroidal axis.
#     My_yy_flange_left = Fy .* Syy_flange_left
#     My_yy_flange_right = Fy .* Syy_flange_right
#     My_yy_flange = min.([My_yy_flange_left  My_yy_flange_right])

#     #CRITICAL ELASTIC LOCAL BUCKLING

#     #Define the purlin critical elastic local buckling moment about the centroidal x-x axis.
#     Mcrℓ_xx = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 7)

#     #Define the purlin critical elastic local buckling moment about the centroidal y-y axis.
#     Mcrℓ_yy = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 8)

#     #Define the purlin critical elastic local buckling moment from torsion.
#     Bcrℓ = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 9)

#     #Define the purlin bottom flange elastic local buckling moment about the centroidal y-y axis.
#     Mcrℓ_yy_flange = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 9)  #check this!!!

#     #CRITICAL ELASTIC DISTORTIONAL BUCKLING

#     #Define the cross-section type, Cee or Zee.
#     CorZ = trunc.(Int, Mesh.create_line_element_property_array(member_definitions, dm, dz, crossSectionDimensions, 3, 6))

#     #Define the stress gradient at the top of the web for bending about the x-x cetroidal axis.
#     f1 = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 9)  #check this!!!

#     #Define the stress gradient at the bottom of the web for bending about the x-x cetroidal axis.
#     f2= Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 9)  #check this!!!

#     #Define the spacing between discrete distortional buckling restraints.
#     Lm = Mesh.create_line_element_property_array(member_definitions, dm, dz, bracingProperties, 6, 3)

#     #Define moment gradient along the purlin line.  Use this to calculate Mcrd.
#     M1=1.0*ones(length(numberOfNodes))  #need to figure this out...
#     M2=1.0*ones(length(numberOfNodes))  #need to figure this out...

#     #Define the rotational stiffness provided by sheathing and connections to the purlin top flange.  Use this to calculate Mcrd.
#     kϕ_dist=0.0*zeros(length(numberOfNodes))  #need to figure this out...

#     #Define the section modulus at the extreme compression fiber.  Use this to calculate Mcrd.
#     Sf_dist = Sxx   #need to figure this out...

#     #Calculate Mcrd along the purlin line.
#     Mcrd=AISIS10016.app23331.(CorZ, t, ho, b, d, θc, E, μ, G, f1, f2, M1, M2, curvatureSign, Lm, kϕ, Sf)

#     #CRITICAL ELASTIC SHEAR BUCKLING STRESS

#     #Define the shear stiffener spacing along the purlin line.
#     a = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, bracingProperties, 6, 4)

#     #Calculate the plate buckling coefficient for pure shear.
#     kv  = AISIS10016.g233.(a, h)

#     #Caclulate the critical elastic shear buckling stress for the purlin web.
#     Fcrv = AISIS10016.g232.(E, μ, kv, h, t)

#     #GLOBAL BUCKLING NOMINAL FLEXURAL STRENGTH

#     #Define purlin global buckling nominal flexural strength about the x-x centroidal axis. Set this equal to the first yield moment because global buckling deformations will be taken into account in the second order analysis.
#     Mne_xx = My_xx

#     #Define purlin global buckling nominal flexural strength about the y-y centroidal axis.
#     Mne_yy = My_yy

#     #Define purlin bottom flange global buckling nominal flexural strength about the y-y centroidal axis.
#     Mne_flange_yy = My_flange_yy

#     purlin_line = PurlinLineProperties(Fy, E, μ, G, ho, h, t, b_top, b_bottom, d, θc, Ixx, Iyy, c_y_top, c_y_bottom, c_x_left, c_x_right, Cw, Wn, Iyy_flange, c_x_flange_left, c_x_flange_right, Sxx_top, Sxx_bottom, Syy_left, Syy_right, Syy_flange_left, Syy_flange_right, My_xx_top, My_xx_bottom, My_xx, My_yy_left, My_yy_right, My_yy, My_yy_flange_left, My_yy_flange_right, My_yy_flange, Mcrℓ_xx, Mcrℓ_yy, Bcrℓ, Mcrℓ_yy_flange, CorZ, f1, f2, Lm, M1, M2, kϕ_dist, Sf_dist, Mcrd, a, kv, Fcrv, Mne_xx, Mne_yy, Mne_flange_yy)

#     return purlin_line

# end

# #Calculate local-global buckling interaction strength along the purlin line.
# function purlin_local_global_buckling_strength(Mne, Mcrℓ, ASDorLRFD)

#     num_nodes = length(Mne)

#     Mnℓ = zeros(Float64, num_nodes)
#     eMnℓ = zeros(Float64, num_nodes)

#     for i in eachindex(Mcrℓxx)
#         Mnℓ[i], eMnℓ[i] =  AISIS10016.f321(Mne[i], Mcrℓ[i], ASDorLRFD)
#     end

#     return Mnℓ, eMnℓ

# end

# #Calculate torsion strength along the purlin line.
# function purlin_torsion_strength(Cw, Fy, Wn, Bcrℓ, ASDorLRFD)

#     num_nodes = length(Cw)

#     Bn = zeros(Float64, numberOfNodes)
#     eBn = zeros(Float64, numberOfNodes)

#     for i in eachindex(Bcrℓ)
#         Bn[i], eBn[i] = AISIS10024.h411(Cw[i], Fy[i], Wn[i], Bcrℓ[i], ASDorLRFD)
#     end

#     return Bn, eBn

# end

# #Calculate distortional buckling strength along the purlin line.
# function purlin_distortional_buckling_strength(My, Mcrd, ASDorLRFD)

#     num_nodes = length(My)

#     Mnd = zeros(Float64, numberOfNodes)
#     eMnd = zeros(Float64, numberOfNodes)

#     for i in eachindex(Mcrd)
#         Mnd[i], eMnd[i] = AISIS10016.f411(My[i], Mcrd[i], ASDorLRFD)
#     end

#     return Mnd, eMnd

# end

# #Calculate shear strength along the purlin line.
# function purlin_shear_strength(E, h, t, Fy, Vcr, ASDorLRFD)

#     Vn = zeros(Float64, numberOfNodes)
#     eVn = zeros(Float64, numberOfNodes)

#     for i in eachindex(Vcr)
#     Vn[i], eVn[i] = AISIS10016.g21(E[i], h[i], t[i], Fy[i], Vcr[i], ASDorLRFD)
#     end

#     return Vn, eVn

# end

# #Calculate the free flange strength along the purlin line.
# function purlin_free_flange_strength(My_yy, Mcrℓ_yyf, ASDorLRFD)

#     Mnℓyy_freeflange = zeros(Float64, numberOfNodes)
#     eMnℓyy_freeflange = zeros(Float64, numberOfNodes)
#     for i in eachindex(Mnℓyy_freeflange)
#         Mnℓyy_freeflange[i], eMnℓyy_freeflange[i] = AISIS10016.f321(My_yyf[i], Mcrℓ_yyf, ASDorLRFD)
#     end

# end

# #Calculated section strengths along a purlin line.
# function section_strengths(ASDorLRFD, dz, dm, member_definitions, section_properties, cross_section_dimensions, material_properties, load_location, bracing_properties, roof_slope, flange_properties)

#     num_nodes = length(dz)+1

#     #calculate strength limit state capacities

#     #strong axis flexure, local-global interaction
#     Fy = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, materialProperties, 4, 3)
#     Ixx = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 1)
#     ho = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 7, 2)
#     ycy = ho./2  #distance from neutral axis to outer fiber
#     Sxx = Ixx./ycy
#     Myxx = Fy.*Sxx
#     Mnexx = Myxx
#     Mcrℓxx = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 7)

#     Mnℓxx = zeros(Float64, numberOfNodes)
#     eMnℓxx = zeros(Float64, numberOfNodes)
#     for i in eachindex(Mcrℓxx)
#         Mnℓxx[i], eMnℓxx[i] =  AISIS10016.f321(Mnexx[i], Mcrℓxx[i], ASDorLRFD)
#     end

#     #weak axis flexure, local-global interaction
#     Iyy = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 2)

#     t = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 7, 1)
#     b = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 3, 3)
#     d = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 3, 4)
#     θc = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 3, 5)

#     #distance from neutral axis to outer fiber
#     ycx = b.+d.*cos.(deg2rad.(θc)) .-t./2
#     Syy = Iyy./ycx
#     Myyy = Fy.*Syy
#     Mneyy = Myyy
#     Mcrℓyy = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 8)

#     Mnℓyy = zeros(Float64, numberOfNodes)
#     eMnℓyy = zeros(Float64, numberOfNodes)
#     for i in eachindex(Mcrℓyy)
#         Mnℓyy[i], eMnℓyy[i] = AISIS10016.f321(Mneyy[i], Mcrℓyy[i], ASDorLRFD)
#     end

#     #torsion
#     Cw = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 5)
#     Fy = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, materialProperties, 4, 3)
#     Wn = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 6)
#     Bcrℓ = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, sectionProperties, 3, 9)

#     Bn = zeros(Float64, numberOfNodes)
#     eBn = zeros(Float64, numberOfNodes)
#     for i in eachindex(Bcrℓ)
#         Bn[i], eBn[i] = AISIS10024.h411(Cw[i], Fy[i], Wn[i], Bcrℓ[i], ASDorLRFD)
#     end

#     #distortional buckling
#     CorZ = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 3, 6)
#     CorZ = trunc.(Int, CorZ)
#     E = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, materialProperties, 4, 1)
#     μ = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, materialProperties, 4, 2)
#     G= E./(2 .*(1 .+ μ))
#     f1=1.0*ones(length(numberOfNodes))
#     f2=-1.0*ones(length(numberOfNodes))
#     Lm = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, bracingProperties, 6, 3)

#     #define moment gradient factor
#     #assume it is 1
#     M1=1.0*ones(length(numberOfNodes))
#     M2=1.0*ones(length(numberOfNodes))
#     curvatureSign = -1.0*ones(length(numberOfNodes))

#     #assume clip stiffness does not help distortional buckling for a standing seam roof
#     #clips are spaced too far apart
#     kϕ=0.0*zeros(length(numberOfNodes))

#     Sf=Sxx

#     Mcrd=AISIS10016.app23331.(CorZ, t, ho, b, d, θc, E, μ, G, f1, f2, M1, M2, curvatureSign, Lm, kϕ, Sf)

#     My=Fy.*Sf

#     Mnd = zeros(Float64, numberOfNodes)
#     eMnd = zeros(Float64, numberOfNodes)
#     for i in eachindex(Mcrd)
#         Mnd[i], eMnd[i] = AISIS10016.f411(My[i], Mcrd[i], ASDorLRFD)
#     end

#     #Web shear strength
#     a = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, bracingProperties, 6, 4)
#     h = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, crossSectionDimensions, 7, 7)
#     kv  = AISIS10016.g233.(a, h)
#     Fcr = AISIS10016.g232.(E, μ, kv, h, t)
#     Vcr = AISIS10016.g231.(h, t, Fcr)

#     Vn = zeros(Float64, numberOfNodes)
#     eVn = zeros(Float64, numberOfNodes)
#     for i in eachindex(Vcr)
#         Vn[i], eVn[i] = AISIS10016.g21(E[i], h[i], t[i], Fy[i], Vcr[i], ASDorLRFD)
#     end

#     #define free flange properties
#     Iyyf = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, FlangeProperties, 3, 2)
#     xcf = Mesh.create_line_element_property_array(memberDefinitions, dm, dz, FlangeProperties, 3, 6)

#     #Free flange strength
#     #Assume yield criterion for now, no local buckling...
#     Sycf_tension = Iyyf ./abs.(xcf)
#     Sycf_comp = Iyyf ./((b .+d .*cos.(deg2rad.(θc)) .- t/2) .-abs.(xcf))
#     Sycf = min.([Sycf_tension Sycf_comp])
#     My_yyf = Fy.*Sycf
#     Mcrℓ_yyf = 999999999999999.0  #no local buckling right now

#     Mnℓyy_freeflange = zeros(Float64, numberOfNodes)
#     eMnℓyy_freeflange = zeros(Float64, numberOfNodes)
#     for i in eachindex(Mnℓyy_freeflange)
#         Mnℓyy_freeflange[i], eMnℓyy_freeflange[i] = AISIS10016.f321(My_yyf[i], Mcrℓ_yyf, ASDorLRFD)
#     end

#     #note nominal capacities are divided by ASD safety factor or
#     #multiplied by LRFD resistance factor here
#     return eMnℓxx, eMnℓyy, eBn, eMnd, eVn, eMnℓyy_freeflange

# end

# #find flexure+torsion D/C
# function bendingTorsionDemandToCapacity(Mxx, Myy, B, Myy_freeflange, eMnℓxx, eMnℓyy, eBn, eMnℓyy_freeflange)

#     #check bending + torsion interaction
#     #ActionM1, ActionM2, ActionB, Interaction
#     interactionCheck = AISIS10024.h42.(Mxx, Myy, B, Myy_freeflange, eMnℓxx, eMnℓyy, eBn, eMnℓyy_freeflange)

#     actionMxx = [x[1] for x in interactionCheck]
#     actionMyy = [x[2] for x in interactionCheck]
#     actionB = [x[3] for x in interactionCheck]
#     actionMyy_freeflange = [x[4] for x in interactionCheck]
#     totalInteraction = [x[5] for x in interactionCheck]

#     demandToCapacity = totalInteraction./1.15

#     return actionMxx, actionMyy, actionB, actionMyy_freeflange, totalInteraction, demandToCapacity

# end

# #find distortional buckling D/C
# function distortionalDemandToCapacity(Mxx,eMnd)

#     #check distortional buckling
#     demandToCapacity=abs.(Mxx./eMnd)

# end

# #find flexure+shear D/C
# function bendingShearDemandToCapacity(Vyy, Mxx, eMnℓxx, eVn)

#     interaction = AISIS10016.h21.(Mxx, Vyy, eMnℓxx, eVn)

#     demandToCapacity = interaction

#     return demandToCapacity

# end

# #find biaxial bending D/C
# function biaxialBendingDemandToCapacity(Mxx, Myy, eMnℓxx, eMnℓyy)

#     #no axial force for now
#     Pbar=zeros(length(Mxx))
#     Pa=ones(length(Mxx))

#     interactionCheck = AISIS10016.h121.(Pbar, Mxx, Myy, Pa, eMnℓxx, eMnℓyy)

#     actionP = [x[1] for x in interactionCheck]
#     actionMxx = [x[2] for x in interactionCheck]
#     actionMyy = [x[3] for x in interactionCheck]
#     totalInteraction = [x[4] for x in interactionCheck]

#     demandToCapacity = totalInteraction

#     return actionP, actionMxx, actionMyy, totalInteraction, demandToCapacity

# end

# #purlin line demand to capacity
# function lineDC(memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad, eMnℓxx, eMnℓyy, eBn, eMnd, eVn, eMnℓyy_freeflange, bridging)

#     #Calculate purlin line deformation and demands.
#     z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad)
#     Mxx = InternalForces.moment(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)
#     Myy = InternalForces.moment(z, beamProperties.dm, -u, beamProperties.E, beamProperties.Iy)
#     Vyy = InternalForces.shear(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)
#     B = InternalForces.bimoment(z, beamProperties.dm, ϕ, beamProperties.E, beamProperties.Cw)

#     #Calculate free flange behavior with shear flow and flexural-torsional buckling.
#     FlangeProperties, Springs, Loads = PurlinDesigner.free_flange_define(memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, bracingProperties, uniformLoad[2], Mxx)
#     uf, vf, ϕf, FreeFlangeProperties = BeamColumn.solve(memberDefinitions, FlangeProperties, materialProperties, Loads, Springs, endBoundaryConditions, bridging)
#     Myyf = InternalForces.moment(z, beamProperties.dm, -uf, FreeFlangeProperties.E, FreeFlangeProperties.Iy)

#     #Get limit state interactions and demand-to-capacities along the purlin line.
#     BTActionMxx, BTActionMyy, BTActionB, BTActionMyy_freeflange, BTTotalInteraction, BTDemandToCapacity = bendingTorsionDemandToCapacity(Mxx, Myy, B, Myyf, eMnℓxx, eMnℓyy, eBn, eMnℓyy_freeflange)
#     distDemandToCapacity = distortionalDemandToCapacity(Mxx,eMnd)
#     MVDemandToCapacity = bendingShearDemandToCapacity(Vyy, Mxx, eMnℓxx, eVn)
#     BBActionP, BBActionM1, BBActionM2, BBTotalInteraction, BBDemandToCapacity = biaxialBendingDemandToCapacity(Mxx, Myy, eMnℓxx, eMnℓyy)

#     #Find worst case demand-to-capacity ratio along the purlin line.
#     purlinDemandToCapacity=maximum([BTDemandToCapacity; distDemandToCapacity; MVDemandToCapacity; BBDemandToCapacity])

#     return purlinDemandToCapacity

# end

# #purlin line load that causes failure
# function rootfinder(q, eps, residual, demandToCapacity, memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad, eMnℓxx, eMnℓyy, eBn, eMnd, eVn, loadAngle, eMnℓyy_freeflange, bridging)

#     #use bisection method of rootfinding to solve for purlin strength

#     maxIterations=10  #hard code this for now

#     for i=1:maxIterations

#         newq=q/demandToCapacity
#         q = q + (newq-q)/2

#         qx = -q*sin(deg2rad(loadAngle))
#         qy = q*cos(deg2rad(loadAngle))
#         uniformLoad=(qx,qy)

#         demandToCapacity=lineDC(memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad, eMnℓxx, eMnℓyy, eBn, eMnd, eVn, eMnℓyy_freeflange, bridging)

#         residual = 1.0 - abs(demandToCapacity)

#         if abs(residual) < eps
#             return q
#         end

#     end

#     error("Maximum rootfinding interation limit exceeded.  A failure load was not found.")

# end

# function line_strength(ASDorLRFD, gravityOrUplift, memberDefinitions, sectionProperties, crossSectionDimensions, materialProperties, loadLocation, bracingProperties, roofSlope, endBoundaryConditions, supports, bridging)

#     dz, z, dm = Mesh.define_line_element(memberDefinitions)
#     numnodes = length(z)

#     #Calculate free flange available strengths.
#     FlangeProperties, Springs, Loads = PurlinDesigner.free_flange_define_line_element(memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, bracingProperties, 0.0, zeros(Float64, numnodes))

#     #Calculate purlin available strengths.
#     eMnℓxx, eMnℓyy, eBn, eMnd,  eVn, eMnℓyy_freeflange  = sectionStrengths(ASDorLRFD, dz, dm, memberDefinitions, sectionProperties, crossSectionDimensions, materialProperties, loadLocation, bracingProperties, roofSlope, FlangeProperties)

#     #Use rootfinding to solve for purlin line available strength.

#     #pick a small load
#     if gravityOrUplift==0
#       q=0.0001
#       loadAngle=roofSlope
#     elseif gravityOrUplift==1
#       q=-0.0001
#       loadAngle=0.0   #for uplift load is always perpendicular to purlin flange
#     end

#     #find DemandToCapacity for this small load
#     qx = -q*sin(deg2rad(loadAngle))
#     qy = q*cos(deg2rad(loadAngle))
#     uniformLoad = (qx,qy)
#     demandToCapacity = lineDC(memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad, eMnℓxx, eMnℓyy, eBn, eMnd, eVn, eMnℓyy_freeflange, bridging)

#     #initialize residual for rootfinding
#     residual=abs(1.0-demandToCapacity)
#     eps=0.01  #residual tolerance, hard coded for now

#     #this spits out the expected failure load
#     eqn = rootfinder(q, eps, residual, demandToCapacity, memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad, eMnℓxx, eMnℓyy, eBn, eMnd, eVn, loadAngle, eMnℓyy_freeflange, bridging)

#     qx = -eqn*sin(deg2rad(loadAngle))
#     qy = eqn*cos(deg2rad(loadAngle))
#     uniformLoad = (qx,qy)

#     #calculate all the final properties and deformations for the purlin line, at failure

#     z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, bracingProperties, endBoundaryConditions, supports, uniformLoad)

#     Mxx = InternalForces.moment(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)
#     Myy = InternalForces.moment(z, beamProperties.dm, -u, beamProperties.E, beamProperties.Iy)
#     Vyy = InternalForces.shear(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)
#     T = InternalForces.torsion(z, beamProperties.dm, ϕ, beamProperties.E, beamProperties.G, beamProperties.J, beamProperties.Cw)
#     B = InternalForces.bimoment(z, beamProperties.dm, ϕ, beamProperties.E, beamProperties.Cw)

#     #Calculate free flange behavior with shear flow and flexural-torsional buckling.
#     FlangeProperties, Springs, Loads = PurlinDesigner.free_flange_define(memberDefinitions, sectionProperties, materialProperties, crossSectionDimensions, bracingProperties, uniformLoad[2], Mxx)
#     uf, vf, ϕf, FreeFlangeProperties = BeamColumn.solve(memberDefinitions, FlangeProperties, materialProperties, Loads, Springs, endBoundaryConditions, bridging)

#     Myyf = InternalForces.moment(z, beamProperties.dm, -uf, FreeFlangeProperties.E, FreeFlangeProperties.Iy)

#     #Get limit state interactions and demand-to-capacities along the purlin line.
#     BTActionMxx, BTActionMyy, BTActionB, BTActionMyy_freeflange, BTTotalInteraction, BTDemandToCapacity = bendingTorsionDemandToCapacity(Mxx, Myy, B, Myyf, eMnℓxx, eMnℓyy, eBn, eMnℓyy_freeflange)
#     distDemandToCapacity = distortionalDemandToCapacity(Mxx,eMnd)
#     MVDemandToCapacity = bendingShearDemandToCapacity(Vyy, Mxx, eMnℓxx, eVn)
#     BBActionP, BBActionMxx, BBActionMyy, BBTotalInteraction, BBDemandToCapacity = biaxialBendingDemandToCapacity(Mxx, Myy, eMnℓxx, eMnℓyy)

#     deformation = NamedTuple{(:u, :v, :ϕ, :uf, :vf, :ϕf)}((u, v, ϕ, uf, vf, ϕf))
#     strengths = NamedTuple{(:eMnℓxx, :eMnℓyy, :eBn, :eMnd, :eVn, :eMnℓyy_freeflange)}((eMnℓxx, eMnℓyy, eBn, eMnd, eVn, eMnℓyy_freeflange))
#     forces = NamedTuple{(:Mxx, :Myy, :Vyy, :T, :B, :Myyf)}((Mxx, Myy, Vyy, T, B, Myyf))
#     interactions = NamedTuple{(:BTMxx, :BTMyy, :BTB, :BTMyy_freeflange, :BTTotal, :BBP, :BBMxx, :BBMyy, :BBTotal)}((BTActionMxx, BTActionMyy, BTActionB, BTActionMyy_freeflange, BTTotalInteraction, BBActionP, BBActionMxx, BBActionMyy, BBTotalInteraction))
#     demand_to_capacity = NamedTuple{(:BT, :dist, :MV, :BB, :envelope)}((BTDemandToCapacity, distDemandToCapacity, MVDemandToCapacity, BBDemandToCapacity, demandToCapacity))

#     return eqn, z, beamProperties, deformation, strengths, forces, interactions, demand_to_capacity

# end

# function free_flange_stiffness(t, E, H, kϕc)

#     Icantilever = 1/12*t^3   #length^4/length for distributed spring

#     #Use Eq. 16 from Gao and Moen (2013) https://ascelibrary.org/doi/abs/10.1061/(ASCE)ST.1943-541X.0000860
#     #kxf = 3*E*Icantilever/H^3
#     kxf = 1/(H^2/kϕc + (H^3/(3*E*Icantilever)))

#     kϕf = E*Icantilever/H

#     return kxf, kϕf

# end

# function free_flange_define(MemberDefinitions, SectionProperties, MaterialProperties, CrossSectionDimensions, BracingProperties, q, Mxx)

#     dz, z, dm = Mesh.define_line_element(MemberDefinitions)

#     numnodes = length(z)

#     CorZ = CrossSectionDimensions[1][6] + 1
#     H = CrossSectionDimensions[1][2]
#     Bc = CrossSectionDimensions[1][3]
#     Dc = CrossSectionDimensions[1][4]
#     θc = CrossSectionDimensions[1][5]
#     t = CrossSectionDimensions[1][1]

#     r = 0.0
#     kipin = 0
#     center = 0
#     n = 5

#     node, elem = CrossSection.CZflange_template(CorZ,H,Bc,Bc,Dc,Dc,r,r,r,r,θc,θc,t,n,n,n,n,n,n,n,n,n,kipin,center)

#     coords = node[:, 2:3]
#     ends = elem[:, 2:4]

#     Af,xcf,ycf,Ixf,Iyf,Ixyf,thetaf,I1f,I2f,Jf,xsf,ysf,Cwf,B1f,B2f,wnf = CrossSection.CUFSMsection_properties(coords, ends)

#     FlangeProperties = [(Af, Ixf, Iyf, Jf, Cwf, xcf, ycf, xsf, ysf)]

#     E = MaterialProperties[1][1]   #consider generalizing this someday

#     kϕc = BracingProperties[1][2]

#     kxf, kϕf = free_flange_stiffness(t, E, H, kϕc)

#     #kx ky kϕ hx hy
#     Springs = [(kxf*ones(numnodes)),(0.0*ones(numnodes)), (kϕf*ones(numnodes)),(0.0*ones(numnodes)),(0.0*ones(numnodes))]

#     #approximate axial force in flange
#     P = -Mxx ./ ((H -t) - 2 * abs(ycf))

#     #There is shear flow in the free flange of a restrained C or Z.
#     #Use Eq. 13b and 15 from Gao and Moen (2013) https://ascelibrary.org/doi/abs/10.1061/(ASCE)ST.1943-541X.0000860

#     Ix = SectionProperties[1][1]
#     c = CrossSectionDimensions[1][8]
#     b = Bc - c

#     if CorZ == 1  #C

#         kH = ((Bc^2*t*H^2)/(4*Ix) + b)/H

#         if xcf > 0   #if flange is oriented left to right from web
#             kH = -kH
#         end

#     elseif CorZ == 2  #Z

#         kH = (H*t*(Bc^2 + 2*Dc*Bc - (2*Dc^2*Bc)/H))/(4*Ix)

#         if xcf > 0   #if flange is oriented left to right from web
#             kH = -kH
#         end

#     end

#     qx = q .* kH

#     #P qx qy ax ay
#     Loads = [P, (qx*ones(numnodes)), (0.0*ones(numnodes)),(0.0*ones(numnodes)),(ycf*ones(numnodes))]

#     return FlangeProperties, Springs, Loads

# end

end #module

ext: jl


text:
## Scatter operations

"""
    scatter_add!(ys, us, xs)

Scatter addition operation.

    ys[:, xs[k]...] = ys[:, xs[k]...] .+ us[:, k...]

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_add! end

"""
    scatter_sub!(ys, us, xs)

Scatter subtraction operation.

    ys[:, xs[k]...] = ys[:, xs[k]...] .- us[:, k...]

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_sub! end

"""
    scatter_mul!(ys, us, xs)

Scatter multiplication operation.

    ys[:, xs[k]...] = ys[:, xs[k]...] .* us[:, k...]

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_mul! end

"""
    scatter_div!(ys, us, xs)

Scatter division operation.

    ys[:, xs[k]...] = ys[:, xs[k]...] ./ us[:, k...]

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_div! end

for op = [:add, :sub, :mul, :div]
    fn = Symbol("scatter_$(op)!")
    @eval function $fn(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
        @simd for k = 1:length(xs)
            k = CartesianIndices(xs)[k]
            ys_v = view(ys, :, xs[k]...)
            us_v = view(us, :, k)
            @inbounds ys_v .= $(name2op[op]).(ys_v, us_v)
        end
        ys
    end
end

"""
    scatter_max!(ys, us, xs)

Scatter maximum operation.

    ys[:, xs[k]...] = max.(ys[:, xs[k]...], us[:, k...])

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_max!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    @simd for k = 1:length(xs)
        k = CartesianIndices(xs)[k]
        ys_v = view(ys, :, xs[k]...)
        us_v = view(us, :, k)
        @inbounds view(ys, :, xs[k]...) .= max.(ys_v, us_v)
    end
    ys
end

"""
    scatter_min!(ys, us, xs)

Scatter minimum operation.

    ys[:, xs[k]...] = min.(ys[:, xs[k]...], us[:, k...])

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_min!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    @simd for k = 1:length(xs)
        k = CartesianIndices(xs)[k]
        ys_v = view(ys, :, xs[k]...)
        us_v = view(us, :, k)
        @inbounds ys_v .= min.(ys_v, us_v)
    end
    ys
end

"""
    scatter_mean!(ys, us, xs)

Scatter mean operation.

    ys[:, xs[k]...] = mean.(ys[:, xs[k]...], us[:, k...])

# Arguments
- `ys`: the destination for `us` to aggregate to. This argument will be mutated.
- `us`: the source data for aggregating.
- `xs`: the mapping for aggregation from source (index) to destination (value).
The index of `xs` is corresponding to the index of `us` and the value of `xs` is
corresponding to the index of `ys`. The value of `xs` can be `Int` or `Tuple` type.

The dimension of `us` must equal to one plus dimension of `xs`. `ys`, `us` and `xs`
must be supported array type and be the same type.`Array`, `StaticArray` and `CuArray`
are currently supported.
"""
function scatter_mean!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    Ns = zero(ys)
    ys_ = zero(ys)
    scatter_add!(Ns, one.(us), xs)
    scatter_add!(ys_, us, xs)
    ys .+= save_div.(ys_, Ns)
    return ys
end

# Support different types of array
for op = ops
    fn = Symbol("scatter_$(op)!")
    @eval function $fn(ys::Array{T}, us::Array{S}, xs::Array{<:IntOrTuple}) where {T<:Real,S<:Real}
        PT = promote_type(T, S)
        $fn(PT.(ys), PT.(us), xs)
    end
    @eval function $fn(ys::StaticArray{<:Tuple,T}, us::StaticArray{<:Tuple,S},
                       xs::StaticArray{<:Tuple,<:IntOrTuple}) where {T<:Real,S<:Real}
        PT = promote_type(T, S)
        $fn(PT.(ys), PT.(us), xs)
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct ZygoteAD

*Experimental feature, not part of stable public API.*

Constructors:

* ```ZygoteAD()```
"""
struct ZygoteAD <: DifferentiationAlgorithm end
export ZygoteAD

function unshaped_gradient!(grad_f_x::AbstractVector{<:Real}, f::Function, x::AbstractVector{<:Real}, diffalg::ZygoteAD)
    value, back = Zygote.pullback(f, x)
    grad_f_x[:] = first(back(Zygote.sensitivity(value)))
    return value
end


text:
export Compose

"""
`Compose(A::AbstractOperator,B::AbstractOperator)`

Shorthand constructor:

`A*B`

Compose different `AbstractOperator`s. Notice that the domain and codomain of the operators `A` and `B` must match, i.e. `size(A,2) == size(B,1)` and `domainType(A) == codomainType(B)`.

```julia
julia> Compose(DFT(16,2),Variation((4,4)))
ℱc*Ʋ  ℝ^(4, 4) -> ℝ^(16, 2)

julia> MatrixOp(randn(20,10))*DCT(10)
▒*ℱc  ℝ^10 -> ℝ^20

```

"""
struct Compose{N, M, L<:NTuple{N,Any}, T<:NTuple{M,Any}} <: AbstractOperator
	A::L
	buf::T       # memory in the bufdle of the operators
end

# Constructors

function Compose(L1::AbstractOperator, L2::AbstractOperator)
	if size(L1,2) != size(L2,1)
		throw(DimensionMismatch("cannot compose operators"))
	end
	if domainType(L1) != codomainType(L2)
		throw(DomainError())
	end
	Compose( L1, L2, Array{domainType(L1)}(undef,size(L2,1)) )
end

Compose(L1::AbstractOperator,L2::AbstractOperator,buf::AbstractArray) =
Compose( (L2,L1), (buf,))

Compose(L1::Compose,       L2::AbstractOperator,buf::AbstractArray) =
Compose( (L2,L1.A...), (buf,L1.buf...))

Compose(L1::AbstractOperator,L2::Compose,       buf::AbstractArray) =
Compose((L2.A...,L1), (L2.buf...,buf))

Compose(L1::Compose,       L2::Compose,       buf::AbstractArray) =
Compose((L2.A...,L1.A...), (L2.buf...,buf,L1.buf...))

#special cases
Scale(coeff,L::Compose) = Compose((L.A[1:end-1]...,Scale(coeff,L.A[end])), L.buf)
Compose(L1::Scale,L2::Eye) = L1
Compose(L1::Eye,L2::Scale) = L2

Compose(L1::AbstractOperator, L2::Eye) = L1
Compose(L1::Eye, L2::AbstractOperator) = L2
Compose(L1::Eye, L2::Eye) = L1

# Mappings

@generated function mul!(y::C, L::Compose{N,M,T1,T2},b::D) where {N,M,T1,T2,C,D}
	ex = :(mul!(L.buf[1],L.A[1],b))
	for i = 2:M
		ex = quote
			$ex
			mul!(L.buf[$i],L.A[$i], L.buf[$i-1])
		end
	end
	ex = quote
		$ex
		mul!(y,L.A[N], L.buf[M])
		return y
	end
end

@generated function mul!(y::D, L::AdjointOperator{Compose{N,M,T1,T2}},b::C) where {N,M,T1,T2,C,D}
	ex = :(mul!(L.A.buf[M],L.A.A[N]',b))
	for i = M:-1:2
		ex = quote
			$ex
			mul!(L.A.buf[$i-1],L.A.A[$i]', L.A.buf[$i])
		end
	end
	ex = quote
		$ex
		mul!(y,L.A.A[1]', L.A.buf[1])
		return y
	end
end

# Properties

size(L::Compose) = ( size(L.A[end],1), size(L.A[1],2) )

fun_name(L::Compose) = length(L.A) == 2 ? fun_name(L.A[2])*"*"*fun_name(L.A[1]) : "Π"

domainType(L::Compose)   = domainType(L.A[1])
codomainType(L::Compose) = codomainType(L.A[end])

is_linear(L::Compose) = all(is_linear.(L.A))
is_diagonal(L::Compose) = is_sliced(L) ? (length(L.A) == 2 && is_diagonal(L.A[2])) : all(is_diagonal.(L.A))
is_invertible(L::Compose) = all(is_invertible.(L.A))

is_sliced(L::Compose) = typeof(L.A[1]) <: GetIndex
is_AAc_diagonal(L::Compose)  = is_sliced(L) && length(L.A) == 2 && is_AAc_diagonal(L.A[2])

diag(L::Compose) = is_sliced(L) ? diag(L.A[2]) : prod(diag.(L.A))
diag_AAc(L::Compose) = is_AAc_diagonal(L) ? diag_AAc(L.A[2]) : error("is_AAc_diagonal( $(typeof(L) ) ) == false")

# utils
function permute(C::Compose, p::AbstractVector{Int})

	i = findfirst( x -> ndoms(x,2) > 1 , C.A)
	P = permute(C.A[i],p)
	AA = (C.A[1:i-1]..., P, C.A[i+1:end]...)
	Compose(AA,C.buf)

end

remove_displacement(C::Compose) = Compose(remove_displacement.(C.A),C.buf)


text:
include("convex_hull.jl")

using FactCheck

# Test 1
# convex_hull function accepts only array of tuples
# If other values was passed, then it raises an exception
facts("Validate convex_hull function arguments") do
  @fact main("STRING") --> "Please provide array of tuples where each element is Int64"
  @fact main(("string","string")) --> "Please provide array of tuples where each element is Int64"
  @fact main([(1.23, 5.55)]) --> "Please provide array of tuples where each element is Int64"
end

# Test 2
# If the number of points is smaller than 3, then return the input
facts("Return input when it contains less than 3 points") do
  @fact convex_hull!([(0,0)]) --> [(0,0)]
  @fact convex_hull!([(0,1),(0,2)]) --> [(0,1),(0,2)]
end

# Test 3
# Detect direction based on the three points
facts("Detect direction based on three points") do
  @fact counter_clockwise((0,1), (1,1),  (2,2))   > 0 --> true
  @fact counter_clockwise((0,1), (-1,-1),(-2,-2)) < 0 --> true
  @fact counter_clockwise((0,1), (0,2),  (0,3))   --> 0
end

# Test 4
facts("Generates array of points and returns it") do
  @fact size(generate_points(10, -30, 30),1) == 10 --> true
  @fact typeof(generate_points(10,-30,30)) == Array{Tuple{Int64,Int64},1} --> true
end

# Test 5
# Check if the convex hull part is build correctly

# general case:
input1 = [(2,3),(3,0),(6,6),(5,2),(5,4),(1,1),(4,5)]
# a case with collinear point:
input2 = [(2,3),(3,0),(6,6),(3,4),(5,2),(5,4),(1,1),(4,5)]

# sorted imputs
sorted_input1 = sort!(union(input1), alg = HeapSort)
sorted_input2 = sort!(union(input2), alg = HeapSort)

# result array
result1 = []
result2 = []

build_convex_hull_part!(result1, sorted_input1, 1:1:length(sorted_input1))
build_convex_hull_part!(result2, sorted_input2, length(sorted_input2):-1:1)

facts("Checking build_convex_hull_part method") do
    @fact result1 --> [(1,1),(2,3),(4,5),(6,6)]
    @fact result2 --> [(6,6),(5,2),(3,0),(1,1)]
end

# Test 6
# Check if convex hull is build correctly
input1 = [(0,3), (1,1), (2,2), (4,4), (0,0), (1,2), (3,1), (3,3)]
input2 = [(3,-3),(-5,-2),(0,-4),(4,3),(-2,0),(2,3),(1,6),(0,-4),(-5,5),(3,3),(-3,0),(6,0),(5,3),(2,2),(6,0),(4,-5),(-2,3),(4,5),(-3,-1),(-3,2)]
input3 = [(0,-4),(-5,-5),(-9,1),(-6,-1),(2,-10),(-4,1),(-1,-4),(2,7),(0,3),(-9,2),(3,-10),(9,5),(8,3),(10,-2),(-6,-7),(0,-9),(-10,-9),(-3,8),(9,4),(2,1),(-8,0),(-6,-9),(-10,3),(6,-7),(2,-2),(4,-5),(0,9),(-2,3),(8,-4),(5,-5)]

result1 = convex_hull!(input1)
result2 = convex_hull!(input2)
result3 = convex_hull!(input3)

facts("Check builing complete convex hull") do
  @fact result1 --> [(0,0),(0,3),(4,4),(3,1)]
  @fact result2 --> [(-5,-2),(-5,5),(1,6),(4,5),(5,3),(6,0),(4,-5),(0,-4)]
  @fact result3 --> [(-10,-9),(-10,3),(-3,8),(0,9),(9,5),(10,-2),(6,-7),(3,-10),(2,-10)]
end

println("All tests passed")

max_stars_repo_name: MrShemek/ConvexHullTDD

max_issues_repo_name: MrShemek/ConvexHullTDD

max_forks_repo_name: MrShemek/ConvexHullTDD


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a PTX-based test compiler, and generate reflection methods for it

function ptx_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false,
                 minthreads=nothing, maxthreads=nothing, blocks_per_sm=nothing,
                 maxregs=nothing, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = PTXCompilerTarget(cap=v"7.0",
                               minthreads=minthreads, maxthreads=maxthreads,
                               blocks_per_sm=blocks_per_sm, maxregs=maxregs)
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function ptx_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function ptx_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function ptx_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function ptx_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    ptx_method = Symbol("ptx_$(method)")
    @eval begin
        $ptx_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $ptx_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function ptx_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end


text:
#!/usr/bin/env julia

# Copyright 2022 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import project1
using Test
using ImageQualityIndexes
using Images, FileIO
using Plots

@testset "project1.jl" begin
    sol = project1.lorenz_solver(0.0, 14)(0.0:20.0:100.0).u
    x = getindex.(sol, 1)
    y = getindex.(sol, 2)
    z = getindex.(sol, 3)
    @test all(isapprox.(x, [0.0, -5.88349, -5.88783, -5.887840, -5.887840, -5.887840], atol=1.0e-3))
    @test all(isapprox.(y, [1.0, -5.87361, -5.88784, -5.88784, -5.887840, -5.887840], atol=1.0e-3))
    @test all(isapprox.(z, [0.0, 13.01934, 12.99998, 12.9999, 13.0000, 13.0000], atol=1.0e-3))
    ENV["GKSwstype"]="nul"
    savefig(project1.lorenz_plot(), "lorenz.png")
    goldfile_path = realpath(dirname(@__FILE__)*"/../images/lorenz_gold.png")
    imgfile_path = realpath(dirname(@__FILE__)*"/../test/lorenz.png")
    img1 = load(goldfile_path)
    img2 = load(imgfile_path)
    @test assess_ssim(img1, img2) >= 0.8
end


text:
## Rendering and plotting functions ##

function render_trace(trace::Trace; kwargs...)
    render_trace!(trace, plot(size=(500, 500)); kwargs...)
end

function render_trace!(trace::Trace, plt=nothing; show_obs=true, kwargs...)
    plt = isnothing(plt) ? plot!() : plt
    points = Array{Float64}(undef, 2, 0)
    for (i, stroke) in enumerate(trace[:points_per_stroke])
        if trace[:strokes => i => :disconnect]
            plot!(points[1, :], points[2,:]; color=:black, xlims=[0, 5], ylims=[0, 5],
                  legend=false, ratio=:equal, kwargs...)
            points = Array{Float64}(undef, 2, 0)
        end
        points = hcat(points, reduce(hcat, stroke))
    end
    plot!(points[1, :], points[2,:]; color=:black, xlims=[0, 5], ylims=[0, 5],
          legend=false, ratio=:equal, kwargs...)
    if show_obs render_obs!(trace[:obs], plt) end
    return plt
end

function render_obs(obs; kwargs...)
    plt = plot(size=(500, 500), ratio=:equal, legend=false, xlims=[0, 5], ylims=[0, 5])
    render_obs!(obs, plt; kwargs...)
end

function render_obs!(obs, plt=nothing; kwargs...)
    plt = isnothing(plt) ? plot!() : plt
    if obs isa AbstractArray && obs[1] isa AbstractVector obs = reduce(hcat, obs) end
    scatter!(obs[1, :], obs[2,:], color=:black, figsize=(500, 500),
             marker=:dot, msize=3, markerstrokewidth=0, kwargs...)
    return plt
end

function render_traces(traces::AbstractVector{<:Trace}, weights::AbstractVector{<:Real},
                       observations=nothing, kwargs...)
    plt = plot(size=(500, 500))
    for (tr, w) in zip(traces, weights)
       if w < 0.001 continue end
       render_trace!(tr, plt; alpha=w, kwargs...)
    end
    if !isnothing(observations)
        plt = render_obs!(observations, plt)
    end
    return plt
end

function render_cb!(t::Int, state, observations;
                    show_plot=true, anim=nothing,
                    plt_buffer=nothing, plt_interval=0.2)
    traces, weights = get_traces(state), get_norm_weights(state)
    plt = render_traces(traces, weights, observations[:,1:t])
    title!(plt, "t = $t")
    if show_plot display(plt) end
    if !isnothing(anim) frame(anim) end
    if !isnothing(plt_buffer)
        T = size(observations)[2]
        plt_interval = plt_interval < 1 ?
            floor(Int, plt_interval * T) : round(Int, plt_interval)
        if mod(t, plt_interval) == 0 || t == T
            push!(plt_buffer, plt)
        end
    end
end


text:
# Load required packages.
using Mimi, ExcelReaders, CSVFiles, DataFrames

# Load helper functions and revenue recycling components being added to the NICE model.
include("helpers.jl")
include("parameters.jl")

include("components/climatedynamics_component.jl")
include("components/co2cycle_component.jl")
include("components/damages_component.jl")
include("components/emissions_component.jl")
include("components/grosseconomy_component.jl")
include("components/neteconomy_component.jl")
include("components/radiativeforcing_component.jl")
include("components/slr_component.jl")
include("components/slrdamages_component.jl")
include("components/welfareNordhaus_component.jl")
include("components/welfareNeutral_component.jl")
#-----------------------#
# ----- Load Data ----- #
#-----------------------#
const model_years = 2005:10:2595

function constructrice(p,objective)

    m = Model()
    set_dimension!(m, :time, model_years)
    set_dimension!(m, :regions, ["US", "EU", "Japan", "Russia", "Eurasia", "China", "India", "MidEast", "Africa", "LatAm", "OHI", "OthAsia"])

    add_comp!(m, grosseconomy, :grosseconomy)
    add_comp!(m, emissions, :emissions)
    add_comp!(m, co2cycle, :co2cycle)
    add_comp!(m, radiativeforcing, :radiativeforcing)
    add_comp!(m, climatedynamics, :climatedynamics)
    add_comp!(m, sealevelrise, :sealevelrise)
    add_comp!(m, sealeveldamages, :sealeveldamages)
    add_comp!(m, damages, :damages)
    add_comp!(m, neteconomy, :neteconomy)
    if objective == "Negishi"
        add_comp!(m, welfareNordhaus, :welfare)
    else
        add_comp!(m, welfareNeutral, :welfare)
    end

    # GROSS ECONOMY COMPONENT
    set_param!(m, :grosseconomy, :al, p[:al])
    set_param!(m, :grosseconomy, :l, p[:l])
    set_param!(m, :grosseconomy, :gama, p[:gama])
    set_param!(m, :grosseconomy, :dk, p[:dk])
    set_param!(m, :grosseconomy, :k0, p[:k0])

    # Note: offset=1 => dependence is on on prior timestep, i.e., not a cycle
    connect_param!(m, :grosseconomy, :I, :neteconomy, :I)

    # EMISSIONS COMPONENT
    set_param!(m, :emissions, :sigma, p[:sigma])
    set_param!(m, :emissions, :MIU, p[:MIU])
    set_param!(m, :emissions, :etree, p[:etree])
    set_param!(m, :emissions, :cost1, p[:cost1])
    set_param!(m, :emissions, :MIU, p[:MIU])
    set_param!(m, :emissions, :expcost2, p[:expcost2])
    set_param!(m, :emissions, :partfract, p[:partfract])
    set_param!(m, :emissions, :pbacktime, p[:pbacktime])

    connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    # CO2 CYCLE COMPONENT
    set_param!(m, :co2cycle, :mat0, p[:mat0])
    set_param!(m, :co2cycle, :mat1, p[:mat1])
    set_param!(m, :co2cycle, :mu0, p[:mu0])
    set_param!(m, :co2cycle, :ml0, p[:ml0])
    set_param!(m, :co2cycle, :b12, p[:b12])
    set_param!(m, :co2cycle, :b23, p[:b23])
    set_param!(m, :co2cycle, :b11, p[:b11])
    set_param!(m, :co2cycle, :b21, p[:b21])
    set_param!(m, :co2cycle, :b22, p[:b22])
    set_param!(m, :co2cycle, :b32, p[:b32])
    set_param!(m, :co2cycle, :b33, p[:b33])

    connect_param!(m, :co2cycle, :E, :emissions, :E)

    # RADIATIVE FORCING COMPONENT
    set_param!(m, :radiativeforcing, :forcoth, p[:forcoth])
    set_param!(m, :radiativeforcing, :fco22x, p[:fco22x])
    set_param!(m, :radiativeforcing, :mat1, p[:mat1])

    connect_param!(m, :radiativeforcing, :MAT, :co2cycle, :MAT)
    connect_param!(m, :radiativeforcing, :MATSUM, :co2cycle, :MATSUM)

    # CLIMATE DYNAMICS COMPONENT
    set_param!(m, :climatedynamics, :fco22x, p[:fco22x])
    set_param!(m, :climatedynamics, :t2xco2, p[:t2xco2])
    set_param!(m, :climatedynamics, :tatm0, p[:tatm0])
    set_param!(m, :climatedynamics, :tatm1, p[:tatm1])
    set_param!(m, :climatedynamics, :tocean0, p[:tocean0])
    set_param!(m, :climatedynamics, :c1, p[:c1])
    set_param!(m, :climatedynamics, :c3, p[:c3])
    set_param!(m, :climatedynamics, :c4, p[:c4])

    connect_param!(m, :climatedynamics, :FORC, :radiativeforcing, :FORC)

    # SEA LEVEL RISE COMPONENT
    set_param!(m, :sealevelrise, :thermeq, p[:thermeq])
    set_param!(m, :sealevelrise, :therm0, p[:therm0])
    set_param!(m, :sealevelrise, :thermadj, p[:thermadj])
    set_param!(m, :sealevelrise, :gsictotal, p[:gsictotal])
    set_param!(m, :sealevelrise, :gsicmelt, p[:gsicmelt])
    set_param!(m, :sealevelrise, :gsicexp, p[:gsicexp])
    set_param!(m, :sealevelrise, :gis0, p[:gis0])
    set_param!(m, :sealevelrise, :gismelt0, p[:gismelt0])
    set_param!(m, :sealevelrise, :gismeltabove, p[:gismeltabove])
    set_param!(m, :sealevelrise, :gismineq, p[:gismineq])
    set_param!(m, :sealevelrise, :gisexp, p[:gisexp])
    set_param!(m, :sealevelrise, :aismelt0, p[:aismelt0])
    set_param!(m, :sealevelrise, :aismeltlow, p[:aismeltlow])
    set_param!(m, :sealevelrise, :aismeltup, p[:aismeltup])
    set_param!(m, :sealevelrise, :aisratio, p[:aisratio])
    set_param!(m, :sealevelrise, :aisinflection, p[:aisinflection])
    set_param!(m, :sealevelrise, :aisintercept, p[:aisintercept])
    set_param!(m, :sealevelrise, :aiswais, p[:aiswais])
    set_param!(m, :sealevelrise, :aisother, p[:aisother])

    connect_param!(m, :sealevelrise, :TATM, :climatedynamics, :TATM)

    set_param!(m, :sealeveldamages, :slrmultiplier, p[:slrmultiplier])
    set_param!(m, :sealeveldamages, :slrelasticity, p[:slrelasticity])
    set_param!(m, :sealeveldamages, :slrdamlinear, p[:slrdamlinear])
    set_param!(m, :sealeveldamages, :slrdamquadratic, p[:slrdamquadratic])

    connect_param!(m, :sealeveldamages, :TOTALSLR, :sealevelrise, :TOTALSLR)
    connect_param!(m, :sealeveldamages, :YGROSS, :grosseconomy, :YGROSS)

    # DAMAGES COMPONENT
    set_param!(m, :damages, :a1, p[:a1])
    set_param!(m, :damages, :a2, p[:a2])
    set_param!(m, :damages, :a3, p[:a3])

    connect_param!(m, :damages, :TATM, :climatedynamics, :TATM)
    # connect_param!(m, :damages, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :damages, :SLRDAMAGES, :sealeveldamages, :SLRDAMAGES)

    # NET ECONOMY COMPONENT
    set_param!(m, :neteconomy, :S, p[:savings])
    set_param!(m, :neteconomy, :l, p[:l])

    connect_param!(m, :neteconomy, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :DAMFRAC, :damages, :DAMFRAC)
    connect_param!(m, :neteconomy, :ABATEFRAC, :emissions, :ABATEFRAC)

    # WELFARE COMPONENT
    connect_param!(m, :welfare, :CPC, :neteconomy, :CPC)
    if objective == "Negishi"
        set_param!(m, :welfare, :l, p[:l])
        set_param!(m, :welfare, :elasmu, p[:elasmu])
        set_param!(m, :welfare, :rr, p[:rr])
        set_param!(m, :welfare, :scale1, p[:scale1])
        set_param!(m, :welfare, :scale2, p[:scale2])
        set_param!(m, :welfare, :alpha, p[:alpha])
    else
        set_param!(m, :welfare, :l, p[:l])
        set_param!(m, :welfare, :eta, 1.5)
        set_param!(m, :welfare, :rho, 0.015)
    end

    return m
end #function that loads parameters into the model and sets connections between components

function get_rice(;objective="neutral",datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))
    params = getrice2010parameters(datafile)
    un_population_data = convert(Matrix, DataFrame(load(joinpath(@__DIR__, "..", "data", "UN_medium_population_scenario.csv"), skiplines_begin=3))[:, 3:end]) ./ 1000
    params[:l] = un_population_data
    m = constructrice(params,objective)

    return m
end #loads parameters and passes it to previous function to return model

max_stars_repo_name: fdennig/RICEupdate

max_issues_repo_name: fdennig/RICEupdate

max_forks_repo_name: fdennig/RICEupdate


text:
# nf data input file format
#=
       1 0 0 1       2 1 0 1       3 1 0 1       4 1 0 1       5 1 0 1
       6 1 0 1       7 1 0 1       8 1 0 1       9 1 0 1      10 1 0 1
      11 1 0 1      12 1 0 1      13 1 0 1      14 1 0 1      15 1 0 1
      ...
=#

function format_nf_line!(l::AbstractString, t::Array{Tuple})
  items=split(l)
  for i in 1:5
    if length(items) < (i-1)*4+1
      break
    end
    tup = (
      parse(Int, items[round.(Int, (i-1)*4+1)]),
      parse(Int, items[round.(Int, (i-1)*4+2)]),
      parse(Int, items[round.(Int, (i-1)*4+3)]),
      parse(Int, items[round.(Int, (i-1)*4+4)])
    )
    push!(t, tup)
  end
  t
end

function read_nf_file(f::AbstractString)
  f = open(f)
  t = Tuple[]
  while !eof(f)
    l = readline(f)
    t = format_nf_line!(l, t)
  end
  close(f)
  t
end

#loads data input file format
#=
    1  0.00000000E+00  0.00000000E+00  0.10416667E-03
    2  0.00000000E+00  0.00000000E+00 -0.41666667E-03
    3  0.00000000E+00  0.00000000E+00  0.20833333E-03
    4  0.00000000E+00  0.00000000E+00 -0.41666667E-03
    ...
=#

function format_loads_line!(l::AbstractString, t::Array{Tuple})
  sl = float(split(l))
  @assert mod(length(sl), 4) == 0
  sl
end

function read_loads_file(f::AbstractString)
  f = open(f)
  t = Tuple[]
  while !eof(f)
    l = readline(f)
    sl = format_loads_line!(l, t)
    for k in 1:length(sl)/4
      ind = round.(Int, (k-1)*4)
      append!(t, [(round.(Int, sl[ind+1]), [sl[ind+2] sl[ind+3] sl[ind+4]])])
    end
  end
  close(f)
  t
end


text:
# You need this for things like Gadfly which are sensitive to the
# default display
using SnoopCompile

snoop_on()


text:
module TestMutate
  using InteractiveUtils # code_native
  include("mutate.jl")
end

module TestPointers
  using InteractiveUtils # code_native
  include("pointers.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ cf8b6a18-2799-11ec-026a-c5d858e1ce1d
begin
	using FileIO
	using ImageIO
	using ImageShow
end

# ╔═╡ ce671147-700a-43a2-b184-69a304db5001
read(`xwininfo`, String) |> Text

# ╔═╡ 2386fe37-2342-4cfd-9077-e6f22495fcce
window_id = "0x9400003"

# ╔═╡ 19ff9e63-c46e-45d5-b1ff-331a0442ae57
cmd = `import -window $window_id png:-`

# ╔═╡ aac40687-e4db-43bd-97fc-7ee1e118fc25
io = read(cmd)

# ╔═╡ 5ab6e1e7-a437-441a-be85-25209a3e5590
img2 = load(IOBuffer(io))

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"

[compat]
FileIO = "~1.11.1"
ImageIO = "~0.5.8"
ImageShow = "~0.3.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "a325370b9dd0e6bf5656a6f1a7ae80755f8ccc46"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.7.2"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "45efb332df2e86f2cb2e992239b6267d97c9e0b6"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.7"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.1"

[[ImageBase]]
deps = ["ImageCore", "Reexport"]
git-tree-sha1 = "c107d7eda71edc2284197e6a2fbb377d38fc4ec1"
uuid = "c817782e-172a-44cc-b673-b171935fbb9e"
version = "0.1.4"

[[ImageCore]]
deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "9a5c62f231e5bba35695a20988fc7cd6de7eeb5a"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.9.3"

[[ImageIO]]
deps = ["FileIO", "Netpbm", "OpenEXR", "PNGFiles", "TiffImages", "UUIDs"]
git-tree-sha1 = "13c826abd23931d909e4c5538643d9691f62a617"
uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
version = "0.5.8"

[[ImageShow]]
deps = ["Base64", "FileIO", "ImageBase", "ImageCore", "OffsetArrays", "StackViews"]
git-tree-sha1 = "e439b5a4e8676da8a29da0b7d2b498f2db6dbce3"
uuid = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
version = "0.3.2"

[[Imath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "87f7662e03a649cffa2e05bf19c303e168732d3e"
uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1"
version = "3.1.2+0"

[[IndirectArrays]]
git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "1.0.0"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MappedArrays]]
git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.1"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[Netpbm]]
deps = ["FileIO", "ImageCore"]
git-tree-sha1 = "18efc06f6ec36a8b801b23f076e3c6ac7c3bf153"
uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
version = "1.0.2"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "c0e9e582987d36d5a61e650e6e543b9e44d9914b"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.7"

[[OpenEXR]]
deps = ["Colors", "FileIO", "OpenEXR_jll"]
git-tree-sha1 = "327f53360fdb54df7ecd01e96ef1983536d1e633"
uuid = "52e1d378-f018-4a11-a4be-720524705ac7"
version = "0.3.2"

[[OpenEXR_jll]]
deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "923319661e9a22712f24596ce81c54fc0366f304"
uuid = "18a262bb-aa17-5467-a713-aee519bc75cb"
version = "3.1.1+0"

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PNGFiles]]
deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
git-tree-sha1 = "e14c485f6beee0c7a8dcf6128bf70b85f1fe201e"
uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
version = "0.3.9"

[[PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "646eed6f6a5d8df6708f15ea7e02a7a2c4fe4800"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.10"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.7.0"

[[StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TiffImages]]
deps = ["ColorTypes", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "OffsetArrays", "OrderedCollections", "PkgVersion", "ProgressMeter"]
git-tree-sha1 = "945b8d87c5e8d5e34e6207ee15edb9d11ae44716"
uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
version = "0.4.3"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═cf8b6a18-2799-11ec-026a-c5d858e1ce1d
# ╠═ce671147-700a-43a2-b184-69a304db5001
# ╠═2386fe37-2342-4cfd-9077-e6f22495fcce
# ╠═19ff9e63-c46e-45d5-b1ff-331a0442ae57
# ╠═aac40687-e4db-43bd-97fc-7ee1e118fc25
# ╠═5ab6e1e7-a437-441a-be85-25209a3e5590
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using StatsBase
using MeshArrays
using DataFrames

"""
    ReadGriddedFields()

Read gridded variables from file using MeshArrays and
return result in uvetc Dictionary.
"""
function ReadGriddedFields()

###### 1) Get gridded variables via MeshArrays.jl

mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write)
nr=8

## Put grid variables in a dictionary:

GridVariables=Dict("XC" => myread(mygrid.path*"XC",MeshArray(mygrid,Float32)),
"YC" => myread(mygrid.path*"YC",MeshArray(mygrid,Float32)),
"XG" => myread(mygrid.path*"XG",MeshArray(mygrid,Float32)),
"YG" => myread(mygrid.path*"YG",MeshArray(mygrid,Float32)),
"dx" => 5000.0)

## Put velocity fields in a dictionary:

t0=0.0 #approximation / simplification
t1=18001.0*3600.0

u0=myread(mygrid.path*"U.0000000001",MeshArray(mygrid,Float32,nr))
u1=myread(mygrid.path*"U.0000018001",MeshArray(mygrid,Float32,nr))
v0=myread(mygrid.path*"V.0000000001",MeshArray(mygrid,Float32,nr))
v1=myread(mygrid.path*"V.0000018001",MeshArray(mygrid,Float32,nr))

kk=3 #3 to match -1406.25 in pkg/flt output
u0=u0[:,kk]; u1=u1[:,kk];
v0=v0[:,kk]; v1=v1[:,kk];

u0=u0./GridVariables["dx"]
u1=u1./GridVariables["dx"]
v0=v0./GridVariables["dx"]
v1=v1./GridVariables["dx"]

## Merge the two dictionaries:

uvetc=Dict("u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1, "t0" => t0, "t1" => t1)

uvetc=merge(uvetc,GridVariables)

## Visualize velocity fields

mskW=myread(mygrid.path*"hFacW",MeshArray(mygrid,Float32,nr))
mskW=1.0 .+ 0.0 * mask(mskW[:,kk],NaN,0.0)
mskS=myread(mygrid.path*"hFacS",MeshArray(mygrid,Float32,nr))
mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)

msk=Dict("mskW" => mskW, "mskS" => mskS)

uvetc=merge(uvetc,msk)

end

"""
    ReadDisplacements(dirIn::String,prec::DataType)

Read displacements from MITgcm output file using MeshArrays
and return as a DataFrame.
"""
function ReadDisplacements(dirIn::String,prec::DataType)

   #load the data into one array
   prec==Float64 ? reclen=8 : reclen=4
   n1=13

   filIn="float_trajectories"
   tmp1=readdir(dirIn)
   tmp1=filter(x -> occursin(filIn,x),tmp1)
   filList=filter(x -> occursin(".data",x),tmp1)
   #hack:
   #filList=filter(x -> occursin("002.002.data",x),tmp1)
   nf=length(filList)

   n2=Array{Int,1}(undef,nf)
   for ff=1:nf
      fil=dirIn*filList[ff]
      #println(fil)
      tmp=stat(fil)
      n2[ff]=Int64(tmp.size/n1/reclen)-1
   end

   arr = Array{prec,2}(undef,(n1+1,sum(n2)));
   ii=0;
   #@softscope for ff=1:nf
   for ff=1:nf
      fil=dirIn*filList[ff]
      fid = open(fil)
      tmp = Array{prec,2}(undef,(n1,n2[ff]+1))
      read!(fid,tmp)
      arr[1:n1,ii+1:ii+n2[ff]] = hton.(tmp[:,2:n2[ff]+1])
      arr[n1+1,ii+1:ii+n2[ff]] .= ff
      ii=ii+n2[ff]
   end

   #sort the whole dataset by time
   jj = sort!([1:ii;], by=i->arr[2,i]); arr=arr[:,jj];
   #arr = sort!(arr, dims=2, by=i->arr[2,i]);

   #nfloats=Int(maximum(arr[1,:]))
   #npoints=counts(Int.(arr[1,:]))

   #reformat data as a DataFrame
   df=DataFrame()
   df.ID=Int.(arr[1,:])
   df.time=Int.(arr[2,:])
   df.lon=arr[3,:]
   df.lat=arr[4,:]
   df.dep=arr[5,:]
   if true
      df.i=arr[6,:]
      df.j=arr[7,:]
      df.k=arr[8,:]
      df.etaN=arr[9,:]
      df.uVel=arr[10,:]
      df.vVel=arr[11,:]
      df.theta=arr[12,:]
      df.salt=arr[13,:]
      df.tile=Int.(arr[14,:])
   end

   nfloats=maximum(df.ID);
   nsteps=maximum(counts(df.ID));

   println("# floats=$nfloats")
   println("# steps=$nsteps")

   return df
end

"""
    myread()

Read a gridded variable from 2x2 tile files. This is used
in `ReadGriddedFields()` with `flt_example/`
"""
function myread(filRoot::String,x::MeshArray)
   prec=eltype(x)
   prec==Float64 ? reclen=8 : reclen=4;

   (n1,n2)=Int64.(x.grid.ioSize ./ 2);
   fil=filRoot*".001.001.data"
   tmp1=stat(fil);
   n3=Int64(tmp1.size/n1/n2/reclen);

   v00=x.grid.write(x)
   for ii=1:2; for jj=1:2;
      fid = open(filRoot*".00$ii.00$jj.data")
      fld = Array{prec,1}(undef,(n1*n2*n3))
      read!(fid,fld)
      fld = hton.(fld)

      n3>1 ? s=(n1,n2,n3) : s=(n1,n2)
      v00[1+(ii-1)*n1:ii*n1,1+(jj-1)*n2:jj*n2,:]=reshape(fld,s)
   end; end;

   return x.grid.read(v00,x)
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Cbc
using JuMP
using MIPLearn
using PyCall
using Test

if is_gurobi_available
    using Gurobi
end

miplearn_tests = pyimport("miplearn.solvers.tests")
traceback = pyimport("traceback")

function _test_solver(optimizer_factory)
    MIPLearn.@python_call miplearn_tests.run_internal_solver_tests(
        JuMPSolver(optimizer_factory),
    )
end

@testset "JuMPSolver" begin
    @testset "Cbc" begin
        _test_solver(Cbc.Optimizer)
    end
    if is_gurobi_available
        @testset "Gurobi" begin
            _test_solver(Gurobi.Optimizer)
        end
    end
end