text:
@testset "Rate type" begin

    default = Yields.Rate{Float64,typeof(Yields.DEFAULT_COMPOUNDING)}

    y = Yields.Constant(0.05)
    @test Yields.__ratetype(y) == Yields.Rate{Float64, Periodic}
    @test Yields.__ratetype(typeof(Periodic(0.05,2))) == Yields.Rate{Float64, Periodic}
    @test Yields.CompoundingFrequency(y) == Periodic(1)

    y = Yields.Constant(Continuous(0.05))
    @test Yields.__ratetype(y) == Yields.Rate{Float64, Continuous}
    @test Yields.__ratetype(typeof(Continuous(0.05))) == Yields.Rate{Float64, Continuous}
    @test Yields.CompoundingFrequency(y) == Continuous()

    y = Yields.Step([0.02,0.05], [1,2])
    @test Yields.__ratetype(y) == Yields.Rate{Float64, Periodic}
    @test Yields.CompoundingFrequency(y) == Periodic(1)

    y = Yields.Forward( [0.01,0.02,0.03] )
    @test Yields.__ratetype(y) == default

    rates =[0.01, 0.01, 0.03, 0.05, 0.07, 0.16, 0.35, 0.92, 1.40, 1.74, 2.31, 2.41] ./ 100
    mats = [1/12, 2/12, 3/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30]

    y = Yields.CMT(rates,mats)
    @test Yields.__ratetype(y) == default
    @test Yields.CompoundingFrequency(y) == Continuous()

    combination = y + y
    @test Yields.__ratetype(combination) == default

end

@testset "type coercion" begin
    @test Yields.__coerce_rate(0.05, Periodic(1)) == Periodic(0.05,1)
    @test Yields.__coerce_rate(Periodic(0.05,12), Periodic(1)) == Periodic(0.05,12)
end


text:
using LinearAlgebraicRepresentation, ViewerGL
Lar = LinearAlgebraicRepresentation
GL = ViewerGL

filename = "./test/svg/Lar.svg"
V,EV = Lar.svg2lar(filename)

GL.VIEW([
	GL.GLLines(V,EV,GL.COLORS[1]),
	GL.GLFrame2
]);

function pointsinout(V,EV, n=2_000_000)
	point_in = []
	point_out = []
	point_on = []
	classify = Lar.pointInPolygonClassification(V,EV)
	for k=1:n
		queryPoint = rand(1,2)
		inOut = classify(queryPoint)
		# println("k = $k, queryPoint = $queryPoint, inOut = $inOut")
		if inOut=="p_in"
			push!(point_in, queryPoint);
		elseif inOut=="p_out"
			push!(point_out, queryPoint);
		elseif inOut=="p_on"
			push!(point_on, queryPoint);
		end
	end
	#GL.GLPoints(queryPoint,GL.COLORS[2])
	return point_in,point_out,point_on
end

using Random
const rnglist = [MersenneTwister() for i in 1:Threads.nthreads()]

function pointsinout_multithreading(V,EV, n=2_000_000,t=Threads.nthreads())

	point_in = [[] for i in 1:t]
	point_out = [[] for i in 1:t]
	point_on = [[] for i in 1:t]

	classify = Lar.pointInPolygonClassification(V,EV)

	Threads.@threads for k=1:n
		queryPoint = rand(rnglist[Threads.threadid()],1,2)
		#queryPoint = rand(1,2)
		inOut = classify(queryPoint)
		if inOut=="p_in"
			push!(point_in[Threads.threadid()], queryPoint);
		elseif inOut=="p_out"
			push!(point_out[Threads.threadid()], queryPoint);
		elseif inOut=="p_on"
			push!(point_on[Threads.threadid()], queryPoint);
		end
	end
	return vcat(point_in...),vcat(point_out...),vcat(point_on...)
end

using Distributed
#nprocs()
#addprocs(1)
#n_cpu = length(Sys.cpu_info()) or Sys.CPU_THREADS

function pointsinout_multiproc(V,EV, n=2_000_000,ws=Distributed.nworkers())

	point_in = [[] for i in 1:ws]
	point_out = [[] for i in 1:ws]
	point_on = [[] for i in 1:ws]
	# al posto di questi canale remoto??

	classify = Lar.pointInPolygonClassification(V,EV)

	@distributed for k=1:n

		# va fatta la chiamata remota
		queryPoint = rand(rnglist[Threads.threadid()],1,2)
		#queryPoint = rand(1,2)
		inOut = classify(queryPoint)
		if inOut=="p_in"
			push!(point_in[Threads.threadid()], queryPoint);
		elseif inOut=="p_out"
			push!(point_out[Threads.threadid()], queryPoint);
		elseif inOut=="p_on"
			push!(point_on[Threads.threadid()], queryPoint);
		end
	end

	#a,b,c, = fetch()

	return vcat(a),vcat(b),vcat(c)
end

#@time points_in, points_out, points_on = pointsinout_multiproc(V,EV)

@time points_in, points_out, points_on = pointsinout_multithreading(V,EV);
@time points_in, points_out, points_on = pointsinout(V,EV);

pointsin = [vcat(points_in...) zeros(length(points_in),1)]
pointsout = [vcat(points_out...) zeros(length(points_out),1)]

polygon = [GL.GLLines(V,EV,GL.COLORS[1])];
in_mesh = [GL.GLPoints(pointsin, GL.COLORS[2])]
out_mesh = [GL.GLPoints(pointsout, GL.COLORS[3])]

result = cat([polygon,in_mesh,out_mesh])
GL.VIEW(result);

#=

n = 2_000_000

with @btime macro:

serial -> 158.054 s (1772000042 allocations: 158.60 GiB)
multithread -> 51.618 s (1772000187 allocations: 158.62 GiB)
multithread with separate RNG for each thread -> 51.141 s (1772000205 allocations: 158.62 GiB)

=#

ext: jl


text:
module JuliaNature

# Write your package code here.

using Statistics

struct Bat_model
    pop_size::Int
    dims::Int
    obj_func
    position
    velocity
    freq_pulse
    pulse_rate
    loudness
    freq_max
    freq_min
    loud_max
    loud_min
    type # 0 - min , 1 - max
end

function initialize(size::Int, dims::Int, obj_func, velocity_range, posit_range, freq_range, loudness_range, type)
    bat_velocity = [rand(velocity_range[1]:.001:velocity_range[2],dims) for i in 1:size]
    bat_position = [rand(posit_range[1]:.001:posit_range[2],dims) for i in 1:size]
    frequencies = [rand(freq_range[1]:.001:freq_range[2],dims) for i in 1:size]
    pulse_rates = [rand(0:.001:1,dims) for i in 1:size]
    loudness = [rand(loudness_range[1]:.001:loudness_range[2],dims) for i in 1:size]
    return Bat_model(size, dims, obj_func, bat_position, bat_velocity,  frequencies, pulse_rates, loudness, [freq_range[1]
    for i in 1:dims], [freq_range[2] for i in 1:dims] , [loudness_range[1] for i in 1:dims], [loudness_range[2] for i in 1:dims], type)
end

function optimize(algo::Bat_model, iter=50, tol= 0.3)
    t = 1
    best_ind = 1
    println(best_ind)
    if algo.type == 0
            best = typemax(typeof(1.0))
        else
            best = typemin(typeof(1.0))
        end
        plus_0 = algo.pulse_rate
        for i in 1: algo.pop_size
            if algo.type == 0
                if algo.obj_func(algo.position[i]...) < best
                    best = algo.obj_func(algo.position[i]...)
                    best_ind = i
                end
            else
                if algo.obj_func(algo.position[i]...) > best
                    best = algo.obj_func(algo.position[i]...)
                    best_ind = i
                end
            end
        end
    println(algo.position[best_ind])
    println(best_ind)
    while t < iter
        alpha = 0.9
        yamma = 0.9
        for i in 1: algo.pop_size
            beta = rand(0:.001:1, algo.dims)
            algo.freq_pulse[i] = algo.freq_min + (algo.freq_max - algo.freq_min) .*  beta
            algo.velocity[i] = algo.velocity[i] + (algo.position[i] - algo.position[best_ind]) .* algo.freq_pulse[i]
            algo.position[i] = algo.position[i] + algo.velocity[i]
            r = rand(0:.01:1, algo.dims)
            if r > algo.pulse_rate[i]
                algo.position[i] = algo.position[best_ind] + 0.01 .* r
            end
            if beta < algo.loudness[i]
                continue
            end
            if algo.type == 0
                if algo.obj_func(algo.position[i]...) < best
                    best = algo.obj_func(algo.position[i]...)
                    best_ind = i
                end
            else
                if algo.obj_func(algo.position[i]...) > best
                    best = algo.obj_func(algo.position[i]...)
                    best_ind = i
                end
            end

            #scaling = 0.01
            #avg = mean(algo.position)
            #e = rand(-1:0.01:1, 1)
            #algo.position[i] = algo.position[i] + scaling .* e .* avg
            #algo.loudness[i] = alpha * algo.loudness[i]
            #algo.pulse_rate[i] = plus_0[i] * (1 - exp(-yamma * t))
        end
        t += 1
        if abs(sum(algo.position[best_ind] - [1, 3])) < tol
            break
        end
    end
    return algo.position[best_ind]
end

end # module

max_stars_repo_name: JuliaNature/JuliaNature

max_issues_repo_name: JuliaNature/JuliaNature

max_forks_repo_name: JuliaNature/JuliaNature


text:
using Revise
using AdFem
using PyPlot

m = 30
n = 30
h = 1/n

xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = [@. (2(x^3/3-x^2/2)); @. (2(y^3/3 - y^2/2))]

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @. x*(1-x)*y*(1-y)

ρ = dt = 1
function step2(u0)
    B = compute_interaction_matrix(m, n, h)
    bc = bcedge("all", m, n, h)
    A, _ = compute_fvm_tpfa_matrix(ones(m*n), bc, zeros(size(bc,1)), m, n, h)
    rhs = ρ / dt * B * u0
    sol = A\rhs
    return sol
end
bc = bcedge("all", m, n, h)

# There are three ways to compute
# int div u dx
# on each element

# rhs = compute_fvm_mechanics_term(u0, m, n, h)
B = compute_interaction_matrix(m, n, h)
rhs = B * u0
# rhs = h^2 * @. -2x*(1-x)-2y*(1-y)

A, _ = compute_fvm_tpfa_matrix(ones(m*n), bc, zeros(size(bc,1)), m, n, h)
sol = A\rhs

sol = step2(u0)

figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fvm_points(sol, m, n, h)
subplot(122)
visualize_scalar_on_fvm_points(p0, m, n, h)


text:
#
# `generate`
#
# Creates a new variable in the DataFrame. Operates row-by-row, therefore operators on scalars are
# expected. Rows for which the `filter` condition is not met are assigned `missing`.
#
# Differences to Stata:
#   - The current observation is denoted by `i` instead of `_n`, e.g. `:varname[i]` instead of `:varname[_n]`
#   - Explicit types are currently not supported. The new variable takes a type that is the result of the
#       evaluated expression in the first row.
#
# Like in Stata, refering to a variable by `:varname` implicitly means the current row: `:varname[i]`.
#
#
# this is the general form of the command
macro generate(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            generate <expression> [if] <expression>
        or
            bysort <varlist> (<varlist>): generate <expression> [if] <expression>
    """)
end

# short form. This is a copy of the above for `gen`.
macro gen(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@generate($t, $by, $sort, $arguments, $filter, $use, $options)
        end
    )
end

macro generate(t::Symbol,
    by::Vector{Symbol},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@generate_byrow!($t, $by, $sort, $arguments, $filter, $options )
        end
    )
end

# this is the specific vesion that leads to generate's that are without `by`/`bysort`
macro generate(t::Symbol,
    by::Nothing,
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Nothing)
    return esc(
        quote
            Douglass.@generate!($t, $sort, $arguments, $filter)
        end
    )
end

macro generate_byrow!(t::Symbol, varlist_by::Vector{Symbol}, varlist_sort::Union{Vector{Symbol}, Nothing}, arguments::Expr, filter::Union{Expr, Nothing}, options::Union{Dict{String,Any}, Nothing})
    # assert that `arguments` is an assignment
    (arguments.head == :(=)) || error("`generate` expects an assignment operation, e.g. :x = :y + :z")
    # replace :varname by :varname[i] if not referenced
    Douglass.ref_quotenodes!(arguments)
    !isnothing(filter) && Douglass.ref_quotenodes!(filter)
    # get the assigned var symbol (note that it's in a QuoteNode)
    assigned_var = arguments.args[1].args[1].value
    # and the QuoteNode
    assigned_var_qn = arguments.args[1].args[1]

    # if the RHS of the assignment expression is currently a symbol, make it an Expr
    #transformation = (typeof(arguments.args[2]) == Symbol) ? Expr(arguments.args[2]) : arguments.args[2]
    transformation = arguments.args[2]
    # return `nothing`s in case the index is < 1
    isexpr(transformation) && replace_invalid_indices!(transformation)

    return esc(
        quote
            # check variable is not present
            ($(assigned_var_qn) ∉ propertynames($t)) || error("Variable $($(assigned_var_qn)) already present in DataFrame.")

            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, vcat($varlist_by, $varlist_sort))
            end
            #determine type of resulting column from the type of the first element
            assigned_var_type = eltype([Douglass.@with($t,$(transformation)) for _n=1])
            $t[!,$(assigned_var_qn)] = missings(assigned_var_type,size($t,1))

            # this is the function that maps every sub-df into its transformed df
            my_f = _df -> Douglass.@with _df begin
                # define _N
                _N = size(_df, 1)
                # fill the new variable, row by row
                for _n in 1:size(_df,1)
                    if (isnothing($filter) ? true : $filter)  # if condition is not satisfied, leave with missing
                        $(assigned_var_qn)[_n] = $(transformation)
                    end
                end
                _df
            end
            # execute, and copy it to another dataframe,
            # otherwise we get copies of the group variables in there as well
            t2 = combine(my_f, groupby($t,$varlist_by, sort=false, skipmissing = false ))

            $t[!,$(assigned_var_qn)] = missings(eltype(t2[!,$(assigned_var_qn)]),size($t,1))
            $t[!,$(assigned_var_qn)] = t2[!,$(assigned_var_qn)]
            $t
        end
    )
end

# version without `by`
macro generate!(t::Symbol,
    varlist_sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing})

    # assert that `arguments` is an assignment
    (arguments.head == :(=)) || error("`generate` expects an assignment operation, e.g. :x = :y + :z")
    # replace :varname by :varname[i] if not referenced
    Douglass.ref_quotenodes!(arguments)
    !isnothing(filter) && Douglass.ref_quotenodes!(filter)
    # get the assigned var symbol and QuoteNode
    if isa(arguments.args[1], Symbol)
        assigned_var_qn = QuoteNode(arguments.args[1])
    elseif isexpr(arguments.args[1]) && (arguments.args[1].head == :$)
        # interpolate
        # Now this becomes a Symbol
        assigned_var_qn = arguments.args[1].args[1] # ")  #QuoteNode(arguments.args[1].args[1]) #

    elseif isa(arguments.args[1].args[1], QuoteNode)
        #assigned_var = arguments.args[1].args[1].value
        # and the QuoteNode
        assigned_var_qn = arguments.args[1].args[1]
    else
        error("Assigned variable must be referenced using the colon syntax, e.g. :x.")
    end

    transformation = arguments.args[2]
    # return `nothing`s in case the index is < 1
    isexpr(transformation) && replace_invalid_indices!(transformation)

    return esc(
        quote
            # check variable is not present
            ($(assigned_var_qn) ∉ propertynames($t)) || error("Variable $($(assigned_var_qn)) already present in DataFrame.")

            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, $varlist_sort)
            end

            #determine type of resulting column from the type of the first element
            assigned_var_type = eltype([Douglass.@with($t,$(transformation)) for _n=1])
            $t[!,$(assigned_var_qn)] = missings(assigned_var_type,size($t,1))

            Douglass.@with $t begin
                # define _N
                _N = size($t, 1)
                # fill the new variable, row by row
                for _n in 1:size($t,1)
                    if (isnothing($filter) ? true : $filter)  # if condition is not satisfied, leave with missing
                        $(assigned_var_qn)[_n] = $(transformation)
                        # This is what needs to be used if the assigned_var_qn is a Symbol
                        #$t[_n, $(assigned_var_qn)] = $(transformation)
                    end
                end
                $t
            end
            $t
        end
    )
end

# # main version, doing stuff by row using @with (which has the advantage that we can use "[i]" syntax)
# # todo:
# #   allow _n
# #   assume that user means [i] if no index is shown (in transformation and filter)
# macro generate_byrow!(t::Symbol, varlist_by::Expr, varlist_sort::Expr, assigned_var, assigned_var_type, transformation::Expr, filter)
#     # @show assigned_var
#     # @show typeof(assigned_var)
#     # @show transformation
#     esc(
#         quote
#             # check that assigned_var_type is a valid type
#             isa($(assigned_var_type),DataType) || error("assigned_var_type must be a DataType")

#             # check variable is not present
#             ($(assigned_var) ∉propertynames($t)) || error("Variable $(assigned_var) already present in DataFrame.")
#             $t[!,$(assigned_var)] = missings($assigned_var_type,size($t,1))

#             # sort, if we need to, (first by-variables, then sort-variables)
#             if !isempty($(varlist_sort))
#                 sort!($t, vcat($varlist_by, $varlist_sort))
#             end

#             # this is the function that maps every sub-df into its transformed df
#             my_f = _df -> Douglass.@with _df begin
#                 # fill the new variable, row by row
#                 for i in 1:size(_df,1)
#                     if ($filter)  # if condition is not satisfied, leave with missing
#                         $(assigned_var)[i] = $(transformation)
#                     end
#                 end
#                 _df
#             end
#             # execute, and copy it to another dataframe,
#             # otherwise we get copies of the group variables in there as well
#             t2 = by($t, $varlist_by, my_f )
#             $t[!,$(assigned_var)] = t2[!,$(assigned_var)]
#             $t
#         end
#     )
# end

# macro generate!(t,varname,ex,filter)
#     esc(
#         quote
#             if $varname ∈ propertynames($t)
#                 error("Table already has a column with this name.")
#             end
#             local x = Douglass.@with($t, ifelse.($filter, $ex , missing))
#             # if we have a scalar, broadcast
#             if size(x,1) == 1
#                 $t[!,$varname] .= x
#             else
#                 $t[!,$varname] = x
#             end
#         end
#     )
# end

# @generate! without if-filter
# macro generate!(t,varname,ex)
#     return esc( :( Douglass.@generate!($t,$varname,$ex,true) )  )
# end

# # this is an alternative version that uses @byrow! Has the advantage that we don't need to use vectorized syntax, but not much else
# macro generate_byrow2!(t::Symbol, varlist_by::Expr, varlist_sort::Expr, assigned_var, assigned_var_type::Expr, transformation::Expr, filter::Expr, arguments::Expr)
#     esc(
#         quote
#             # create arguments in a local scope
#             args = $arguments

#             # check variable is not present
#             ($(assigned_var) ∉ propertynames($t)) || error("Variable $(assigned_var) already present in DataFrame.")

#             # this is the function that maps every sub-df into its transformed df
#             my_f = _df -> @byrow! _df begin
#                 @newcol $assigned_var::Array{Union{$assigned_var_type, Missing},1}
#                 $(assigned_var) = $(transformation)
#             end
#             $t = by($t, $varlist_by, my_f )
#             $t
#         end
#     )
# end


text:
"""
The geometry consists of a “hook” in the form of a curved strip rigidly clamped
at one end and loaded with a unit in-plane shear along the width at the other
end. It has two circular segments that are connected at the tangent point. The
smaller segment has a mean radius of 0.3556 m (14 inches) and spans 60° from
the clamped end to the tangent point. The larger segment spans 150° from the
tangent point to the free end and has a mean radius of 1.1684 m (46 inches).
The hook is 0.0508 m (2 inches) thick and 0.508 m (20 inches) wide, modeled as
linear elastic with an elastic modulus of 22.77 MPa (3300 psi) and a Poisson's
ratio of 0.35. In most tests the shear force is applied through the use of a
distributing coupling constraint. The coupling constraint provides coupling
between a reference node on which the load is prescribed and the nodes located
on the free end. The distributed nodal loads on the free end are equivalent to
a uniformly distributed load of 8.7563 N/m (0.05 lb/in). In two of the tests an
equivalent shear force is applied as a distributed shear traction instead.

"""
module raasch_examples

using LinearAlgebra, Statistics
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule: num_normals
using TestT3FF.FEMMShellT3DSGMTModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

using Infiltrator

function _execute(input = "raasch_s4_1x9.inp", drilling_stiffness_scale = 1.0, visualize = true, nL = 9, nW = 1)
    E = 3300.0;
    nu = 0.35;
    thickness  =  2.0;
    tolerance = thickness/2
    # analytical solution for the vertical deflection under the load
    bench_sol = 5.022012648671993; # 20-node hex results
    # bench_sol = 4.82482; # Ko et al, Performance ...
    R = 46.0;
    formul = FEMMShellT3FFModule
    # formul = FEMMShellT3DSGMTModule
    @show formul

    if input == ""
        fens, fes = T3block(210.0, 20.0, nL, nW)
        fens.xyz = xyz3(fens)
        for k in 1:count(fens)
            a, z = fens.xyz[k, :]
            if a <= 60.0
                fens.xyz[k, :] .= (14*sin(a*pi/180), 14*(1-cos(a*pi/180)), z)
            else
                a = a - 60.0 + 30.0
                xc = 14*sin(60*pi/180) + 46*cos(30*pi/180)
                yc = 14*(1-cos(60*pi/180)) - 46*sin(30*pi/180)
                fens.xyz[k, :] .= (xc-46*cos(a*pi/180), yc+46*sin(a*pi/180), z)
            end
        end
    else
        output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
        fens = output["fens"]
        fes = output["fesets"][1]

        connected = findunconnnodes(fens, fes);
        fens, new_numbering = compactnodes(fens, connected);
        fes = renumberconn!(fes, new_numbering);

        fens, fes = Q4toT3(fens, fes)
    # fens, fes = T3refine(fens, fes)# .
    # fens, fes = T3refine(fens, fes)# .
    # fens, fes = T3refine(fens, fes)# .
    end

    @show count(fens), count(fes)

    # plots = cat(plot_space_box([[0 0 -R/2]; [R/2 R/2 R/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

        # Report

    @info "Mesh: $input"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.drilling_stiffness_scale = drilling_stiffness_scale
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped end
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    @show num_normals(femm)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = [97.96152422706632 97.96152422706632 -16 -16 0 20], inflate = 1.0e-6)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, 0.05, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);
    @assert  isapprox(sum(F), 1.0)

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    nl = selectnode(fens; box = Float64[97.96152422706632 97.96152422706632 -16 -16 0 20], inflate = 1.0e-6)
    targetu =  mean(dchi.values[nl, 3])
    @info "Solution: $(round(targetu, digits=8)),  $(round(targetu/bench_sol, digits = 4)*100)%"

    # Visualization
    if visualize
        scattersysvec!(dchi, (R/2)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -R]; [R R R]]),
            plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end
    return targetu/bench_sol
end

function test_convergence()

    @info "Raasch hook"

    # for m in ["1x9", "3x18", "5x36", "10x72", "20x144"]
        # _execute("raasch_s4_" * m * ".inp", 1.0, false)
    for n in 1:7
        _execute("", 1.0, false, 9*2^(n-1), 1*2^(n-1))
    end
    return true
end

function test_dep_drilling_stiffness_scale()

    @info "Raasch hook"

    all_results = []
    all_drilling_stiffness_scale = [1000.0, 1.0, 0.1, 0.0001, 0.000001]
    for drilling_stiffness_scale in all_drilling_stiffness_scale
        results = Float64[]
        # for m in ["1x9", "3x18", "5x36", "10x72", "20x144"]
            # v = _execute("raasch_s4_" * m * ".inp", drilling_stiffness_scale, false)
        for n in 3:7
            v = _execute("", drilling_stiffness_scale, false, 9*2^(n-1), 1*2^(n-1))
                    push!(results, v)
        end
        push!(all_results, results)
    end
    return all_drilling_stiffness_scale, all_results
end

end # module

using .raasch_examples
raasch_examples.test_convergence()

all_drilling_stiffness_scale, all_results = raasch_examples.test_dep_drilling_stiffness_scale()

using PGFPlotsX

objects = []

ns = [20, 76, 222, 803]
es = [18, 57, 72, 288, 1152, 4608, 18432, 73728, 294912, 1179648]

styles = ["solid", "dashed", "dotted", "dashdotted", "densely dotted"]

let
    for (drilling_stiffness_scale, results, style) in zip(all_drilling_stiffness_scale, all_results, styles)
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        line_width  = 0.7,
        style = style,
        },
        Coordinates([v for v in  zip(es, results)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("$drilling_stiffness_scale"))
    end

    @pgf ax = Axis(
    {
    xlabel = "Number of Elements [ND]",
    ylabel = "Normalized Displacement [ND]",
    ymin = 0.95,
    ymax = 1.0,
    xmode = "log",
    ymode = "linear",
    yminorgrids = "true",
    grid = "both",
    legend_style = {
    at = Coordinate(0.5, 1.05),
    anchor = "south",
    legend_columns = -1
    },
    },
    objects...
    )

    display(ax)
    pgfsave("raasch_examples-dependence-on-drilling_stiffness_scale.pdf", ax)
end


text:
using LinearAlgebra: sorteig!, BlasReal, BlasComplex, StridedMatrix, eigsortby

const liblapack = Base.liblapack_name

import LinearAlgebra.BLAS.@blasfunc
using Base: iszero, require_one_based_indexing
using LinearAlgebra.LAPACK: chklapackerror
using LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare

# function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
    # alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
    # return GeneralizedEigen(sorteig!(alpha./beta, vr, sortby)...)
# end

for (ggev, elty, relty) in
    ((:zggev_,:ComplexF64,:Float64),
     (:cggev_,:ComplexF32,:Float32))
    @eval begin
        # SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
        #      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBVL, JOBVR
        #       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
        # *     ..
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION   RWORK( * )
        #       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
        #      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
        #      $                   WORK( * )
	function lapack_ggev!(A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, alpha::AbstractVector{$elty}, beta::AbstractVector{$elty}, vl::AbstractMatrix{$elty}, vr::AbstractMatrix{$elty}, work::AbstractVector{$elty}, rwork::AbstractVector{$relty})
            jobvl, jobvr, = 'N', 'V'
            require_one_based_indexing(A, B)
            chkstride1(A, B)
            n, m = checksquare(A, B)
            if n != m
                throw(DimensionMismatch("A has dimensions $(size(A)), and B has dimensions $(size(B)), but A and B must have the same size"))
            end
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            # alpha = similar(A, $elty, n)
            # beta = similar(A, $elty, n)
	    if ((size(alpha, 1) != n) || (size(beta, 1) != n))
               throw(DimensionMismatch("alpha has dimensions $(size(alpha)), and beta has dimensions $(size(beta)), but A and B must have size $(n)"))
            end
            ldvl = 0
            if jobvl == 'V'
                ldvl = n
            elseif jobvl == 'N'
                ldvl = 1
            else
                throw(ArgumentError("jobvl must be 'V' or 'N', but $jobvl was passed"))
            end
	    if (size(vl) != (ldvl, n))
                throw(DimensionMismatch("vl has dimensions $(size(vl)),  must have size ($(ldvl), $(n)"))
            end
            # vl = similar(A, $elty, ldvl, n)
            ldvr = 0
            if jobvr == 'V'
                ldvr = n
            elseif jobvr == 'N'
                ldvr = 1
            else
                throw(ArgumentError("jobvr must be 'V' or 'N', but $jobvr was passed"))
            end
	    if (size(vl) != (ldvl, n))
                throw(DimensionMismatch("vl has dimensions $(size(vl)),  must have size ($(ldvl), $(n)"))
            end
            # vr = similar(A, $elty, ldvr, n)
            # work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            # rwork = Vector{$relty}(undef, 8n)
	    if (size(rwork, 1) != 8*n)
                throw(DimensionMismatch("rwork has dimensions $(size(rwork)),  must have size 8n $(8*n)"))
            end
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($ggev), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                     Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                     Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                     Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$relty},
                     Ptr{BlasInt}),
                    jobvl, jobvr, n, A,
                    lda, B, ldb, alpha,
                    beta, vl, ldvl, vr,
                    ldvr, work, lwork, rwork,
                    info)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                end
            end
            alpha, beta, vl, vr
        end
    end
end

N = 3
A, B = rand(ComplexF64, N, N), rand(ComplexF64, N, N)
alpha, beta = zeros(ComplexF64, N), zeros(ComplexF64, N)
vl, vr = zeros(ComplexF64, 1, N), zeros(ComplexF64, N, N)
work, rwork = zeros(ComplexF64, 1), zeros(Float64, 8*N)
display(lapack_ggev!(A, B, alpha, beta, vl, vr, work, rwork))


text:
const _SUPPORTED_FORMATS = ("csv",)

mutable struct SimulationResultsExport
    problems::Dict{String, ProblemResultsExport}
    start_time::Dates.DateTime
    end_time::Dates.DateTime
    path::Union{Nothing, String}
    format::String
end

function SimulationResultsExport(
    problems::Vector{ProblemResultsExport},
    params::SimulationStoreParams;
    start_time = nothing,
    end_time = nothing,
    path = nothing,
    format = "csv",
)
    # This end time is outside the bounds of the simulation.
    sim_end_time = params.initial_time + params.step_resolution * params.num_steps

    if start_time === nothing
        start_time = params.initial_time
    elseif start_time < params.initial_time || start_time >= sim_end_time
        throw(IS.InvalidValue("invalid start_time: $start_time"))
    end

    if end_time === nothing
        # Reduce the end_time to be within the simulation.
        end_time = sim_end_time - Dates.Second(1)
    elseif end_time < params.initial_time || end_time >= sim_end_time
        throw(IS.InvalidValue("invalid end_time: $end_time"))
    end

    if !(format in list_supported_formats(SimulationResultsExport))
        throw(IS.InvalidValue("format = $format is not supported"))
    end

    return SimulationResultsExport(
        Dict(x.name => x for x in problems),
        start_time,
        end_time,
        path,
        format,
    )
end

function SimulationResultsExport(filename::AbstractString, params::SimulationStoreParams)
    if splitext(filename)[2] != ".json"
        throw(IS.InvalidValue("only JSON files are supported: $filename"))
    end

    return SimulationResultsExport(read_json(filename), params)
end

function SimulationResultsExport(data::AbstractDict, params::SimulationStoreParams)
    problems = Vector{ProblemResultsExport}()
    for problem in get(data, "problems", [])
        if !haskey(problem, "name")
            throw(IS.InvalidValue("problem data does not define 'name'"))
        end

        problem_export = ProblemResultsExport(
            problem["name"],
            Set(Symbol(x) for x in get(problem, "duals", Set{String}())),
            Set(Symbol(x) for x in get(problem, "parameters", Set{String}())),
            Set(Symbol(x) for x in get(problem, "variables", Set{String}())),
            get(problem, "optimizer_stats", false),
        )
        push!(problems, problem_export)
    end

    start_time = get(data, "start_time", nothing)
    if start_time isa AbstractString
        start_time = Dates.DateTime(start_time)
    end

    end_time = get(data, "end_time", nothing)
    if end_time isa AbstractString
        end_time = Dates.DateTime(end_time)
    end

    return SimulationResultsExport(
        problems,
        params;
        start_time = start_time,
        end_time = end_time,
        path = get(data, "path", nothing),
        format = get(data, "format", "csv"),
    )
end

function get_problem_exports(x::SimulationResultsExport, problem_name)
    if !haskey(x.problems, problem_name)
        throw(IS.InvalidValue("problem $problem_name is not stored"))
    end

    return x.problems[problem_name]
end

function get_export_file_type(exports::SimulationResultsExport)
    if exports.format == "csv"
        return CSV.File
    end

    throw(IS.InvalidValue("format not supported: $(exports.format)"))
end

list_supported_formats(::Type{SimulationResultsExport}) = ("csv",)

function should_export(exports::SimulationResultsExport, tstamp::Dates.DateTime)
    return tstamp >= exports.start_time && tstamp <= exports.end_time
end

function should_export_dual(exports::SimulationResultsExport, tstamp, problem, name)
    return _should_export(exports, tstamp, problem, :duals, name)
end

function should_export_parameter(exports::SimulationResultsExport, tstamp, problem, name)
    return _should_export(exports, tstamp, problem, :parameters, name)
end

function should_export_variable(exports::SimulationResultsExport, tstamp, problem, name)
    return _should_export(exports, tstamp, problem, :variables, name)
end

function _should_export(exports::SimulationResultsExport, tstamp, problem, field_name, name)
    if tstamp < exports.start_time || tstamp >= exports.end_time
        return false
    end

    problem_exports = get_problem_exports(exports, problem)
    return _should_export(problem_exports, field_name, name)
end


text:
## Physics:
##
## We simulate 2d cavity flow for compressible fluid,
## more precisely, for an ideal gas at constant temperature.
## Therefore, from the equation of state, we know that
## P / ρ = 1 (if we choose the units such that the constants
## terms are reduced to 1).
##
## The equations used are the continuity and the momentum equations
##
## The velocity in the x-direction on the top is constant,
## there are walls on the bottom, left and right edges.
##
## →→→→→→→→
## |      |
## |      |
## --------
##
## Initial conditions:
## u = u0 = 0.05 at y = 2
## v⃗ = (u, v) = 0⃗ everywhere else
## ρ = p = p0 everywhere (Atmospheric pressure)
##
## Boundary conditions:
## u = 1 at y = 2
## ∂p / ∂y = 0 at y = 0
## ρ = p = p0 at y = 2
## ∂p / ∂x = 0 at x = 0, 2

using PyPlot

## Time steps
duration = 10.0
Δt = 0.001
nsteps = convert(Int64, duration / Δt)

## Space grid [0, 2] × [0, 2]
Δx = 0.05
Δy = Δx
nx = convert(Int64, 2.0 / Δx)
ny = convert(Int64, 2.0 / Δy)

## Viscocity
ν = 0.1

## Atmospheric pressure
p0 = 1.0

## Velocity at the top
u0 = 0.05

u = zeros(nx, ny)
v = zeros(nx, ny)
ρ = zeros(nx, ny)

ρ[:, :] .= p0
u[:, ny] .= u0

function update!(ρnext, unext, vnext, ρwork, uwork, vwork,
        Δt, Δx, Δy, u0, ν, nx, ny)
    for j = 2 : ny - 1
        for i = 2 : nx - 1
            ## Auxiliary variables to avoid typos when
            ## working with partial derivatives
            ##
            ## NW--N--NE
            ##     |
            ## W---C---E
            ##     |
            ## SW--S--SE

            ρc = ρwork[i, j]
            ρe = ρwork[i + 1, j]
            ρw = ρwork[i - 1, j]
            ρn = ρwork[i, j + 1]
            ρs = ρwork[i, j - 1]

            uc = uwork[i, j]
            ue = uwork[i + 1, j]
            uw = uwork[i - 1, j]

            vc = vwork[i, j]
            vn = vwork[i, j + 1]
            vs = vwork[i, j - 1]

            ## Continuity equation to calculate ρ
            ## ∂ρ / ∂t + div(ρv⃗) = 0
            ρnext[i, j] = ρc - Δt * (uc * (ρe - ρw) / (2.0 * Δx) +
                                     vc * (ρn - ρs) / (2.0 * Δy) +
                                     ρc * (ue - uw) / (2.0 * Δx) +
                                     ρc * (vn - vs) / (2.0 * Δy))

        end
    end

    ## Question about view: https://stackoverflow.com/questions/65979923/what-is-happening-behind-the-scenes-in-julias-view-function-a3-viewa

    ## Enforcing boundary conditions for ρ = P
    ρnext[1, :] = @view ρnext[2, :]
    ρnext[nx, :] = @view ρnext[nx - 1, :]

    ρnext[:,  1] = @view ρnext[:, 2]
    ρnext[:, ny] .= p0

    ## Here we use a central difference scheme to discretize the second
    ## order derivatives, so we iterate over the 'inner' grid
    for j = 2 : ny - 1
        for i = 2 : nx - 1
            ## Auxiliar variables
            ρc = ρwork[i, j]
            ρe = ρwork[i + 1, j]
            ρw = ρwork[i - 1, j]
            ρn = ρwork[i, j + 1]
            ρs = ρwork[i, j - 1]

            uc = uwork[i, j]
            ue = uwork[i + 1, j]
            uw = uwork[i - 1, j]
            un = uwork[i, j + 1]
            us = uwork[i, j - 1]

            vc = vwork[i, j]
            ve = vwork[i + 1, j]
            vw = vwork[i - 1, j]
            vn = vwork[i, j + 1]
            vs = vwork[i, j - 1]

            ## Momentum equation
            ## ∂v⃗ / ∂t + (v⃗ ⋅∇)v⃗ = - (1 / ρ) ∇P + ν∇²v⃗
            ## -> we later replace this by Darcy's law
            viscous_term_u = ν * ((uw - 2.0 * uc + ue) / Δx^2 +
                                  (un - 2.0 * uc + us) / Δy^2)

            viscous_term_v = ν * ((vw - 2.0 * vc + ve) / Δx^2 +
                                  (vn - 2.0 * vc + vs) / Δy^2)

            pressure_term_u = -(1 / ρc) * (ρe - ρw) / (2.0 * Δx)

            pressure_term_v = -(1 / ρc) * (ρn - ρs) / (2.0 * Δy)

            inertia_term_u = uc * (ue - uw) / (2.0 * Δx) +
                             vc * (un - us) / (2.0 * Δy)

                             inertia_term_v = uc * (ve - vw) / (2.0 * Δx) +
                             vc * (vn - vs) / (2.0 * Δy- vw) / (2.0 * Δx)

            unext[i, j] = uc + Δt * (-inertia_term_u + pressure_term_u +
                                     viscous_term_u)

            vnext[i, j] = vc + Δt * (-inertia_term_v + pressure_term_v +
                                     viscous_term_v)
        end
    end
end

function cavity_flow!(ρ, u, v, Δx, Δy, Δt, nsteps, nx, ny, u0, ν)
    ## Helper arrays to store values from previous iteration
    ρwork, uwork, vwork = ρ, u, v
    unext = copy(u)
    vnext = copy(v)
    ρnext = copy(ρ)

    for t = 1 : nsteps
        update!(ρnext, unext, vnext, ρwork, uwork, vwork,
                Δt, Δx, Δy, u0, ν, nx, ny)

        ## Swap next and current and continue iteration
        ρnext, ρwork = ρwork, ρnext
        unext, uwork = uwork, unext
        vnext, vwork = vwork, vnext
    end

    ## If in the end ρwork is not the original ρ,
    ## copy the most recent data into ρ, u, v
    if ρwork !== ρ
        ρ .= ρwork
        u .= uwork
        v .= vwork
    end

    return ρ, u, v
end

function plot(u, v, ρ)
    fig = PyPlot.figure(figsize=(10, 10))
    ax = PyPlot.axes()
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    #legend(loc="upper right",fancybox="true")

    rgx = range(0, 2.0, length=nx)
    rgy = range(0, 2.0, length=ny)
    x = [rgx;]'
    y = [rgy;]'
    X = repeat([rgx;]', length(x))
    Y = repeat([rgy;],1, length(y))

    ## Velocity vector field
    ## scale = 100.0
    quiver(x, y, u', v')

    ## Streamplot of velocity vector field
    #streamplot(x, y, u', v')

    ## Density contour map
    pos = ax.contourf(X, Y, ρ', alpha=0.5, cmap=matplotlib.cm.viridis)
    fig.colorbar(pos, ax=ax)
    cp = contour(X, Y, ρ', cmap=matplotlib.cm.viridis)

    PyPlot.title("Cavity Flow for Ideal Gas")

    # PyPlot.svg(true)
    # savefig("plot.svg")

    savefig("img/ideal-gas.png")
    savefig("img/ideal-gas.svg")
end

## Run our solver
ρ, u, v = cavity_flow!(ρ, u, v, Δx, Δy, Δt, nsteps, nx, ny, u0, ν)

## Plotting
plot(u, v, ρ)

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
# Autor:        Juergen Vorloeper
# Organization: Hochschule Ruhr West
# Date:         Summer 2017

module HRW

export glpk
export solve_dgl
export GnuPlotOutput_Direct, GnuPlotOutput, GnuPlotOutput_Close

# load files
include("HRW_glpk.jl")
include("HRW_ode.jl")
include("HRW_gnuplot.jl")

# set up global variables, that stores gnuplot's state
gnuplot_state = GnuplotState(false,0,Any[])

# Register function "GnuPlotOutput_Close" to be called,
# when gnuplot_state goes out of scope
finalizer(gnuplot_state,GnuPlotOutput_Close)

end

max_stars_repo_name: j-vorloeper/HRW

max_issues_repo_name: j-vorloeper/HRW

max_forks_repo_name: j-vorloeper/HRW


text:
# Bar geometry summarizes data as vertical bars.
immutable BarGeometry <: Gadfly.GeometryElement
    # How bars should be positioned if they are grouped by color.
    # Valid options are:
    #   :stack -> place bars on top of each other (default)
    #   :dodge -> place bar next to each other
    position::Symbol

    # :vertical or :horizontal
    orientation::Symbol

    default_statistic::Gadfly.StatisticElement

    function BarGeometry(default_statistic=Gadfly.Stat.identity();
                         position::Symbol=:stack, orientation::Symbol=:vertical)
        new(position, orientation, default_statistic)
    end
end

const bar = BarGeometry

function histogram(; position=:stack, bincount=nothing,
                   minbincount=3, maxbincount=150,
                   orientation::Symbol=:vertical,
                   density::Bool=false)
    BarGeometry(Gadfly.Stat.histogram(bincount=bincount,
                                      minbincount=minbincount,
                                      maxbincount=maxbincount,
                                      position=position,
                                      orientation=orientation,
                                      density=density),
                position=position,
                orientation=orientation)
end

function element_aesthetics(::BarGeometry)
    [:x, :xmin, :xmax, :y, :ymin, :ymax, :color]
end

function default_statistic(geom::BarGeometry)
    return geom.default_statistic
end

# Render a single color bar chart
function render_colorless_bar(geom::BarGeometry,
                              theme::Gadfly.Theme,
                              aes::Gadfly.Aesthetics,
                              orientation::Symbol)
    if orientation == :horizontal
        ctx = compose!(
            context(),
            rectangle([0.0],
                      [Measure(cy=ymin) + theme.bar_spacing/2 for ymin in aes.ymin],
                      aes.x,
                      [Measure(cy=(ymax - ymin)) - theme.bar_spacing
                       for (ymin, ymax) in zip(aes.ymin, aes.ymax)]),
            svgclass("geometry"))
    else
        ctx = compose!(
            context(),
            rectangle([Measure(cx=xmin) + theme.bar_spacing/2 for xmin in aes.xmin],
                      [0.0],
                      [Measure(cx=(xmax - xmin)) - theme.bar_spacing
                       for (xmin, xmax) in zip(aes.xmin, aes.xmax)],
                      aes.y),
            svgclass("geometry"))
    end

    compose!(ctx, fill(theme.default_color))
    if isa(theme.bar_highlight, Function)
        compose!(ctx, stroke(theme.bar_highlight(theme.default_color)))
    elseif isa(theme.bar_highlight, ColorValue)
        compose!(ctx, stroke(theme.bar_highlight))
    else
        compose!(ctx, stroke(nothing))
    end
    return ctx
end

# Render a bar chart grouped by discrete colors and stacked.
function render_colorful_stacked_bar(geom::BarGeometry,
                                     theme::Gadfly.Theme,
                                     aes::Gadfly.Aesthetics,
                                     orientation::Symbol)
    # preserve factor orders of pooled data arrays
    if isa(aes.color, PooledDataArray)
        idxs = sortperm(aes.color.refs, rev=true)
    else
        idxs = 1:length(aes.color)
    end

    ctx = context()
    if orientation == :horizontal
        stack_height_dict = Dict()
        for y in aes.ymin
            stack_height_dict[y] = zero(eltype(aes.x))
        end
        stack_height = zeros(eltype(aes.x), length(idxs))

        for (i, j) in enumerate(idxs)
            stack_height[i] = stack_height_dict[aes.ymin[j]]
            stack_height_dict[aes.ymin[j]] += aes.x[j]
        end

        compose!(
            ctx,
            rectangle(
                stack_height,
                [aes.ymin[i]*cy - theme.bar_spacing/2 for i in idxs],
                [aes.x[i] for i in idxs],
                [(aes.ymax[i] - aes.ymin[i])*cy + theme.bar_spacing for i in idxs]))
    elseif orientation == :vertical
        stack_height_dict = Dict()
        for x in aes.xmin
            stack_height_dict[x] = zero(eltype(aes.y))
        end
        stack_height = zeros(eltype(aes.y), length(idxs))

        for (i, j) in enumerate(idxs)
            stack_height[i] = stack_height_dict[aes.xmin[j]]
            stack_height_dict[aes.xmin[j]] += aes.y[j]
        end

        compose!(
            ctx,
            rectangle(
                [aes.xmin[i]*cx + theme.bar_spacing/2 for i in idxs],
                stack_height,
                [(aes.xmax[i] - aes.xmin[i])*cx - theme.bar_spacing for i in idxs],
                [aes.y[i] for i in idxs]))
    else
        error("Orientation must be :horizontal or :vertical")
    end

    cs = [aes.color[i] for i in idxs]
    compose!(ctx, fill(cs), svgclass("geometry"))
    if isa(theme.bar_highlight, Function)
        compose!(ctx, stroke([theme.bar_highlight(c) for c in cs]))
    elseif isa(theme.bar_highlight, ColorValue)
        compose!(ctx, stroke(theme.bar_highlight))
    else
        compose!(ctx, stroke(nothing))
    end
    return ctx
end

# Render a bar chart grouped by discrete colors and stuck next to each other.
function render_colorful_dodged_bar(geom::BarGeometry,
                                    theme::Gadfly.Theme,
                                    aes::Gadfly.Aesthetics,
                                    orientation::Symbol)
    # preserve factor orders of pooled data arrays
    if isa(aes.color, PooledDataArray)
        idxs = sortperm(aes.color.refs, rev=true)
    else
        idxs = 1:length(aes.color)
    end

    ctx = context()
    if orientation == :horizontal
        dodge_count = DefaultDict(() -> 0)
        for i in idxs
            dodge_count[aes.ymin[i]] += 1
        end

        dodge_height = Dict()
        dodge_pos_dict = DefaultDict(() -> 0cy)
        for i in idxs
            dodge_height[aes.ymin[i]] =
                ((aes.ymax[i] - aes.ymin[i]) / dodge_count[aes.ymin[i]]) * cy
            dodge_pos_dict[aes.ymin[i]] = aes.ymin[i]*cy
        end

        dodge_pos = Array(Measure, length(idxs))
        for (i, j) in enumerate(idxs)
            dodge_pos[i] = dodge_pos_dict[aes.ymin[j]] - theme.bar_spacing/2
            dodge_pos_dict[aes.ymin[j]] += dodge_height[aes.ymin[j]]
        end

        compose!(
            ctx,
            rectangle(
                [0.0],
                dodge_pos,
                [aes.x[i] for i in idxs],
                [((aes.ymax[i] - aes.ymin[i])*cy + theme.bar_spacing) / dodge_count[aes.ymin[i]]
                 for i in idxs]))
    elseif orientation == :vertical
        dodge_count = DefaultDict(() -> 0)
        for i in idxs
            dodge_count[aes.xmin[i]] += 1
        end

        dodge_width = Dict()
        dodge_pos_dict = DefaultDict(() -> 0cx)
        for i in idxs
            dodge_width[aes.xmin[i]] =
                ((aes.xmax[i] - aes.xmin[i]) / dodge_count[aes.xmin[i]]) * cx
            dodge_pos_dict[aes.xmin[i]] = aes.xmin[i]*cx
        end

        dodge_pos = Array(Measure, length(idxs))
        for (i, j) in enumerate(idxs)
            dodge_pos[i] = dodge_pos_dict[aes.xmin[j]] + theme.bar_spacing/2
            dodge_pos_dict[aes.xmin[j]] += dodge_width[aes.xmin[j]]
        end

        compose!(
            ctx,
            rectangle(
                dodge_pos,
                [0.0],
                [((aes.xmax[i] - aes.xmin[i])*cx - theme.bar_spacing) / dodge_count[aes.xmin[i]]
                 for i in idxs],
                [aes.y[i] for i in idxs]))
    else
        error("Orientation must be :horizontal or :vertical")
    end

    cs = [aes.color[i] for i in idxs]
    compose!(ctx, fill(cs), svgclass("geometry"))
    if isa(theme.bar_highlight, Function)
        compose!(ctx, stroke([theme.bar_highlight(c) for c in cs]))
    elseif isa(theme.bar_highlight, ColorValue)
        compose!(ctx, stroke(theme.bar_highlight))
    else
        compose!(ctx, stroke(nothing))
    end
    return ctx
end

# Render bar geometry.
#
# Args:
#   geom: bar geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::BarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    if geom.orientation == :horizontal
        if (is(aes.ymin, nothing) || is(aes.ymax, nothing)) && is(aes.y, nothing)
            error("Geom.bar required \"y\" to be bound or both \"y_min\" and \"y_max\".")
        end
        if (aes.y != nothing && length(aes.y) != length(aes.x)) ||
           (aes.ymin != nothing && (length(aes.ymin) != length(aes.x) || length(aes.ymax) != length(aes.x)))
            error("Geom.bar requires x and y to be of equal length.")
        end

        var = :y
        minvar = :ymin
        maxvar = :ymax
    else
        if (is(aes.xmin, nothing) || is(aes.xmax, nothing)) && is(aes.x, nothing)
            error("Geom.bar required \"x\" to be bound or both \"x_min\" and \"x_max\".")
        end
        if (aes.x != nothing && length(aes.x) != length(aes.y)) &&
           (aes.xmin != nothing && (length(aes.xmin) != length(aes.y) || length(aes.xmax) != length(aes.y)))
            error("Geom.bar requires x and y to be of equal length.")
        end
        var = :x
        minvar = :xmin
        maxvar = :xmax
    end

    if getfield(aes, minvar) === nothing
        aes2 = Gadfly.Aesthetics()
        values = getfield(aes, var)
        minvalue, maxvalue = minimum(values), maximum(values)
        T = typeof((maxvalue - minvalue) / 1.0)

        span = zero(T)
        unique_count = length(Set(values))
        if unique_count > 1
            span = (maximum(values) - minimum(values)) / convert(Float64, (unique_count - 1))
        end

        if span == zero(T)
            span = one(T)
        end

        T = promote_type(eltype(values), typeof(span/2.0))
        setfield!(aes2, minvar, Array(T, length(values)))
        setfield!(aes2, maxvar, Array(T, length(values)))

        for (i, x) in enumerate(values)
            getfield(aes2, minvar)[i] = x - span/2.0
            getfield(aes2, maxvar)[i] = x + span/2.0
        end
        aes = inherit(aes, aes2)
    end

    if aes.color === nothing
        ctx = render_colorless_bar(geom, theme, aes, geom.orientation)
    elseif geom.position == :stack
        ctx = render_colorful_stacked_bar(geom, theme, aes, geom.orientation)
    elseif geom.position == :dodge
        ctx = render_colorful_dodged_bar(geom, theme, aes, geom.orientation)
    else
        error("$(geom.position) is not a valid position for the bar geometry.")
    end

    return compose!(
        context(),
        ctx,
        linewidth(theme.highlight_width),
        svgattribute("shape-rendering", "crispEdges"))
end


text:
module VSCodeServer

export vscodedisplay, @vscodedisplay, @enter, @run
export view_profile, @profview

using REPL, Sockets, Base64, Pkg, UUIDs, Dates, Profile
import Base: display, redisplay
import Dates
import Profile
import Logging
import InteractiveUtils

function __init__()
    atreplinit() do repl
        @async try
            hook_repl(repl)
        catch err
            Base.display_error(err, catch_backtrace())
        end
    end

    push!(Base.package_callbacks, on_pkg_load)

    for pkgid in keys(Base.loaded_modules)
        on_pkg_load(pkgid)
    end

    if VERSION >= v"1.4" && isdefined(InteractiveUtils, :EDITOR_CALLBACKS)
        pushfirst!(InteractiveUtils.EDITOR_CALLBACKS, function (cmd::Cmd, path::AbstractString, line::Integer)
            cmd == `code` || return false
            JSONRPC.send(conn_endpoint[], repl_open_file_notification_type, (; path = String(path), line = Int(line)))
            return true
        end)
    end
end

include("../../JSON/src/JSON.jl")
include("../../CodeTracking/src/CodeTracking.jl")

module IJuliaCore
using ..JSON
using Printf
import Base64

include("../../IJuliaCore/src/packagedef.jl")
end

module JSONRPC
import ..JSON
import UUIDs

include("../../JSONRPC/src/packagedef.jl")
end

module JuliaInterpreter
using ..CodeTracking

@static if VERSION >= v"1.6.0"
    include("../../JuliaInterpreter/src/packagedef.jl")
else
    include("../../../packages-old/JuliaInterpreter/src/packagedef.jl")
end
end

module DebugAdapter
import ..JuliaInterpreter
import ..JSON
import ..JSONRPC
import ..JSONRPC: @dict_readable, Outbound

include("../../DebugAdapter/src/packagedef.jl")
end

const conn_endpoint = Ref{Union{Nothing,JSONRPC.JSONRPCEndpoint}}(nothing)

include("../../../error_handler.jl")
include("repl_protocol.jl")
include("misc.jl")
include("trees.jl")
include("module.jl")
include("progress.jl")
include("eval.jl")
include("completions.jl")
include("repl.jl")
include("./tables/tableviewer.jl")
include("display.jl")
include("profiler.jl")
include("debugger.jl")
include("notebookdisplay.jl")
include("serve_notebook.jl")

function dispatch_msg(conn_endpoint, msg_dispatcher, msg, is_dev)
    if is_dev
        try
            JSONRPC.dispatch_msg(conn_endpoint[], msg_dispatcher, msg)
        catch err
            Base.display_error(err, catch_backtrace())
        end
    else
        JSONRPC.dispatch_msg(conn_endpoint[], msg_dispatcher, msg)
    end
end

function serve(args...; is_dev = false, crashreporting_pipename::Union{AbstractString,Nothing} = nothing)
    if !HAS_REPL_TRANSFORM[] && isdefined(Base, :active_repl)
        hook_repl(Base.active_repl)
    end

    @debug "connecting to pipe"
    conn = connect(args...)
    conn_endpoint[] = JSONRPC.JSONRPCEndpoint(conn, conn)
    @debug "connected"
    if EVAL_BACKEND_TASK[] === nothing
        start_eval_backend()
    end
    run(conn_endpoint[])
    @debug "running"

    @async try
        msg_dispatcher = JSONRPC.MsgDispatcher()

        msg_dispatcher[repl_runcode_request_type] = repl_runcode_request
        msg_dispatcher[repl_interrupt_notification_type] = repl_interrupt_request
        msg_dispatcher[repl_getvariables_request_type] = repl_getvariables_request
        msg_dispatcher[repl_getlazy_request_type] = repl_getlazy_request
        msg_dispatcher[repl_showingrid_notification_type] = repl_showingrid_notification
        msg_dispatcher[repl_loadedModules_request_type] = repl_loadedModules_request
        msg_dispatcher[repl_isModuleLoaded_request_type] = repl_isModuleLoaded_request
        msg_dispatcher[repl_getcompletions_request_type] = repl_getcompletions_request
        msg_dispatcher[repl_resolvecompletion_request_type] = repl_resolvecompletion_request
        msg_dispatcher[repl_startdebugger_notification_type] = (conn, params) -> repl_startdebugger_request(conn, params, crashreporting_pipename)
        msg_dispatcher[repl_startdebugger_notification_type] = (conn, params) -> repl_startdebugger_request(conn, params, crashreporting_pipename)
        msg_dispatcher[repl_toggle_plot_pane_notification_type] = toggle_plot_pane
        msg_dispatcher[repl_toggle_diagnostics_notification_type] = toggle_diagnostics
        msg_dispatcher[repl_toggle_progress_notification_type] = toggle_progress
        msg_dispatcher[cd_notification_type] = cd_to_uri
        msg_dispatcher[activate_project_notification_type] = activate_uri
        msg_dispatcher[repl_getdebugitems_request_type] = debugger_getdebugitems_request
        msg_dispatcher[repl_gettabledata_request_type] = get_table_data
        msg_dispatcher[repl_clearlazytable_notification_type] = clear_lazy_table

        @sync while conn_endpoint[] isa JSONRPC.JSONRPCEndpoint && isopen(conn)
            msg = JSONRPC.get_next_message(conn_endpoint[])

            if msg["method"] == repl_runcode_request_type.method
                @async dispatch_msg(conn_endpoint, msg_dispatcher, msg, is_dev)
            else
                dispatch_msg(conn_endpoint, msg_dispatcher, msg, is_dev)
            end
        end
    catch err
        if !isopen(conn) && is_disconnected_exception(err)
            # expected error
            @debug "remote closed the connection"
        else
            try
                global_err_handler(err, catch_backtrace(), crashreporting_pipename, "REPL")
            catch err
                @error "Error handler threw an error." exception = (err, catch_backtrace())
            end
        end
    finally
        @debug "JSONRPC dispatcher task finished"
    end
end

end  # module

max_stars_repo_name: TanookiToad/julia-vscode

max_issues_repo_name: TanookiToad/julia-vscode

max_forks_repo_name: TanookiToad/julia-vscode


text:
module Graph

using DocStringExtensions

import ..Types
import ..Delayed

"""
$(TYPEDSIGNATURES)

Given an delayed object, return an array of argument values.
"""
function delayedproperties(item::Types.AbstractDelayed)::Vector

    return [
        getproperty(item, propertyname)
            for propertyname in propertynames(item)
    ]

end
delayedproperties(a::Delayed.AbstractSwitch) = [a.choice, a.options...]
delayedproperties(a::Delayed.CategoricalIndex) = a.probabilities

"""
$(TYPEDSIGNATURES)

Depth-first search
Unrolls a graph into an array of all the nodes.
"""
dfs(space::Types.SPACE_TYPE)::Vector = dfs!(Types.AbstractDelayed[], space)

"""
$(TYPEDSIGNATURES)

Unrolls a nested space into a vector of all the Types.AbstractDelayed nodes.
"""
function dfs!(seq::Vector{Types.AbstractDelayed}, item::Types.AbstractDelayed)::Vector
    # For every input (arg) of the object, add those to the list too.
    for prop in delayedproperties(item)
        dfs!(seq, prop)
    end
    push!(seq, item)

    return seq
end
function dfs!(seq::Vector{Types.AbstractDelayed}, item::Dict)::Vector
    for v in values(item)
        dfs!(seq, v)
    end

    return seq
end
function dfs!(seq::Vector{Types.AbstractDelayed}, item::Union{Tuple, Vector})::Vector
    for v in item
        dfs!(seq, v)
    end

    return seq
end
dfs!(seq::Vector{Types.AbstractDelayed}, item::Any)::Vector = seq

function checklabel!(labels::Vector{Symbol}, node::Delayed.AbstractParam)::Nothing
    if node.label in labels
        throw(DomainError(node.label, "Your space has a duplicate label :$(node.label)"))
    end
    push!(labels, node.label)

    return nothing
end
checklabel!(labels::Vector{Symbol}, node::T where T)::Nothing = nothing

"""
$(TYPEDSIGNATURES)

Ensure the user hasn't submitted any duplicate labels in their space.
"""
function checkspace(space::Types.SPACE_TYPE)::Nothing
    labels = Symbol[]
    for node in Graph.dfs(space)
        checklabel!(labels, node)
    end

    return nothing
end

end # module Graph


text:
type CuDNNConvState
  inputs_desc   :: Vector{CuDNN.Tensor4dDescriptor}
  outputs_desc  :: Vector{CuDNN.Tensor4dDescriptor}
  conv_desc     :: Vector{CuDNN.ConvolutionDescriptor}
  filter_desc   :: CuDNN.FilterDescriptor
  bias_desc     :: CuDNN.Tensor4dDescriptor

  fwd_algorithm :: Int
  workspace     :: CUDA.CuPtr
  workspace_size:: Int

  bottom_offset :: Int
  top_offset    :: Int
  weight_offset :: Int
  bias_offset   :: Int
end

function setup_etc(backend::GPUBackend, layer::ConvolutionLayer, dtype, width, height, channels,
    batch_size, width_out, height_out, inputs)

  filter_desc = CuDNN.create_filter_descriptor(dtype, (layer.kernel[1], layer.kernel[2],
      div(channels,layer.n_group), div(layer.n_filter,layer.n_group)))
  bias_desc = CuDNN.create_tensor4d_descriptor(dtype, (1,1,div(layer.n_filter,layer.n_group),1))

  inputs_desc = Array(CuDNN.Tensor4dDescriptor, length(inputs))
  outputs_desc = Array(CuDNN.Tensor4dDescriptor, length(inputs))
  conv_desc = Array(CuDNN.ConvolutionDescriptor, length(inputs))
  for i = 1:length(inputs)
    inputs_desc[i] = CuDNN.create_tensor4d_descriptor(dtype, (width,height,div(channels,layer.n_group),batch_size),
        (1, width, width*height, width*height*channels))
    outputs_desc[i] = CuDNN.create_tensor4d_descriptor(dtype, (width_out, height_out, div(layer.n_filter,layer.n_group), batch_size),
        (1, width_out, width_out*height_out, width_out*height_out*layer.n_filter))
    conv_desc[i] = CuDNN.create_convolution_descriptor(inputs_desc[i], filter_desc, layer.pad,
        layer.stride, (1,1), CuDNN.CUDNN_CROSS_CORRELATION)
  end

  # NOTE: in ConvolutionLayer, we require the shape for all inputs to be the same
  fwd_algorithm = CuDNN.get_convolution_forward_algorithm(backend.cudnn_ctx, inputs_desc[1], filter_desc, conv_desc[1], outputs_desc[1],
      CuDNN.CUDNN_CONVOLUTION_FWD_PREFER_FASTEST, 0)
  workspace_size = CuDNN.get_convolution_forward_workspace_size(backend.cudnn_ctx, inputs_desc[1], filter_desc, conv_desc[1], outputs_desc[1],
      fwd_algorithm)
  if workspace_size == 0
    workspace = CUDA.CuPtr()
  else
    workspace = CUDA.cualloc(UInt8, workspace_size) # workspace_size is in bytes
  end

  bottom_offset = div(channels,layer.n_group) * height * width * sizeof(dtype)
  top_offset = div(layer.n_filter,layer.n_group) * height_out * width_out * sizeof(dtype)
  weight_offset = div(layer.n_filter,layer.n_group) * div(channels,layer.n_group) * layer.kernel[1] * layer.kernel[2] * sizeof(dtype)
  bias_offset = div(layer.n_filter,layer.n_group) * sizeof(dtype)

  etc = CuDNNConvState(inputs_desc, outputs_desc, conv_desc, filter_desc, bias_desc,
      fwd_algorithm, workspace, workspace_size,
      bottom_offset, top_offset, weight_offset, bias_offset)
  return etc
end
function shutdown_etc(backend::GPUBackend, state::ConvolutionLayerState)
  etc = state.etc
  CuDNN.destroy_filter_descriptor(etc.filter_desc)
  CuDNN.destroy_tensor4d_descriptor(etc.bias_desc)
  map(CuDNN.destroy_tensor4d_descriptor, etc.inputs_desc)
  map(CuDNN.destroy_tensor4d_descriptor, etc.outputs_desc)
  map(CuDNN.destroy_convolution_descriptor, etc.conv_desc)
  if etc.workspace_size != 0
    CUDA.free(etc.workspace)
  end
end

function forward(backend::GPUBackend, state::ConvolutionLayerState, inputs::Vector{Blob})
  alpha = one(eltype(inputs[1]))
  beta_accumulate = one(eltype(inputs[1]))
  beta_dont_accumulate = zero(eltype(inputs[1]))
  fwd_algorithm = state.etc.fwd_algorithm
  workspace_ptr = state.etc.workspace
  workspace_size = state.etc.workspace_size

  for i = 1:length(inputs)
    for g = 1:state.layer.n_group
      input_ptr = CuPtr(inputs[i].ptr.p + state.etc.bottom_offset * (g-1))
      output_ptr = CuPtr(state.blobs[i].ptr.p + state.etc.top_offset * (g-1))
      filter_ptr = CuPtr(state.filter.ptr.p + state.etc.weight_offset * (g-1))
      CuDNN.convolution_forward(backend.cudnn_ctx, alpha, state.etc.inputs_desc[i], input_ptr,
          state.etc.filter_desc, filter_ptr, state.etc.conv_desc[i],
          state.etc.outputs_desc[i], output_ptr, workspace_ptr, workspace_size, fwd_algorithm,
          beta_dont_accumulate)

      # bias
      CuDNN.add_tensor4d(backend.cudnn_ctx, CuDNN.CUDNN_ADD_SAME_C, alpha,
          state.etc.bias_desc, CuPtr(state.bias.ptr.p + state.etc.bias_offset * (g-1)),
          beta_accumulate, state.etc.outputs_desc[i], output_ptr)
    end
  end
end

function backward(backend::GPUBackend, state::ConvolutionLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  erase!(state.∇filter)
  erase!(state.∇bias)
  alpha = one(eltype(inputs[1]))
  beta_accumulate = one(eltype(inputs[1]))
  beta_dont_accumulate = zero(eltype(inputs[1]))
  for i = 1:length(inputs)
    bottom = inputs[i]
    top_diff = state.blobs_diff[i]

    for g = 1:state.layer.n_group
      if !state.frozen
        # gradient w.r.t. bias
        CuDNN.convolution_backward_bias(backend.cudnn_ctx, alpha,
            state.etc.outputs_desc[i], CuPtr(top_diff.ptr.p + state.etc.top_offset * (g-1)),
            beta_accumulate, state.etc.bias_desc, CuPtr(state.∇bias.ptr.p + state.etc.bias_offset * (g-1)))

        # gradient w.r.t. weights
        CuDNN.convolution_backward_filter(backend.cudnn_ctx, alpha,
            state.etc.inputs_desc[i], CuPtr(bottom.ptr.p + state.etc.bottom_offset * (g-1)),
            state.etc.outputs_desc[i], CuPtr(top_diff.ptr.p + state.etc.top_offset * (g-1)),
            state.etc.conv_desc[i],
            beta_accumulate, state.etc.filter_desc, CuPtr(state.∇filter.ptr.p + state.etc.weight_offset * (g-1)))
      end

      # gradient w.r.t. bottom data
      if isa(diffs[i], CuTensorBlob)
        CuDNN.convolution_backward_data(backend.cudnn_ctx, alpha,
            state.etc.filter_desc, CuPtr(state.filter.ptr.p + state.etc.weight_offset * (g-1)),
            state.etc.outputs_desc[i], CuPtr(top_diff.ptr.p + state.etc.top_offset * (g-1)),
            state.etc.conv_desc[i],
            beta_dont_accumulate, state.etc.inputs_desc[i], CuPtr(diffs[i].ptr.p + state.etc.bottom_offset * (g-1)))
      end
    end
  end
end


text:
__precompile__()

module Snappy

using snappy_jll

export compress, uncompress

# snappy status
const SnappyOK             = Cint(0)
const SnappyInvalidInput   = Cint(1)
const SnappyBufferTooSmall = Cint(2)

# High-level Interfaces

function compress(input::Vector{UInt8})
    ilen = length(input)
    maxlen = snappy_max_compressed_length(UInt(ilen))
    compressed = Array{UInt8}(undef, maxlen)
    olen, st = snappy_compress(input, compressed)
    if st != SnappyOK
        error("failed to compress the data")
    end
    resize!(compressed, olen)
    compressed
end

function uncompress(input::Array{UInt8})
    ilen = length(input)
    explen, st = snappy_uncompressed_length(input)
    if st != SnappyOK
        error("faield to guess the length of the uncompressed data (the compressed data may be broken?)")
    end
    uncompressed = Array{UInt8}(undef, explen)
    olen, st = snappy_uncompress(input, uncompressed)
    if st != SnappyOK
        error("failed to uncompress the data")
    end
    @assert explen == olen
    resize!(uncompressed, olen)
    uncompressed
end

# Low-level Interfaces

function snappy_compress(input::Vector{UInt8}, compressed::Vector{UInt8})
    ilen = length(input)
    olen = Ref{Csize_t}(length(compressed))
    status = ccall(
        (:snappy_compress, libsnappy),
        Cint,
        (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Ref{Csize_t}),
        input, ilen, compressed, olen
    )
    olen[], status
end

function snappy_uncompress(compressed::Vector{UInt8}, uncompressed::Vector{UInt8})
    ilen = length(compressed)
    olen = Ref{Csize_t}(length(uncompressed))
    status = ccall(
        (:snappy_uncompress, libsnappy),
        Cint,
        (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Ref{Csize_t}),
        compressed, ilen, uncompressed, olen
    )
    olen[], status
end

function snappy_max_compressed_length(source_length::UInt)
    ccall((:snappy_max_compressed_length, libsnappy), Csize_t, (Csize_t,), source_length)
end

function snappy_uncompressed_length(compressed::Vector{UInt8})
    len = length(compressed)
    result = Ref{Csize_t}(0)
    status = ccall((:snappy_uncompressed_length, libsnappy), Cint, (Ptr{UInt8}, Csize_t, Ref{Csize_t}), compressed, len, result)
    result[], status
end

function snappy_validate_compressed_buffer(compressed::Vector{UInt8})
    ilen = length(compressed)
    ccall((:snappy_validate_compressed_buffer, libsnappy), Cint, (Ptr{UInt8}, Csize_t), compressed, ilen)
end

end # module


text:
using LinearAlgebra

export canonicalcor
function canonicalcor(x::AbstractMatrix, y::AbstractMatrix)
    ma = inv(cov(x))*cov(x, y)*inv(cov(y))*cov(y,x)
    mb = inv(cov(y))*cov(y, x)*inv(cov(x))*cov(x,y)
    evx = eigvecs(ma)
    evy = eigvecs(mb)
    abs(cor(x*evx[:, end], y*evy[:, end]))
    #[-cor(x*evx, y*evy) for (evx, evy) in zip(eachcol(evx), eachcol(evy))]
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
#####
##### `AbstractDifferential`
#####

"""
The subtypes of `AbstractDifferential` define a custom \"algebra\" for chain
rule evaluation that attempts to factor various features like complex derivative
support, broadcast fusion, zero-elision, etc. into nicely separated parts.

All subtypes of `AbstractDifferential` implement the following operations:

`+(a, b)`: linearly combine differential `a` and differential `b`

`*(a, b)`: multiply the differential `a` by the differential `b`

`Base.conj(x)`: complex conjugate of the differential `x`

`extern(x)`: convert `x` into an appropriate non-`AbstractDifferential` type for
use outside of `ChainContext`.

Valid arguments to these operations are `T` where `T<:AbstractDifferential`, or
where `T` has proper `+` and `*` implementations.

Additionally, all subtypes of `AbstractDifferential` support `Base.iterate` and
`Base.Broadcast.broadcastable(x)`.
"""
abstract type AbstractDifferential end

Base.:+(x::AbstractDifferential) = x

"""
    extern(x)

Return `x` converted to an appropriate non-`AbstractDifferential` type, for use
with external packages that might not handle `AbstractDifferential` types.

Note that this function may return an alias (not necessarily a copy) to data
wrapped by `x`, such that mutating `extern(x)` might mutate `x` itself.
"""
@inline extern(x) = x

@inline Base.conj(x::AbstractDifferential) = x

"""
    refine_differential(𝒟::Type, der)

Converts, if required, a differential object `der`
(e.g. a `Number`, `AbstractDifferential`, `Matrix`, etc.),
to another  differential that is more suited for the domain given by the type 𝒟.
Often this will behave as the identity function on `der`.
"""
refine_differential(::Any, der) = der  # most of the time leave it alone.


text:
@testset "Maximum bipartite matching" begin
    @testset "Interface" begin
        n = 5

        # Error: the graph is not bipartite.
        graph = complete_graph(n) # Requires n > 2
        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:n
            for j in 1:n
                if i != j
                    rewards[Edge(i, j)] = 1.0
                end
            end
        end
        @test_throws ErrorException BipartiteMatchingInstance(graph, rewards)

        # Error: edges to the same side of the bipartite graph.
        # Both in the left part.
        graph = complete_bipartite_graph(n, n)
        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:n
            for j in 1:n
                rewards[Edge(i, j)] = 1.0
            end
        end
        i = BipartiteMatchingInstance(graph, rewards)
        @test_throws ErrorException matching_hungarian(i)

        # Both in the right part.
        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:n
            for j in 1:n
                rewards[Edge(i + n, j + n)] = 1.0
            end
        end
        i = BipartiteMatchingInstance(graph, rewards)
        @test_throws ErrorException matching_hungarian(i)

        # Automatic handling: edges are not from V1 to V2 in the bipartite graph.
        rewards = Dict{Edge{Int}, Float64}()
        weights = Dict{Edge{Int}, Int}()
        for i in 1:n
            for j in 1:n
                rewards[Edge(j, i + n)] = 1.0
            end
        end
        i = BipartiteMatchingInstance(graph, rewards)
        @test matching_hungarian(i) != π # Just a @test_nothrows: https://github.com/JuliaLang/julia/issues/18780

        # Test the partial-copy constructor.
        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:n
            for j in 1:n
                rewards[Edge(j, i + n)] = 2.0
            end
        end
        i2 = BipartiteMatchingInstance(i, rewards)
        @test i.graph == i2.graph
        @test i.reward != i2.reward # Only difference to be found.
        @test i2.reward == rewards
        @test i.partition == i2.partition
        @test i.n_left == i2.n_left
        @test i.n_right == i2.n_right
        @test i.vertex_left == i2.vertex_left
        @test i.vertex_right == i2.vertex_right
    end

    @testset "Basic: 3×3" begin
        # Three nodes on each side (1-2-3 to the left, 4-5-6 to the right),
        # optimum is matching them in order (1-4, 2-5, 3-6)
        graph = complete_bipartite_graph(3, 3)
        rewards = Dict(
            Edge(1, 4) => 1.0,
            Edge(1, 5) => 0.0,
            Edge(1, 6) => 0.0,
            Edge(2, 4) => 0.0,
            Edge(2, 5) => 1.0,
            Edge(2, 6) => 0.0,
            Edge(3, 4) => 0.0,
            Edge(3, 5) => 0.0,
            Edge(3, 6) => 1.0,
        )

        i = BipartiteMatchingInstance(graph, rewards)

        s1 = matching_hungarian(i)
        s2 = matching_dp_imperfect(i)
        s3 = matching_dp(i)

        for s in [s1, s2, s3]
            @test s.instance == i
            @test length(s.solution) == 3
            @test Edge(1, 4) in s.solution
            @test Edge(2, 5) in s.solution
            @test Edge(3, 6) in s.solution
            @test s.value ≈ 3.0
        end
    end

    @testset "Conformity: 5×5" begin
        # Just a larger graph. Size 5 is what is used for conformity tests in budgeted matchings.
        n = 5
        graph = complete_bipartite_graph(n, n)

        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:n
            for j in 1:n
                rewards[Edge(i, j + n)] = 1.0
            end
        end

        i = BipartiteMatchingInstance(graph, rewards)

        s1 = matching_hungarian(i)
        s2 = matching_dp_imperfect(i)
        s3 = matching_dp(i)

        for s in [s1, s2, s3]
            @test s.instance == i
            @test length(s.solution) == n
            for i in 1:n
                @test Edge(i, i + n) in s.solution
            end
            @test s.value ≈ n
        end

        # Ensure one node has no matching in V1.
        graph = complete_bipartite_graph(n - 1, n)

        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:(n - 1)
            for j in 1:n
                rewards[Edge(i, j + n - 1)] = 1.0
            end
        end

        i = BipartiteMatchingInstance(graph, rewards)

        s1 = matching_hungarian(i)
        s2 = matching_dp_imperfect(i)
        s3 = matching_dp(i)

        for s in [s1, s2, s3]
            @test s.instance == i
            @test length(s.solution) == n - 1
            @test s.value ≈ n - 1.0
        end

        # Ensure one node has no matching in V2.
        graph = complete_bipartite_graph(n, n - 1)

        rewards = Dict{Edge{Int}, Float64}()
        for i in 1:n
            for j in 1:(n - 1)
                rewards[Edge(i, j + n)] = 1.0
            end
        end

        i = BipartiteMatchingInstance(graph, rewards)

        s1 = matching_hungarian(i)
        s2 = matching_dp_imperfect(i)
        s3 = matching_dp(i)

        for s in [s1, s2, s3]
            @test s.instance == i
            @test length(s.solution) == n - 1
            @test s.value ≈ n - 1.0
        end
    end
end

@testset "Budgeted maximum bipartite matching" begin
    @testset "Interface" begin
        # Graph must be bipartite.
        graph = complete_graph(3)
        rewards = Dict(Edge(1, 2) => 1.0, Edge(1, 3) => 0.5, Edge(2, 3) => 3.0)
        weights = Dict(Edge(1, 2) => 0, Edge(1, 3) => 2, Edge(2, 3) => 0)
        @test_throws ErrorException BudgetedBipartiteMatchingInstance(
            graph,
            rewards,
            weights,
            0,
        )

        # Values of ζ.
        graph = complete_bipartite_graph(2, 2)
        rewards = Dict(
            Edge(1, 3) => 1.0,
            Edge(1, 4) => 0.0,
            Edge(2, 3) => 0.0,
            Edge(2, 4) => 1.0,
        )
        weights = Dict(
            Edge(1, 3) => 0,
            Edge(1, 4) => 1,
            Edge(2, 3) => 1,
            Edge(2, 4) => 0,
        )
        i = BudgetedBipartiteMatchingInstance(graph, rewards, weights, 0)

        @test_throws ErrorException matching_hungarian_budgeted_lagrangian_refinement(
            i,
            ζ⁻=1.0,
        )
        @test_throws ErrorException matching_hungarian_budgeted_lagrangian_refinement(
            i,
            ζ⁻=2.0,
        )
        @test_throws ErrorException matching_hungarian_budgeted_lagrangian_refinement(
            i,
            ζ⁺=1.0,
        )
        @test_throws ErrorException matching_hungarian_budgeted_lagrangian_refinement(
            i,
            ζ⁺=0.2,
        )

        # Infeasible solutions.
        s1 =
            SimpleBudgetedBipartiteMatchingSolution(i) ==
            SimpleBudgetedBipartiteMatchingSolution(i, edgetype(graph)[])
    end

    @testset "Helpers" begin
        @test Kombinator._edge_any_end_match(Edge(1, 2), Edge(1, 2))
        @test Kombinator._edge_any_end_match(Edge(1, 2), Edge(1, 3))
        @test Kombinator._edge_any_end_match(Edge(1, 2), Edge(4, 2))
        @test Kombinator._edge_any_end_match(Edge(1, 3), Edge(1, 2))
        @test Kombinator._edge_any_end_match(Edge(4, 2), Edge(1, 2))
        @test !Kombinator._edge_any_end_match(Edge(1, 2), Edge(3, 4))
    end

    @testset "Basic: 2×2" begin
        # Two nodes on each side (1-2 to the left, 3-4 to the right),
        # unbudgeted optimum is matching them in order (1-3, 2-4).
        # Two solutions: 1-3 and 2-4 (reward: 2; weight: 0) or 1-4 and 2-3 (reward: 0; weight: 2).
        graph = complete_bipartite_graph(2, 2)
        rewards = Dict(
            Edge(1, 3) => 1.0,
            Edge(1, 4) => 0.0,
            Edge(2, 3) => 0.0,
            Edge(2, 4) => 1.0,
        )
        weights = Dict(
            Edge(1, 3) => 0,
            Edge(1, 4) => 1,
            Edge(2, 3) => 1,
            Edge(2, 4) => 0,
        )
        ε = 0.0001

        @testset "No budget constraint" begin
            # TODO: for a zero budget, always rely on standard matching algorithm? Could be much more efficient.
            budget = 0
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value ≈ 2.0 atol = ε # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 2
                @test Edge(1, 3) in s.solution
                @test Edge(2, 4) in s.solution
                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 2.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
            end
        end

        @testset "Loose budget constraint" begin
            budget = 1 # The only solution respecting this budget has a total weight of 2.
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value ≈ 1.0 atol = ε # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 2
                @test Edge(1, 4) in s.solution
                @test Edge(2, 3) in s.solution
                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 2.0 atol = ε
            end
        end

        @testset "Tight budget constraint" begin
            budget = 2
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value ≈ 0.0 atol = ε # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 2
                @test Edge(1, 4) in s.solution
                @test Edge(2, 3) in s.solution
                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 2.0 atol = ε
            end
        end

        @testset "Impossibly tight budget constraint" begin
            budget = 3 # Infeasible, but the Lagrangian relaxation does not really care.
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value <= 0.0 # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
            end

            # The Lagrangian relaxation does not care.
            @test s1.instance == i
            @test length(s1.solution) == 2

            # The other approximation algorithms correctly return an infeasible solution.
            for s in [s2, s3, s4]
                @test length(s.solution) == 0
                @test s.value == -Inf
            end
        end
    end

    @testset "Conformity: 3×3" begin
        # Three nodes on each side (1-2-3 to the left, 4-5-6 to the right),
        # optimum is matching them in order (1-4, 2-5, 3-6)
        graph = complete_bipartite_graph(3, 3)
        rewards = Dict{Edge{Int}, Float64}()
        weights = Dict{Edge{Int}, Int}()
        for i in 1:3
            for j in 1:3
                k = j + 3
                rewards[Edge(i, k)] = (i == j) ? 1.0 : 0.0
                weights[Edge(i, k)] = (i == j) ? 0 : 1
            end
        end
        ε = 0.0001

        @testset "No budget contraint" begin
            budget = 0
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value ≈ 3.0 atol = ε # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 3
                @test Edge(1, 4) in s.solution
                @test Edge(2, 5) in s.solution
                @test Edge(3, 6) in s.solution
                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 3.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
            end
        end

        @testset "Loose budget constraint" begin
            budget = 2
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value ≈ 1.0 atol = ε # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 3
                @test !(Edge(1, 4) in s.solution)
                @test !(Edge(2, 5) in s.solution)
                @test !(Edge(3, 6) in s.solution)
                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 3.0 atol = ε
            end
        end
    end

    # Not really useful to test 4×4, as everything would be fixed by the major loops in the 1/2 scheme.

    @testset "Conformity: 5×5" begin
        # Larger graph, to test more parts of the approximation scheme (fixing four edges in the 1/2 scheme, so take
        # at least one more so that the underlying additive approximation scheme has some work to do).
        graph = complete_bipartite_graph(5, 5)
        rewards = Dict{Edge{Int}, Float64}()
        weights = Dict{Edge{Int}, Int}()
        for i in 1:5
            for j in 1:5
                k = j + 5
                rewards[Edge(i, k)] = (i == j) ? 1.0 : 0.0
                weights[Edge(i, k)] = (i == j) ? 0 : 1
            end
        end
        ε = 0.0001
        max_reward_edges =
            [Edge(1, 6), Edge(2, 7), Edge(3, 8), Edge(4, 9), Edge(5, 10)]

        @testset "No budget constraint" begin
            budget = 0
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            s1 = matching_hungarian_budgeted_lagrangian_search(i, ε)
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            @test s1.value ≈ 5.0 atol = ε # Lagrangian value.

            for s in [s1, s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 5

                for e in max_reward_edges
                    @test e in s.solution
                end

                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 5.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
            end
        end

        @testset "Loose budget constraint" begin
            budget = 1
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            # Don't test Lagrangian relaxation now, constraint is not respected.
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            for s in [s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 5

                n_max_reward_edges =
                    sum(e in s.solution for e in max_reward_edges)
                @test n_max_reward_edges <= 4

                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 3.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 2.0 atol = ε
            end
        end

        @testset "Tight budget constraint" begin
            budget = 5
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            # Don't test Lagrangian relaxation now, constraint is not respected.
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            for s in [s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 5

                for e in max_reward_edges
                    @test !(e in s.solution)
                end

                @test Kombinator._budgeted_bipartite_matching_compute_value(
                    i,
                    s.solution,
                ) ≈ 0.0 atol = ε
                @test Kombinator._budgeted_bipartite_matching_compute_weight(
                    i,
                    s.solution,
                ) ≈ 5.0 atol = ε
            end
        end

        @testset "Impossibly tight budget constraint" begin
            budget = 6
            i = BudgetedBipartiteMatchingInstance(
                graph,
                rewards,
                weights,
                budget,
            )

            # Don't test Lagrangian relaxation now, constraint is not respected.
            s2 = matching_hungarian_budgeted_lagrangian_refinement(i)
            s3 = matching_hungarian_budgeted_lagrangian_approx_half(i)
            s4 = matching_dp_budgeted(i)

            for s in [s2, s3, s4]
                @test s.instance == i
                @test length(s.solution) == 0
            end
        end
    end
end


text:
"""
types.jl
Type definitions for the structured stock-level and agent-level model components
"""

type AgentAssumptions
  """
  Assumptions regarding mortality, movement, and growth
  """
  naturalmortality::Array
  extramortality::Vector
  growth::Vector
  movement::Array
  autonomy::Vector
end

type EnvironmentAssumptions
  """
  A specialized type which contains layers of information to indicate spawning area, habitat type, and additional risks. location id should be specified as NaN when a valid location does not exist.
  """
  spawning::Array
  habitat::Array
  risk::Array
end

type StockDB
  """
  A database which contains population size data for each time step and adult class
  """
  population::DataFrame
  harvest::DataFrame
end

type StockAssumptions
  """
  Age specific survivorship (survivorship at carrying capacity if density depedence occurs)
  Age at 50% mature (Binomial cdf assumed)
  Age specific fecundity (i.e. mean quantity of eggs each spawning female will produce)
  Compensatory fecundity - compensatory strength for changes in fecundity. Compensatory strength is a divisor of K which will result in a 68% change in fecundity - smaller values indicate lower compensation strength. Compensation function based on Normal CDF. Use NaN if compensation is assumed to not occur.
  Compensatory sexual maturity - compensatory strength for changes in age of sexual maturity. Compensatory strength is a divisor of K which will result in a 68% change in age of sexual maturity - smaller values indicate lower compensation strength. Compensation function based on Normal CDF. Use NaN if compensation is assumed to not occur.
  Compensatory mortality - compensatory strength for changes in age of sexual maturity. Compensatory strength is a divisor of K which will result in a 68% change in natural mortality - smaller values indicate lower compensation strength. Compensation function based on Normal CDF. Use NaN if compensation is assumed to not occur.
Age specific catachability
  """
  naturalmortality::Vector
  halfmature::Float64
  broodsize::Vector
  fecunditycompensation::Float64
  maturitycompensation::Float64
  mortalitycompensation::Float64
  catchability::Vector
end


text:
"""
$(SIGNATURES)

Checks if upcoming integration is a new bit for GPSL5.
"""
function is_upcoming_integration_new_bit(gpsl5::GPSL5, prns, num_prns)
    num_prns < 10 && return false
    masked_bit_synchronizer = prns & 0x3ff # First 10
    xored_bit_synchronizer = masked_bit_synchronizer ⊻ 0x35 # 0x35 == 0000110101
    # If xored_bit_synchronizer == 0 -> bit -1 and if xored_bit_synchronizer == 0x3ff -> bit 1
    xored_bit_synchronizer == 0 || xored_bit_synchronizer == 0x3ff
end

function get_default_correlator(gpsl5::GPSL5, num_ants::NumAnts{N}) where N
    EarlyPromptLateCorrelator(num_ants)
end


text:
(::Base.Colon)(start::Epoch{S}, stop::Epoch{S}) where {S} = (:)(start, 1.0seconds, stop)

function (::Base.Colon)(start::Epoch{S}, step::AstroPeriod, stop::Epoch{S}) where S
    step = seconds(step)
    step = start < stop ? step : -step
    StepRangeLen(start, step, floor(Int, value(stop-start)/value(step))+1)
end

Base.step(r::StepRangeLen{T}) where {T<:Epoch} = r.step

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime