text:
"""
Represents an object detection as ``(x_{\\mathrm{midpoint}},
y_{\\mathrm{midpoint}}, \\mathrm{category}``
"""
const Detection2D = Tuple{Float64, Float64, Int64}

"""
Represents an object in the scene as ``(\\mathrm{location}, category)``
"""
const Object3D = Tuple{Float64, Float64, Float64, Int64} #x, y, z, category

export Detection2D
export Object3D

##############################################################################################
#
# struct Multinomial <: Gen.Distribution{Vector{Int}} end
#
# const multinomial = Multinomial()
#
# function Gen.random(::Multinomial, n::Int, probs::AbstractArray{U,1})  where {U <: Real}
#     collect(Int64, map(i -> rand(Distributions.Categorical(probs)), 1:n))
# end
#
# #n is the number of times sampled
# function Gen.logpdf(::Multinomial, x::AbstractArray{Int}, n::Int, probs::AbstractArray{U,1}) where {U <: Real}
#     length(x) != n && return -Inf
#     sum(map(i -> Gen.logpdf(categorical, i, probs), x))
# end
#
# function Gen.logpdf_grad(::Multinomial, x::AbstractArray{Int}, n::Int, probs::AbstractArray{U,1})  where {U <: Real}
# 	gerror("Not implemented")
# 	(nothing, nothing)
# end
#
# (::Multinomial)(n, probs) = Gen.random(Multinomial(), n, probs)
# is_discrete(::Multinomial) = true
#
# has_output_grad(::Multinomial) = false
# has_argument_grads(::Multinomial) = (false,)
#
# export multinomial

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

"Object distribution for 3D detection"
struct Object_Distribution_Present <: Gen.Distribution{Object3D} end

const object_distribution_present = Object_Distribution_Present()

function Gen.random(::Object_Distribution_Present, mu::AbstractVector{V},
                cov::AbstractMatrix{V}, cat::Int64) where {V <: Float64}

    #Detection(rand(Distributions.MvNormal(mu, cov))..., cat)
    rand(Distributions.MvNormal(mu, cov))..., cat
end

function Gen.logpdf(::Object_Distribution_Present, x::Object3D, mu::AbstractVector{V},
                cov::AbstractMatrix{V}, cat::Int64) where {V <: Float64}
    n = length(x)
    #if category mismatch
    cat!=x[n] && return -Inf
    dist = Distributions.MvNormal(mu, cov)
    Distributions.logpdf(dist, collect(x[1:n-1]))
end

function Gen.logpdf_grad(::Object_Distribution_Present, x::AbstractVector{V}, mu::AbstractVector{V},
    cov::AbstractMatrix{V}, cat::Int64) where {V <: Float64}
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution_Present)(mu, cov, cat) = Gen.random(Object_Distribution_Present(), mu, cov, cat)

has_output_grad(::Object_Distribution_Present) = false
has_argument_grads(::Object_Distribution_Present) = (false,)

export object_distribution_present

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

"Object distribution for 2D detection"
struct Object_Distribution_Image <: Gen.Distribution{Detection2D} end

const object_distribution_image = Object_Distribution_Image()

function Gen.random(::Object_Distribution_Image, mu::AbstractVector{V},
                cov::AbstractMatrix{V}, cat::Int64) where {V <: Float64}

    #Detection(rand(Distributions.MvNormal(mu, cov))..., cat)
    rand(Distributions.MvNormal(mu, cov))..., cat
end

function Gen.logpdf(::Object_Distribution_Image, x::Detection2D, mu::AbstractVector{V},
                cov::AbstractMatrix{V}, cat::Int64) where {V <: Float64}
    #if category mismatch
    cat!=x[3] && return -Inf
    dist = Distributions.MvNormal(mu, cov)
    Distributions.logpdf(dist, collect(x[1:2]))
end

function Gen.logpdf_grad(::Object_Distribution_Image, x::AbstractVector{V}, mu::AbstractVector{V},
    cov::AbstractMatrix{V}, cat::Int64) where {V <: Float64}
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution_Image)(mu, cov, cat) = Gen.random(Object_Distribution_Image(), mu, cov, cat)

has_output_grad(::Object_Distribution_Image) = false
has_argument_grads(::Object_Distribution_Image) = (false,)

export object_distribution_image

##############################################################################################
"For a new 3-D object placed anywhere"
struct Object_Distribution <: Gen.Distribution{Object3D} end

const object_distribution = Object_Distribution()

function Gen.random(::Object_Distribution, params::Video_Params)

    c = categorical(params.probs_possible_objects)
    objects_3D = construct_3D(c, params)
end

function Gen.logpdf(::Object_Distribution, object_3D::Object3D, params::Video_Params)
    #categorical
    cat = object_3D[4] #grabbing the category type
    p_categorical = Gen.logpdf(categorical, cat, params.probs_possible_objects)

    #x-coordinate
    #all the x_coordinates
    x = object_3D[1]
    p_x = Gen.logpdf(uniform, x, params.x_min, params.x_max)

    #y-coordinate
    y = object_3D[2]
    p_y = Gen.logpdf(uniform, y, params.y_min, params.y_max)

    #y-coordinate
    z = object_3D[3]
    p_z = Gen.logpdf(uniform, z, params.z_min, params.z_max)

    p_categorical + p_x + p_y + p_z
end

function Gen.logpdf_grad(::Object_Distribution, objects_3D::Object3D, params::Video_Params)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution)(params) = Gen.random(Object_Distribution(), params)

has_output_grad(::Object_Distribution) = false
has_argument_grads(::Object_Distribution) = (false,)

export object_distribution

"Little helper function for constructing 3D objects"
function construct_3D(cat::Int64, params::Video_Params)
    x = uniform(params.x_min, params.x_max)
    y = uniform(params.y_min, params.y_max)
    z = uniform(params.z_min, params.z_max)
    return (x, y, z, cat)
end

##############################################################################################
"For populating the scene with objects spaced out from one another with a sigmoidal prior"
struct Object_Distribution_Logistic <: Gen.Distribution{Array{Object3D}} end

const object_distribution_logistic = Object_Distribution_Logistic()

function Gen.random(::Object_Distribution_Logistic, params::Video_Params)
    num_objects = geometric(params.p_objects)
    placed_objects = Array{Object3D}(undef, 0)
    for i = 1:num_objects
        cat = categorical(params.probs_possible_objects)
        candidate_obj = construct_3D(cat, params) #placed uniformly randomly
        j = 1 #indexes through the objects that have already been placed and accepted
        while j <= length(placed_objects)
            #resample with probabilty given by Logistic with distance
            d = dist(candidate_obj, placed_objects[j])
            p = 1/(1 + exp(-params.k*(d - delta)))
            if bernoulli(p) #if too close, resample
                candidate_obj = construct_3D(cat, params) #placed uniformly randomly
                j = 1 #restart loop over objects already placed
            else
                j = j + 1
            end
        end
        push!(placed_objects, candidate_obj)
    end
    return placed_objects
end

function Gen.logpdf(::Object_Distribution_Logistic, placed_objects::Array{Object3D}, params::Video_Params)
    num_objects = length(placed_objects)
    p_geometric = Gen.logpdf(geometric, num_objects, params.p_objects)

    if num_objects == 0
        return p_geometric
    end

    #iterate over all pairwise combindations of objects and get distances
    p_distance = 0
    for c in combinations(placed_objects, 2)
        d = dist(c[1], c[2])
        p = 1/(1 + exp(-params.k*(d - delta))) #probability of resampling
        inv = 1-p #probability of sticking with it
        log_inv = log(inv)
        p_distance = p_distance + log_inv
    end

    #categorical
    cats = last.(placed_objects) #grabbing the category types
    p_cats = sum(map(i -> Gen.logpdf(categorical, i, params.probs_possible_objects), cats))

    xs = first.(placed_objects)
    p_xs = sum(map(x -> Gen.logpdf(uniform, x, params.x_min, params.x_max), xs))

    ys = (a->a[2]).(placed_objects)
    p_ys = sum(map(y -> Gen.logpdf(uniform, y, params.y_min, params.y_max), ys))

    zs = (a->a[3]).(placed_objects)
    p_zs = sum(map(z -> Gen.logpdf(uniform, z, params.z_min, params.z_max), zs))

    return p_geometric + p_cats + p_xs + p_ys + p_zs + p_distance
end

function Gen.logpdf_grad(::Object_Distribution_Logistic, objects_3D::Object3D, params::Video_Params)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution_Logistic)(params) = Gen.random(Object_Distribution_Logistic(), params)

has_output_grad(::Object_Distribution_Logistic) = false
has_argument_grads(::Object_Distribution_Logistic) = (false,)

export object_distribution_logistic

##############################################################################################
"For populating the scene with objects spaced out from one another with a graded, Gaussian prior"
struct Object_Distribution_Gaussian <: Gen.Distribution{Array{Object3D}} end

const object_distribution_gaussian = Object_Distribution_Gaussian()

function Gen.random(::Object_Distribution_Gaussian, params::Video_Params)
    num_objects = geometric(params.p_objects)
    placed_objects = Array{Object3D}(undef, 0)
    for i = 1:num_objects
        cat = categorical(params.probs_possible_objects)
        candidate_obj = construct_3D(cat, params) #placed uniformly randomly
        j = 1 #indexes through the objects that have already been placed and accepted
        while j <= length(placed_objects)
            #resample with probabilty given by Gaussian with distance
            d = dist(candidate_obj, placed_objects[j])
            p = exp(-0.5*(d/params.delta)^2) #one when dist==0
            if bernoulli(p) #if too close, resample
                candidate_obj = construct_3D(cat, params) #placed uniformly randomly
                j = 1 #restart loop over objects already placed
            else
                j = j + 1
            end
        end
        push!(placed_objects, candidate_obj)
    end
    return placed_objects
end

function Gen.logpdf(::Object_Distribution_Gaussian, placed_objects::Array{Object3D}, params::Video_Params)
    num_objects = length(placed_objects)
    p_geometric = Gen.logpdf(geometric, num_objects, params.p_objects)

    if num_objects == 0
        return p_geometric
    end

    #iterate over all pairwise combindations of objects and get distances
    p_distance = 0
    for c in combinations(placed_objects, 2)
        d = dist(c[1], c[2])
        p = exp(-0.5*(d/params.delta)^2) #probability of resampling
        inv = 1-p #probability of sticking with it
        log_inv = log(inv)
        p_distance = p_distance + log_inv
    end

    #categorical
    cats = last.(placed_objects) #grabbing the category types
    p_cats = sum(map(i -> Gen.logpdf(categorical, i, params.probs_possible_objects), cats))

    xs = first.(placed_objects)
    p_xs = sum(map(x -> Gen.logpdf(uniform, x, params.x_min, params.x_max), xs))

    ys = (a->a[2]).(placed_objects)
    p_ys = sum(map(y -> Gen.logpdf(uniform, y, params.y_min, params.y_max), ys))

    zs = (a->a[3]).(placed_objects)
    p_zs = sum(map(z -> Gen.logpdf(uniform, z, params.z_min, params.z_max), zs))

    return p_geometric + p_cats + p_xs + p_ys + p_zs + p_distance
end

function Gen.logpdf_grad(::Object_Distribution_Gaussian, objects_3D::Object3D, params::Video_Params)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution_Gaussian)(params) = Gen.random(Object_Distribution_Gaussian(), params)

has_output_grad(::Object_Distribution_Gaussian) = false
has_argument_grads(::Object_Distribution_Gaussian) = (false,)

export object_distribution_gaussian
##############################################################################################
"For populating the scene with objects spaced out from one another. no two object can be within distance delta of each other"
struct Object_Distribution_Delta <: Gen.Distribution{Array{Object3D}} end

const object_distribution_delta = Object_Distribution_Delta()

function Gen.random(::Object_Distribution_Delta, params::Video_Params)
    num_objects = geometric(params.p_objects)
    placed_objects = Array{Object3D}(undef, 0)
    for i = 1:num_objects
        cat = categorical(params.probs_possible_objects)
        candidate_obj = construct_3D(cat, params) #placed uniformly randomly
        j = 1 #indexes through the objects that have already been placed and accepted
        while j <= length(placed_objects)
            if dist(candidate_obj, placed_objects[j]) < params.delta #if too close, resample
                candidate_obj = construct_3D(cat, params) #placed uniformly randomly
                j = 1 #restart loop over objects already placed
            else
                j = j + 1
            end
        end
        push!(placed_objects, candidate_obj)
    end
    return placed_objects
end

function Gen.logpdf(::Object_Distribution_Delta, placed_objects::Array{Object3D}, params::Video_Params)
    #either the objects are spaced out enough or they aren't.
    #iterate over all pairwise combindations of objects
    for c in combinations(placed_objects, 2)
        if dist(c[1], c[2]) < params.delta
            #println("c ", c)
            #println("c[1] ", c[1])
            #println("c[2] ", c[2])
            return -Inf #if the objects aren't far enough apart, log probability is -Inf
        end
    end

    num_objects = length(placed_objects)
    p_geometric = Gen.logpdf(geometric, num_objects, params.p_objects)

    if num_objects == 0
        return p_geometric
    end

    #categorical
    cats = last.(placed_objects) #grabbing the category types
    p_cats = sum(map(i -> Gen.logpdf(categorical, i, params.probs_possible_objects), cats))

    xs = first.(placed_objects)
    p_xs = sum(map(x -> Gen.logpdf(uniform, x, params.x_min, params.x_max), xs))

    ys = (a->a[2]).(placed_objects)
    p_ys = sum(map(y -> Gen.logpdf(uniform, y, params.y_min, params.y_max), ys))

    zs = (a->a[3]).(placed_objects)
    p_zs = sum(map(z -> Gen.logpdf(uniform, z, params.z_min, params.z_max), zs))

    return p_geometric + p_cats + p_xs + p_ys + p_zs
end

function Gen.logpdf_grad(::Object_Distribution_Delta, objects_3D::Object3D, params::Video_Params)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution_Delta)(params) = Gen.random(Object_Distribution_Delta(), params)

has_output_grad(::Object_Distribution_Delta) = false
has_argument_grads(::Object_Distribution_Delta) = (false,)

export object_distribution_delta

#tiny helper function for getting distance between two objects
function dist(a::Object3D, b::Object3D)
    sqrt((a[1]-b[1])^2 + (a[2]-b[2])^2 + (a[3]-b[3])^2)
end

##############################################################################################
"For a new 3-D object placed anywhere. category already determined"
struct Location_Distribution_Uniform <: Gen.Distribution{Object3D} end

const location_distribution_uniform = Location_Distribution_Uniform()

function Gen.random(::Location_Distribution_Uniform, cat::Int64, params::Video_Params)
    objects_3D = construct_3D(cat, params)
end

function Gen.logpdf(::Location_Distribution_Uniform, object_3D::Object3D, cat::Int64, params::Video_Params)
    #could add something making sure category is right.

    #x-coordinate
    #all the x_coordinates
    x = object_3D[1]
    p_x = Gen.logpdf(uniform, x, params.x_min, params.x_max)

    #y-coordinate
    y = object_3D[2]
    p_y = Gen.logpdf(uniform, y, params.y_min, params.y_max)

    #y-coordinate
    z = object_3D[3]
    p_z = Gen.logpdf(uniform, z, params.z_min, params.z_max)

    p_x + p_y + p_z
end

function Gen.logpdf_grad(::Location_Distribution_Uniform, objects_3D::Object3D, cat::Int64, params::Video_Params)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Location_Distribution_Uniform)(cat, params) = Gen.random(Location_Distribution_Uniform(), cat, params)

has_output_grad(::Location_Distribution_Uniform) = false
has_argument_grads(::Location_Distribution_Uniform) = (false,)

export location_distribution_uniform

##############################################################################################
"For a new 3-D object placed along certain line segments with Gaussian noise, category has already been chosen"
struct New_Location_Distribution_Noisy <: Gen.Distribution{Object3D} end

const new_location_distribution_noisy = New_Location_Distribution_Noisy()

function Gen.random(::New_Location_Distribution_Noisy, cat::Int64, params::Video_Params, line_segments::Array{Line_Segment})
    n = size(line_segments)[1]
    if n < 0
        println("problem. sampling location from new_location_distribution_noisy when we shouldn't be")
    end
    i = categorical(fill(1/n, n)) #which line segment
    line_segment = line_segments[i]
    #println("line_segment ", line_segment)
    #length of the line segment
    length = sqrt((line_segment.start.x-line_segment.endpoint.x)^2 + (line_segment.start.y-line_segment.endpoint.y)^2 + (line_segment.start.z-line_segment.endpoint.z)^2)
    #println("length ", length)
    d = uniform(0, length) #sample a distance along the line segment
    #println("d ", d)
    #plug in d/length as t
    x = line_segment.start.x + line_segment.a*(d/length)
    y = line_segment.start.y + line_segment.b*(d/length)
    z = line_segment.start.z + line_segment.c*(d/length)

    #add noise
    dist = trunc_normal(0., 0.1, 0., 10.) #mean noise is 0, sd is 0.1?
    #println("dist ", dist)
    angle = uniform(0, 2*pi)

    (x_noisy, y_noisy, z_noisy) = sample_point(line_segment, Coordinate(x,y,z), dist, angle)

    return (x_noisy, y_noisy, z_noisy, cat)
end

function Gen.logpdf(::New_Location_Distribution_Noisy, object_3D::Object3D, cat::Int64, params::Video_Params, line_segments::Array{Line_Segment})
    #check how many line segments have this point on it
    n = length(line_segments)
    point = Coordinate(object_3D[1], object_3D[2], object_3D[3])
    ps = map(prob_given_line_segment, fill(point, n), line_segments)
    p_point = log(1/n) + log(sum(ps)) #could do something with categorical distibution, but don't know how
    return p_point
end

function Gen.logpdf_grad(::New_Location_Distribution_Noisy, objects_3D::Object3D, cat::Int64, params::Video_Params, line_segments::Array{Line_Segment})
    gerror("Not implemented")
    (nothing, nothing)
end

(::New_Location_Distribution_Noisy)(cat, params, line_segments) = Gen.random(New_Location_Distribution_Noisy(), cat, params, line_segments)

has_output_grad(::New_Location_Distribution_Noisy) = false
has_argument_grads(::New_Location_Distribution_Noisy) = (false,)

function sample_point(line_segment::Line_Segment, point_on_line::Coordinate, dist::Float64, angle::Float64)

    #first, sample a vector on the plane. the plane is defined by being normal to the line segment and including the point
    #point_on_line. describe this plane in point normal form. pick any two values of x and y and then solve for z to get
    #another point on this plane. then a vector on this plane is given by (x-point_on_line_x, y-point_on_line_y, z-point_on_line_z)
    x = 1
    y = 1
    z = (-line_segment.a*(x-point_on_line.x) - line_segment.b*(y-point_on_line.x))/line_segment.c + point_on_line.y
    v1 = [x-point_on_line.x, y-point_on_line.y, z-point_on_line.z]

    #get v2, which is perpendicular to the line segement and to v1, bu taking the cross product.
    v2 = cross(v1, [line_segment.a, line_segment.b, line_segment.c])

    #make unit vectors
    v1 = v1./sqrt(sum(v1.^2))
    v2 = v2./sqrt(sum(v2.^2))

    point = [point_on_line.x, point_on_line.y, point_on_line.z] + dist * sin(angle) .*v1 + dist * cos(angle) .*v2
    return (point[1], point[2], point[3])
end

#gives the probability of a point given a segment
function prob_given_line_segment(point::Coordinate, line_segment::Line_Segment)

    #point is on a plane normal to line segment. that's given by point normal form. finding intersection between
    #this plane and the line segment by solving for t
    t = (line_segment.a*(point.x - line_segment.start.x) + line_segment.b*(point.y - line_segment.start.y) + line_segment.c*(point.z - line_segment.start.z))/(line_segment.a^2 + line_segment.b^2 + line_segment.c^2)
    if t < 0 || t > 1 #would mean that this is off the line segment
        return 0
    end

    line_seg_coef = [line_segment.a, line_segment.b, line_segment.c]
    line_seg_start = [line_segment.start.x, line_segment.start.y, line_segment.start.z]
    point_on_line = t.*line_seg_coef .+ line_seg_start

    #get distance between point_on_line and point
    point = [point.x, point.y, point.z]
    dist = sqrt(sum((point - point_on_line).^2)) #norm might do the same

    #actually calculate probabilities
    len = sqrt((line_segment.start.x-line_segment.endpoint.x)^2 + (line_segment.start.y-line_segment.endpoint.y)^2 + (line_segment.start.z-line_segment.endpoint.z)^2)
    p1 = Gen.logpdf(uniform, 0, 0, len)#doesn't matter what the actual value is because it's uniform, so plug in 0
    p2 = Gen.logpdf(trunc_normal, dist, 0., 0.1, 0., 10.)#make sure it matches above
    p3 = Gen.logpdf(uniform, 0, 0, 2*pi)

    return MathConstants.e^(p1 + p2 + p3) #transform it back to probability space
end

export new_location_distribution_noisy

##############################################################################################
"For changing the location of an object: either by Gaussian noise based on previous location, or resampled along certain line segments with Gaussian noise, category has already been chosen"
struct New_Location_Distribution_Noisy_Or_Gaussian <: Gen.Distribution{Object3D} end

const new_location_distribution_noisy_or_gaussian = New_Location_Distribution_Noisy_Or_Gaussian()

function Gen.random(::New_Location_Distribution_Noisy_Or_Gaussian, mu::Vector{Float64},
    cov::Matrix{Float64}, cat::Int64, params::Video_Params, line_segments::Array{Array{Line_Segment,1},1})

    if length(line_segments[cat]) > 0
        #coin flip.
        if bernoulli(0.5)
            #println("from data-driven distribution")
            to_return = new_location_distribution_noisy(cat, params, line_segments[cat])
        else
            #println("from gaussian around previous location")
            to_return = object_distribution_present(mu::AbstractVector{Float64}, cov::AbstractMatrix{Float64}, cat::Int64)
        end
    else
        to_return = object_distribution_present(mu::AbstractVector{Float64}, cov::AbstractMatrix{Float64}, cat::Int64)
    end

    return to_return
end

function Gen.logpdf(::New_Location_Distribution_Noisy_Or_Gaussian, x::Object3D, mu::AbstractVector{Float64},
    cov::AbstractMatrix{Float64}, cat::Int64, params::Video_Params, line_segments::Array{Array{Line_Segment,1},1})

    p_from_gaussian = Gen.logpdf(object_distribution_present, x, mu, cov, cat)

    if length(line_segments) < 0
        return p_from_gaussian
    else
        p_from_noisy = Gen.logpdf(new_location_distribution_noisy, x, cat, params, line_segments[cat])
        return log(0.5*(MathConstants.e^p_from_noisy + MathConstants.e^p_from_gaussian))
    end
end

function Gen.logpdf_grad(::New_Location_Distribution_Noisy_Or_Gaussian, mu::AbstractVector{Float64},
    cov::AbstractMatrix{Float64}, cat::Int64, params::Video_Params, line_segments::Array{Line_Segment})
    gerror("Not implemented")
    (nothing, nothing)
end

(::New_Location_Distribution_Noisy_Or_Gaussian)(mu, cov, cat, params, line_segments) = Gen.random(New_Location_Distribution_Noisy_or_Gaussian(), mu, cov, cat, params, line_segments)

has_output_grad(::New_Location_Distribution_Noisy_Or_Gaussian) = false
has_argument_grads(::New_Location_Distribution_Noisy_Or_Gaussian) = (false,)

export new_location_distribution_noisy_or_gaussian

##############################################################################################
#Sample the category. If the object category has never been observed, sample location from a uniform distribution
#If it has been observed, sample location with 10% chance from a uniform or from new_object_distribution_noisy
struct New_Object_Distribution_Noisy_Or_Uniform <: Gen.Distribution{Object3D} end

const new_object_distribution_noisy_or_uniform = New_Object_Distribution_Noisy_Or_Uniform()

function Gen.random(::New_Object_Distribution_Noisy_Or_Uniform, params::Video_Params, line_segments_per_category::Array{Array{Line_Segment,1},1})
    cat = categorical(params.probs_possible_objects)
    line_segments = line_segments_per_category[cat]

    if length(line_segments) > 0
        coin_flip = bernoulli(0.9)
        if coin_flip
            #println("from data-driven distribution")
            to_return = new_location_distribution_noisy(cat, params, line_segments)
        else
            #println("from uniform distribution")
            to_return = location_distribution_uniform(cat, params)
        end
    else
        #println("from uniform distribution")
        to_return = location_distribution_uniform(cat, params)
    end
    return to_return
end

function Gen.logpdf(::New_Object_Distribution_Noisy_Or_Uniform, object_3D::Object3D, params::Video_Params, line_segments_per_category::Array{Array{Line_Segment,1},1})
    cat = object_3D[4]
    line_segments = line_segments_per_category[cat]
    p_cat = Gen.logpdf(categorical, cat, params.probs_possible_objects)

    if length(line_segments) > 0
        p_from_noisy = Gen.logpdf(new_location_distribution_noisy, object_3D, cat, params, line_segments)
        p_from_uniform = Gen.logpdf(location_distribution_uniform, object_3D, cat, params)
        return p_cat + log(0.9*(MathConstants.e^p_from_noisy + MathConstants.e^p_from_uniform)) #0.5 based on p in coin_flip bernoulli
    else #no chance from new_location_distribution_noisy
        p_from_uniform = Gen.logpdf(location_distribution_uniform, object_3D, cat, params)
        return p_cat + p_from_uniform
    end
end

function Gen.logpdf_grad(::New_Object_Distribution_Noisy_Or_Uniform, objects_3D::Object3D, params::Video_Params, line_segments_per_category::Array{Array{Line_Segment,1},1})
    gerror("Not implemented")
    (nothing, nothing)
end

(::New_Object_Distribution_Noisy_Or_Uniform)(params, line_segments_per_category) = Gen.random(New_Object_Distribution_Noisy_Or_Uniform(), params, line_segments_per_category)

has_output_grad(::New_Object_Distribution_Noisy_Or_Uniform) = false
has_argument_grads(::New_Object_Distribution_Noisy_Or_Uniform) = (false,)

export new_object_distribution_noisy_or_uniform

##############################################################################################
#For a new 3-D object with the same location but new category
struct Object_Distribution_Category <: Gen.Distribution{Object3D} end

const object_distribution_category = Object_Distribution_Category()

function Gen.random(::Object_Distribution_Category, x::Float64, y::Float64, z::Float64, perturb_params::Perturb_Params)

    cat = categorical(perturb_params.probs_possible_objects)

    return (x, y, z, cat)
end

function Gen.logpdf(::Object_Distribution_Category, object_3D::Object3D, x::Float64, y::Float64, z::Float64, perturb_params::Perturb_Params)
    #categorical
    cat = object_3D[4] #grabbing the category type
    Gen.logpdf(categorical, cat, perturb_params.probs_possible_objects)
end

function Gen.logpdf_grad(::Object_Distribution_Category, objects_3D::Object3D, perturb_params::Perturb_Params)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Object_Distribution_Category)(x, y, z, perturb_params) = Gen.random(Object_Distribution_Category(), x, y, z, perturb_params)

has_output_grad(::Object_Distribution_Category) = false
has_argument_grads(::Object_Distribution_Category) = (false,)

export new_object_distribution_category

##############################################################################################
struct Hallucination_Distribution <: Gen.Distribution{Detection2D} end

const hallucination_distribution = Hallucination_Distribution()

function Gen.random(::Hallucination_Distribution, params::Video_Params, v::Matrix{Real}, rec_field::Receptive_Field)

    #saying there must be no fewer than 0 objects per scene, and at most 100
    #numObjects = @trace(trunc_poisson(sum(params.v[:,1]), -1.0, 100.0), (:numObjects)) #may want to truncate so 0 objects isn't possible
    #objects = @trace(multinomial_objects(numObjects,[0.2,0.2,0.2,0.2,0.2], ), (:objects))
    ws = @inbounds 1.0/sum(v[:,1]) * v[:,1]
    c = categorical(ws)
    detection_2D = construct_2D(c, params, rec_field)
end

function Gen.logpdf(::Hallucination_Distribution, detection_2D::Detection2D, params::Video_Params, v::Matrix{Real}, rec_field::Receptive_Field)

    #multinomial
    cat = detection_2D[3] #grabbing the category type
    ws = @inbounds 1.0/sum(v[:,1]) * v[:,1]
    p_categorical = Gen.logpdf(categorical, cat, ws) #prob_vec is the hallucination lambdas normalized

    #x-coordinate
    #all the x_coordinates
    x = detection_2D[1]
    p_x = Gen.logpdf(uniform, x, rec_field.p1[1], rec_field.p2[1])

    #y-coordinate
    y = detection_2D[2]
    p_y = Gen.logpdf(uniform, y, rec_field.p1[2], rec_field.p2[2])

    p_categorical + p_x + p_y
end

function Gen.logpdf_grad(::Hallucination_Distribution, detection_2D::Detection2D, params::Video_Params, v::Matrix{Real}, rec_field::Receptive_Field)
    gerror("Not implemented")
    (nothing, nothing)
end

(::Hallucination_Distribution)(params, v, rec_field) = Gen.random(Hallucination_Distribution(), params, v, rec_field)

has_output_grad(::Hallucination_Distribution) = false
has_argument_grads(::Hallucination_Distribution) = (false,)

export hallucination_distribution

#little helper function for constructing 2D objects
function construct_2D(cat::Int64, params::Video_Params, rec_field::Receptive_Field)
    x = uniform(rec_field.p1[1], rec_field.p2[1])
    y = uniform(rec_field.p1[2], rec_field.p2[2])
    return (x, y, cat)
end

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

#Distributions
#TruncatedPoisson
export trunc_poisson

struct TruncatedPoisson <: Gen.Distribution{Int} end

const trunc_poisson = TruncatedPoisson()

function Gen.logpdf(::TruncatedPoisson, x::Int, lambda::U, low::U, high::U) where {U <: Real}
	d = Distributions.Poisson(lambda)
	td = Distributions.Truncated(d, low, high)
	Distributions.logpdf(td, x)
end

function Gen.logpdf_grad(::TruncatedPoisson, x::Int, lambda::U, low::U, high::U)  where {U <: Real}
	gerror("Not implemented")
	(nothing, nothing)
end

function Gen.random(::TruncatedPoisson, lambda::U, low::U, high::U)  where {U <: Real}
	d = Distributions.Poisson(lambda)
	rand(Distributions.Truncated(d, low, high))
end

(::TruncatedPoisson)(lambda, low, high) = random(TruncatedPoisson(), lambda, low, high)
is_discrete(::TruncatedPoisson) = true

has_output_grad(::TruncatedPoisson) = false
has_argument_grads(::TruncatedPoisson) = (false,)

export trunc_poisson

##############################################################################################
#TruncatedNormal
export trunc_normal

#small issue where all of the inputs need to by Float64. Doesn't accept Int64s
struct TruncatedNormal <: Gen.Distribution{Float64} end

const trunc_normal = TruncatedNormal()

function Gen.random(::TruncatedNormal, mu::U, std::U, low::U, high::U)  where {U <: Real}
	n = Distributions.Normal(mu, std)
	rand(Distributions.Truncated(n, low, high))
end

function Gen.logpdf(::TruncatedNormal, x::U, mu::U, std::U, low::U, high::U) where {U <: Real}
	n = Distributions.Normal(mu, std)
	tn = Distributions.Truncated(n, low, high)
	Distributions.logpdf(tn, x)
end

function Gen.logpdf_grad(::TruncatedNormal, x::U, mu::U, std::U, low::U, high::U)  where {U <: Real}
    precision = 1. / (std * std)
    diff = mu - x
    deriv_x = diff * precision
    deriv_mu = -deriv_x
    deriv_std = -1. / std + (diff * diff) / (std * std * std)

    if x<=low
        deriv_x = log(0.001) #trying to have a very small positive gradient
    elseif x>=high
        deriv_x = log(-0.001) #trying to have a very small negative gradient
    end

    (deriv_x, deriv_mu, deriv_std)
end

(::TruncatedNormal)(mu, std, low, high) = random(TruncatedNormal(), mu, std, low, high)
is_discrete(::TruncatedNormal) = false
has_output_grad(::TruncatedNormal) = true
has_argument_grads(::TruncatedNormal) = (true, true, false, false) #just has output gradients for mu and std

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
# visibility mask = [
#     1 1 1;  # x_11 x_12 x_13
#     1 0 1;  # x_21      x_23
#     0 1 1;  #      x_32 x_33
#     1 1 0;  # x_41 x_42
# ]

viewpoint_indices = [1, 2, 3, 1, 3, 2, 3, 1, 2]
point_indices = [1, 1, 1, 2, 2, 3, 3, 4, 4]
indices = SBA.Indices(viewpoint_indices, point_indices)

@test SBA.n_points(indices) == 4
@test SBA.n_viewpoints(indices) == 3

# get array indices of X for x_*1
# where x_*1 are projections of all visible points in the 1st viewpoint
@test SBA.points_by_viewpoint(indices, 1) == [1, 4, 8]

# get array indices of X for x_*3
@test SBA.points_by_viewpoint(indices, 3) == [3, 5, 7]

# get array indices of X for x_1*
# where x_1* are projections of the 1st 3D point in all observable viewpoints
@test SBA.viewpoints_by_point(indices, 1) == [1, 2, 3]

# get array indices of X for x_4*
@test SBA.viewpoints_by_point(indices, 4) == [8, 9]

# [x_12, x_32] and [x_13, x_33] are shared
@test SBA.shared_point_indices(indices, 2, 3) == ([2, 6], [3, 7])

# [x_11, x_21] and [x_13, x_23] are shared
@test SBA.shared_point_indices(indices, 1, 3) == ([1, 4], [3, 5])

# another case
# n_points = 2, n_viewpoints = 3,
# indices     1    2    3
#   X =   [x_13 x_21 x_22]

# then the corresponding mask should be
mask = [
    0 0 1;  #           x_13
    1 1 0;  # x_21 x_22
]

viewpoint_indices = [3, 1, 2]
point_indices = [1, 2, 2]
indices = SBA.Indices(viewpoint_indices, point_indices)

@test SBA.n_points(indices) == 2
@test SBA.n_viewpoints(indices) == 3

# get array indices of X for x_*1
@test SBA.points_by_viewpoint(indices, 1) == [2]

# get array indices of X for x_*3
@test SBA.points_by_viewpoint(indices, 3) == [1]

# get array indices of X for x_1*
@test SBA.viewpoints_by_point(indices, 1) == [1]

# get array indices of X for x_2*
@test SBA.viewpoints_by_point(indices, 2) == [2, 3]

# [x_21] and [x_22] are shared
@test SBA.shared_point_indices(indices, 1, 2) == ([2], [3])

# no points are shared
@test isnothing(SBA.shared_point_indices(indices, 1, 3))

# second row has only zero elements
# visibility mask = [
#     1 0 1 0;
#     0 0 0 0;
#     0 1 1 1
# ]
viewpoint_indices = [1, 3, 2, 3, 4]
point_indices = [1, 1, 3, 3, 3]
@test_throws AssertionError SBA.Indices(viewpoint_indices, point_indices)

# third column has only zero elements
# visibility mask = [
#     1 0 0 1;
#     0 1 0 0;
#     0 1 0 1
# ]
viewpoint_indices = [1, 4, 2, 2, 4]
point_indices = [1, 1, 2, 3, 3]
@test_throws AssertionError SBA.Indices(viewpoint_indices, point_indices)

ext: jl


text:
# File: soem.jl
# Project: example
# Created Date: 30/12/2020
# Author: Shun Suzuki
# -----
# Last Modified: 17/12/2021
# Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp)
# -----
# Copyright (c) 2020 Hapis Lab. All rights reserved.
#

using Printf

using AUTD3
using StaticArrays
using StaticArrays: size

function simple(autd::AUTD)
    set_silent_mode(autd, true)

    g = focal_point_gain(SVector(90.0, 80.0, 150.0))
    freq::Int32 = 150
    m = sine_modulation(freq)

    send(autd, m, g)

    println("press enter to exit...")
    readline()

    dispose(g)
    dispose(m)
    stop(autd)
end

function bessel(autd::AUTD)
    set_silent_mode(autd, true)

    g = bessel_beam_gain(SVector(90.0, 80.0, 150.0), SVector(0.0, 0.0, 1.0), 13.0 / 180.0 * pi)
    freq::Int32 = 150
    m = sine_modulation(freq)

    send(autd, m, g)

    println("press enter to exit...")
    readline()

    dispose(g)
    dispose(m)
    stop(autd)
end

function holo(autd::AUTD)
    set_silent_mode(autd, true)

    foci = [
        SVector(120.0, 80.0, 150.0),
        SVector(60.0, 80.0, 150.0)
    ]
    amps = [1.0, 1.0]
    g = holo_gain_sdp(foci, amps)

    freq::Int32 = 150
    m = sine_modulation(freq)

    send(autd, m, g)

    println("press enter to exit...")
    readline()

    dispose(g)
    dispose(m)
    stop(autd)
end

function stm_test(autd::AUTD)
    set_silent_mode(autd, false)

    m = static_modulation()
    send(autd, m, Nothing)

    center = SVector(90.0, 80.0, 150.0)

    radius = 30.0
    size = 200
    stm_cnt = stm(autd)
    for i = 1:size
        theta::Float64 = 2pi * i / size
        r = center + radius * SVector(cos(theta), sin(theta), 0)
        f = focal_point_gain(r)
        add_stm_gain(stm_cnt, f)
        dispose(f)
    end

    freq::Float64 = 1.0
    start_stm(stm_cnt, freq)

    println("press enter to exit...")
    readline()

    stop_stm(stm_cnt)
    finish_stm(stm_cnt)

    dispose(m)
    stop(autd)
end

function seq(autd::AUTD)
    set_silent_mode(autd, false)

    m = static_modulation()

    seq = sequence()

    center = SVector(90.0, 80.0, 150.0)
    radius = 30.0
    size::UInt64 = 200
    for i = 1:size
        theta = 2.0 * pi * i / size
        p = radius * SVector(cos(theta), sin(theta), 0.0)
        add_point(seq, center + p)
    end
    freq::Float64 = 200.0
    set_freq(seq, freq)

    send(autd, m, seq)

    println("press enter to exit...")
    readline()

    dispose(seq)
    dispose(m)
    stop(autd)
end

function seq_gain(autd::AUTD)
    set_silent_mode(autd, false)

    m = static_modulation()

    center = SVector(90.0, 80.0, 150.0)
    radius = 30.0
    size::UInt64 = 200
    seq = gain_sequence(autd, duty_phase_full)

    for i = 1:size
        theta = 2.0 * pi * i / size
        p = radius * SVector(cos(theta), sin(theta), 0.0)
        g = focal_point_gain(center + p)
        add_gain(seq, g)
    end

    freq::Float64 = 1
    set_freq(seq, freq)

    send(autd, m, seq)

    println("press enter to exit...")
    readline()

    dispose(seq)
    dispose(m)
    stop(autd)
end

function run(autd::AUTD)
    samples = [
        (simple, "Single Focal Point Sample"),
        (bessel, "Bessel beam Sample"),
        (holo, "Multiple Focal Points Sample"),
        (stm_test, "Spatio-Temporal Modulation Sample"),
        (seq, "PointSequence (Hardware STM) Sample"),
        (seq_gain, "GainSequence (Hardware STM with arbitrary Gain) Sample")
    ]

    firm_info_list = firmware_info_list(autd)
    for (i, firm_info) in enumerate(firm_info_list)
        @printf("AUTD[%d]: CPU: %s, FPGA: %s\n", i, firm_info[1], firm_info[2])
    end

    clear(autd)

    while true
        for (i, (_, name)) in enumerate(samples)
            @printf("[%d]: %s\n", i, name)
        end
        println("[Other]: finish")
        print("Choose number: ")

        idx = tryparse(Int64, readline())
        if idx === nothing || idx > length(samples) || idx < 1
            break
        end

        (fn, _) = samples[idx]
        fn(autd)
    end

    clear(autd)
    dispose(autd)
end

function get_adapter()
    adapters = enumerate_adapters()
    for (i, adapter) in enumerate(adapters)
        @printf("[%d]: %s, %s\n", i, adapter[2], adapter[1])
    end

    print("Input number: ")
    idx = tryparse(Int64, readline())
    if idx === nothing || idx > length(adapters) || idx < 1
        println("choose correct number!")
        return ""
    end

    adapters[idx][2]
end

function main()
    autd = AUTD()

    add_device(autd, SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0))

    adapter = get_adapter()
    link = soem_link(adapter, num_devices(autd))
    if !open_autd(autd, link)
        println(last_error())
        return
    end

    run(autd)
end

main()


text:
#
# Dimension: translates sequences of Symbols, Strings, or Ints to
# ordinal values, so the Dimension keys can be used for lookup in
# Arrays in Mimi. A Dimension can be declared one of several ways:
# dim = Dimension(:foo, :bar, :baz)      # varargs
# dim = Dimension([:foo, :bar, :baz])    # Vector
# dim = Dimension(2010:2100)             # Range
# dim = Dimension(4)                     # Same as 1:4
#
# Similarly, values can be referenced several ways:
# value  = Dimension[:foo]               # single lookup
# values = Dimension[(:foo, :bar)]       # Tuple of values
# values = Dimension[[2010, 2020, 2030]] # Vector of values
# values = Dimension[2010, 2020, 2030]   # Varargs
# values = Dimension[2010:10:2030]       # AbstractRange
#
using DataStructures

key_type(dim::Dimension{T}) where {T <: DimensionKeyTypes} = T

#
# Iteration and basic dictionary methods are delegated to the internal dict
#
Base.length(dim::Dimension) = length(dim.dict)
Base.iterate(dim::Dimension, state...) = iterate(dim.dict, state...)

Base.keys(dim::Dimension)   = keys(dim.dict)
Base.values(dim::Dimension) = values(dim.dict)

Base.get(dim::Dimension, key::Union{Number, Symbol, String}, default::Any) = get(dim.dict, key, default)

Base.getindex(dim::Dimension, key::Colon) = collect(values(dim.dict))

Base.getindex(dim::AbstractDimension, key::Union{Number, Symbol, String}) = getindex(dim.dict, key)

# Support dim[[:foo, :bar, :baz]], dim[(:foo, :bar, :baz)], and dim[2010:2020]
Base.getindex(dim::AbstractDimension, keys::Union{Vector{T} where T, Tuple, AbstractRange}) = [getindex(dim.dict, key) for key in keys]

Base.getindex(dim::AbstractDimension, keys...) = getindex(dim, keys)

Base.length(dim::RangeDimension) = length(dim.range)
Base.iterate(dim::RangeDimension, state...)  = iterate(dim.range, state...)

Base.keys(dim::RangeDimension)   = collect(dim.range)
Base.values(dim::RangeDimension) = collect(1:length(dim.range))

# Get last value from OrderedDict of keys
Base.lastindex(dim::AbstractDimension) = dim.dict.keys[length(dim)]

#
# Compute the index of a "key" (e.g., a year) in the range.
#
function Base.get(dim::RangeDimension, key::Int, default::Any=0)
    r = dim.range
    i = key - r.start
    return i == 0 ? 1 : (i % r.step != 0 ? default : 1 + div(i, r.step))
end

# Support dim[[2010, 2020, 2030]], dim[(:foo, :bar, :baz)], and dim[2010:2050]
Base.getindex(dim::RangeDimension, keys::Union{Vector{Int}, Tuple, AbstractRange}) = [get(dim, key, 0) for key in keys]


text:
# module for object estimated by slam

using Plots
pyplot()

include(joinpath(split(@__FILE__, "src")[1], "src/common/covariance_ellipse/covariance_ellipse.jl"))

mutable struct EstimatedObject
  pose
  id
  shape
  color
  size
  cov

  # init
  function EstimatedObject(x, y; shape=:star, color=:blue,
                           size=15, id=-1, cov=nothing)
    self = new()
    self.pose = [x, y]
    self.id = id
    self.shape = shape
    self.color = color
    self.size = size
    self.cov = cov
    return self
  end
end

function draw!(self::EstimatedObject)
  if !(self.cov === nothing)
    scatter!([self.pose[1]], [self.pose[2]], markershape=self.shape,
             markercolor=self.color, markersize=self.size)
    if !(self.id === -1)
      annotate!(self.pose[1], self.pose[2],
                text("id:$(self.id)", :black, :left, 10))
    end
    draw_covariance_ellipse!(self.pose, self.cov, 3)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# ================================================================
# Reach set
# ================================================================

"""
    ReachSet{N, ST<:LazySet{N}} <: AbstractLazyReachSet{N}

Type that wraps a reach-set using a `LazySet` as underlying representation.

### Fields

- `X`  -- set
- `Δt` -- time interval

### Notes

A `ReachSet` is a struct representing (an approximation of) the reachable states
for a given time interval. The type of the representation is `ST`, which may be
any subtype LazySet. For efficiency reasons, `ST` should be concretely typed.

By assumption the coordinates in this reach-set are associated to the integers
`1, …, n`. The function `vars` returns such tuple.
"""
struct ReachSet{N, ST<:LazySet{N}} <: AbstractLazyReachSet{N}
    X::ST
    Δt::TimeInterval
end

# abstract reach set interface functions
set(R::ReachSet) = R.X
setrep(R::ReachSet{N, ST}) where {N, ST<:LazySet{N}} = ST
setrep(::Type{ReachSet{N, ST}}) where {N, ST<:LazySet{N}} = ST
setrep(::AbstractVector{RT}) where {RT<:AbstractReachSet} = setrep(RT)
tstart(R::ReachSet) = inf(R.Δt)
tend(R::ReachSet) = sup(R.Δt)
tspan(R::ReachSet) = R.Δt
dim(R::ReachSet) = dim(R.X)
vars(R::ReachSet) = Tuple(Base.OneTo(dim(R.X)),)

# no-op
function Base.convert(T::Type{ReachSet{N, ST}}, R::ReachSet{N, ST}) where {N, ST<:LazySet{N}}
    return R
end

# convert to a reach-set with another underlying set representation
function Base.convert(T::Type{ReachSet{N, LT}}, R::ReachSet{N, ST}) where {N, ST<:LazySet{N}, LT<:LazySet{N}}
    X = convert(LT, set(R))
    return ReachSet(X, tspan(R))
end

function shift(R::ReachSet, t0::Number)
    return ReachSet(set(R), tspan(R) + t0)
end

function reconstruct(R::ReachSet, Y::LazySet)
    return ReachSet(Y, tspan(R))
end

# constructor with a time point
function ReachSet(X::ST, t::Real) where {N, ST<:LazySet{N}}
    return ReachSet(X, interval(t))
end

#=
# specialized concrete projection of zonotopic set by removing zero generators
function linear_map(R::ReachSet{N, <:AbstractZonotope{N}}, variables::NTuple{D, M}; check_vars::Bool=true) where {D, M<:Integer}
    πR = linear(R, variables, check_vars)
    X = remove_zero_generators(set(πR))
    SparseReachSet(proj, tspan(R), variables)
end

function project(R::ReachSet{N, <:AbstractZonotope{N}}, vars) where {N}
    πR = linear_map(M, R)
    dt = tspan(R)
    if isnothing(vars)
        X = remove_zero_generators(set(πR))
        v = ReachabilityAnalysis.vars(R)
    else
        Y = project(πR, vars)
        X = remove_zero_generators(set(Y))
        v = ReachabilityAnalysis.vars(Y)
    end
    return SparseReachSet(X, dt, v)
end
=#


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import LibGit2: AbstractCredential, UserPasswordCredential, SSHCredential,
    CachedCredentials, CredentialPayload, Payload
using Base: coalesce

const DEFAULT_PAYLOAD = CredentialPayload(allow_ssh_agent=false, allow_git_helpers=false)

"""
Emulates the LibGit2 credential loop to allows testing of the credential_callback function
without having to authenticate against a real server.
"""
function credential_loop(
        valid_credential::AbstractCredential,
        url::AbstractString,
        user::Union{AbstractString, Nothing},
        allowed_types::UInt32,
        payload::CredentialPayload;
        shred::Bool=true)
    cb = LibGit2.credentials_cb()
    libgitcred_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)

    # Number of times credentials were authenticated against. With the real LibGit2
    # credential loop this would be how many times we sent credentials to the remote.
    num_authentications = 0

    # Emulate how LibGit2 uses the credential callback by repeatedly calling the function
    # until we find valid credentials or an exception is raised.
    err = Cint(0)
    while err == 0
        err = ccall(cb, Cint, (Ptr{Ptr{Cvoid}}, Cstring, Cstring, Cuint, Any),
                    libgitcred_ptr_ptr, url, coalesce(user, C_NULL),
                    allowed_types, payload)
        num_authentications += 1

        # Check if the callback provided us with valid credentials
        if payload.credential !== nothing && payload.credential == valid_credential
            LibGit2.approve(payload, shred=shred)
            break
        end

        if num_authentications > 50
            error("Credential callback seems to be caught in an infinite loop")
        end
    end

    # Note: LibGit2.GitError(0) will not work if an error message has been set.
    git_error = if err == 0
        LibGit2.GitError(LibGit2.Error.None, LibGit2.Error.GIT_OK, "No errors")
    else
        LibGit2.GitError(err)
    end

    # Reject and shred the credential when an authentication error occurs
    if git_error.code == LibGit2.Error.EAUTH
        LibGit2.reject(payload, shred=shred)
    end

    return git_error, num_authentications, payload
end

function credential_loop(
        valid_credential::UserPasswordCredential,
        url::AbstractString,
        user::Union{AbstractString, Nothing}=nothing,
        payload::CredentialPayload=DEFAULT_PAYLOAD;
        shred::Bool=true)
    credential_loop(valid_credential, url, user, 0x000001, payload, shred=shred)
end

function credential_loop(
        valid_credential::SSHCredential,
        url::AbstractString,
        user::Union{AbstractString, Nothing}=nothing,
        payload::CredentialPayload=DEFAULT_PAYLOAD;
        shred::Bool=true)
    credential_loop(valid_credential, url, user, 0x000046, payload, shred=shred)
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
using CompScienceMeshes
using Test

m = meshsegment(1.0, 1.0, 3)

p1 = point(0.5, 0, 0)
p2 = point(0.5, 1e-4, 0)
p3 = point(0.0, 0, 0)
p4 = point(-1e-4, 0, 0)

f = inclosure_gpredicate(m)

##
@test f(p1) == true
@test f(p2) == false
@test f(p3) == true
@test f(p4) == false


text:
using ARules
using StatsBase
using Test

# write your own tests here
itemlist = randstr(100, 16);

n = 100_000
m = 20               # NOTE: most impactful for runtime complexity (num. items in transactions)
mx_depth = 10        # max depth of itemset tree (max size of transactions explored)
t = [sample(itemlist, m, replace = false) for _ in 1:n];

unq = unique_items(t);
@test typeof(unq) == Array{String, 1}

occ = occurrence(t, unq);
@test typeof(occ) == BitArray{2}

xtree1 = frequent_item_tree(t, unq, round(Int, 0.01*n), mx_depth);
@test typeof(xtree1) == Node

rules = apriori(t, supp = 0.01, conf = 0.1, maxlen = mx_depth);

t2 = [["a", "b"],
     ["b", "c", "d"],
     ["a", "c"],
     ["e", "b"],
     ["a", "c", "d"],
     ["a", "e"],
     ["a", "b", "c"],
     ["c", "b", "e", "f"]]

unq2 = unique_items(t2)
@test length(unq2) == 6

xtree2 = frequent_item_tree(t2, unq2, 1, 4);
xsup = gen_support_dict(xtree2, length(t2))
@test length(xsup) == 27
@test typeof(xsup) == Dict{Array{Int16,1}, Int64}

xrules = gen_node_rules(xtree2.children[1].children[1].children[1], xsup, 3, 8, 0.1)
@test length(xrules) == 3

rule_arr = Array{Rule, 1}(undef, 0)
gen_rules!(rule_arr, xtree2.children[1], xsup, 2, 8, 0.1)
@test eltype(rule_arr) == Rule
@test length(rule_arr) == 14

# A full run through
a_list = [
    ["a", "b"],
    ["a", "c"],
    ["a", "b", "c"],
    ["a", "b", "d"],
    ["a", "c", "d"],
    ["a", "b", "c", "d"],
    ["a", "b", "c", "e"],
    ["b", "d", "e", "f"],
    ["a", "c", "e", "f"],
    ["b", "c", "d", "e", "f"],
    ["a", "c", "d", "e", "f"],
    ["b", "c", "d", "e", "f"]
]

xunq3 = unique_items(a_list)
xtree3 = frequent_item_tree(a_list, xunq3, 1, 6);
xsup3 = gen_support_dict(xtree3, length(a_list))

rules = apriori(a_list, supp = 0.01, conf = 0.01, maxlen = 6)

transactions = [["milk", "eggs", "bread"],
			    ["butter", "milk", "sugar", "flour", "eggs"],
			    ["bacon", "eggs", "milk", "beer"],
			    ["bread", "ham", "turkey"],
			    ["cheese", "ham", "bread", "ketchup"],
			    ["mustard", "hot dogs", "buns", "hamburger", "cheese", "beer"]]

freq = frequent(transactions, 1, 8)

rules = apriori(transactions, supp = 0.01, conf = 0.01, maxlen = 6)
@test size(rules) == (329, 5)

unq = unique_items(transactions);
occ = occurrence(transactions, unq);
rules = apriori(occ, supp = 0.01, conf = 0.01, maxlen = 6)
@test length(rules) == 329


text:
function lp1_test()
    #=
        min -4x_2 - 1
    s.t.
        x_1 + 2x_2 <= 3
        x_1 >= 3
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0], X), 0.0),
            MOI.LessThan(3.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
            MOI.GreaterThan(3.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0], [X[2]]), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp2_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2 + 1 <= 4
        x1 + 2x2 + 1 <= 4
        x1 >= 1
        x2 >= 0
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2.0, 1.0], X), 0.0),
             MOI.LessThan(3.0))

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0], X), 0.0),
             MOI.LessThan(3.0))

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(1.0, X[1])], 0.0),
             MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(1.0, X[2])], 0.0),
             MOI.GreaterThan(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], X), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp3_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2 + 1 <= 4
        x1 + 2x2 + 1 <= 4
        x1 >= 1
        x2 >= 0
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2.0, 1.0], X), 0.0),
             MOI.LessThan(3.0))

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0], X), 0.0),
             MOI.LessThan(3.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
             MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
             MOI.GreaterThan(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], X), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp4_test()
    #=
        max 4x1 + 3x2
      s.t.
        x1 >= 1
        x2 >= 0
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
             MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
             MOI.GreaterThan(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0], X), 0.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    return model
end

function lp5_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2  == 3
        x1 + 2x2  == 3
        x1 >= 1
        x2 == 0
    =#
    model= TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2.0, 1.0], X), 0.0),
            MOI.EqualTo(3.0))

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0], X), 0.0),
            MOI.EqualTo(3.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
            MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
            MOI.EqualTo(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], X), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp6_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2 - 3 <= 0
        x1 + 2x2 - 3 <= 0
        x1 >= 1
        x2 >= 0
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    c1 = MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.([1, 1, 2, 2], MOI.ScalarAffineTerm.([2.0, 1.0, 1.0, 2.0], [X; X])),
                [-3.0, -3.0])

    MOI.add_constraint(model, c1, MOI.Nonpositives(2))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
            MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
            MOI.GreaterThan(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], X), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp7_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2 - 3 <= 0
        x1 + 2x2 - 3 <= 0
        x1 >= 1
        x2 >= 0
    =#
    model= TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    c1 = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1, 1, 2, 2], MOI.ScalarAffineTerm.([2.0, 1.0, 1.0, 2.0], [X; X])),
            [-3.0, -3.0])

    MOI.add_constraint(model, c1, MOI.Nonpositives(2))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
            MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
            MOI.GreaterThan(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], X), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp8_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2 - 3 <= 0
        x1 + 2x2 - 3 <= 0
        x1 >= 1
        x2 >= 0
    =#
    model = TestModel{Int64}()

    X = MOI.add_variables(model, 2)

    c1 = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1, 1, 2, 2], MOI.ScalarAffineTerm.([2, 1, 1, 2], [X; X])),
            [-3, -3])

    MOI.add_constraint(model, c1, MOI.Nonpositives(2))

    c2 = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([1, 1], X)),
            [1, 0])

    MOI.add_constraint(model, c2, MOI.Nonnegatives(2))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4, -3], X), -1))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp9_test()
    #=
        min x + y
    s.t.
        -1 <= x + y <= 10
        x,  y >= 0
    =#
    model = TestModel{Float64}()

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)

    vc = MOI.add_constraints(model,
        [MOI.SingleVariable(x), MOI.SingleVariable(y)],
        [MOI.GreaterThan(0.0), MOI.GreaterThan(0.0)])

    c = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x,y]), 0.0), MOI.Interval(-1.0, 10.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp10_test()
    #=
        min x1
    s.t.
        2x1 + x2  == 3
        x1 + 2x2  == 3
        x1 >= 1
        x2 == 0
    =#
    model= TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2.0, 1.0], X), 0.0),
            MOI.EqualTo(3.0))

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0], X), 0.0),
            MOI.EqualTo(3.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
            MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
            MOI.EqualTo(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(X[1]))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp11_test()
    #=
        Feasibility
    s.t.
        x1 >= 1
        x2 >= 0
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
             MOI.GreaterThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
             MOI.GreaterThan(0.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0], X), 0.0))

    return model
end

function lp12_test()
    #=
        min 4x_3 + 5
    s.t.
        x_1 + 2x_2 + x_3 <= 20
        x_1 <= 1
        x_2 <= 3
    =#
    model = TestModel{Float64}()

    X = MOI.add_variables(model, 3)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0, 1.0], X), 0.0),
            MOI.LessThan(20.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[1]),
            MOI.LessThan(1.0))

    MOI.add_constraint(model,
        MOI.SingleVariable(X[2]),
            MOI.LessThan(3.0))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0], [X[3]]), 5.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function lp13_test()
    #=
        min -4x1 -3x2 -1
    s.t.
        2x1 + x2 - 3 <= 0
        x1 + 2x2 - 3 <= 0
        x1 >= 0
        x2 >= 0
    =#
    model= TestModel{Float64}()

    X = MOI.add_variables(model, 2)

    c1 = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1, 1, 2, 2], MOI.ScalarAffineTerm.([2.0, 1.0, 1.0, 2.0], [X; X])),
            [-3.0, -3.0])

    MOI.add_constraint(model, c1, MOI.Nonpositives(2))

    MOI.add_constraint(model, MOI.VectorOfVariables(X), MOI.Nonnegatives(2))

    MOI.set(model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], X), -1.0))

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end


text:
function NakagamiDistribution(arg0::jdouble, arg1::jdouble)
    return NakagamiDistribution((jdouble, jdouble), arg0, arg1)
end

function NakagamiDistribution(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return NakagamiDistribution((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function cumulative_probability(obj::NakagamiDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::NakagamiDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_numerical_mean(obj::NakagamiDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::NakagamiDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_scale(obj::NakagamiDistribution)
    return jcall(obj, "getScale", jdouble, ())
end

function get_shape(obj::NakagamiDistribution)
    return jcall(obj, "getShape", jdouble, ())
end

function get_support_lower_bound(obj::NakagamiDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::NakagamiDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function is_support_connected(obj::NakagamiDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end


text:
function trs_small(P::AbstractMatrix{T}, q::AbstractVector{T}, r::T; kwargs...) where {T}
	output = trs_boundary_small(P, q, r; kwargs...)
	return check_interior!(output..., P, q; direct=true)
end

function trs_small(P::AbstractMatrix{T}, q::AbstractVector{T}, r::T, C::AbstractMatrix{T}; kwargs...) where {T}
	output = trs_boundary_small(P, q, r, C; kwargs...)
	return check_interior!(output..., P, q; direct=true)
end

function trs_boundary_small(P::AbstractMatrix{T}, q::AbstractVector{T}, r::T, C::AbstractMatrix{T}; kwargs...) where T
	check_inputs(P, q, r, C)
	return trs_boundary((nev; kw...) -> gen_eigenproblem_small(P, q, r, C, nev; kw...),
		   (λ, V; kw...) -> pop_solution!(λ, V, P, q, r, C; kw..., direct=true); kwargs...)
end

function trs_boundary_small(P::AbstractMatrix{T}, q::AbstractVector{T}, r::T; kwargs...) where {T}
	check_inputs(P, q, r)
	return trs_boundary((nev; kw...) -> eigenproblem_small(P, q, r, nev; kw...),
		   (λ, V; kw...) -> pop_solution!(λ, V, P, q, r, I; kw..., direct=true); kwargs...)
end


text:
@testset "Sanity checks" begin
    @test isempty(detect_unbound_args(Simulation73))
end


text:
export windprofile_logLL, windprofile_logL, windprofile_DH,
      windprofile_PL_α, windprofile_PL_α_P, windprofile_PL_n_JM,
      windprofile_PL_n_SR

include("constant.jl")
#using .constant

# Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
# Renewable and Sustainable Energy Reviews 102 (2019): 215-233.

# Monin-obukhov length
"""

    Monin_Obukhov_length(u✶,H)

    Monin_Obukhov_length(u✶,H,T)

calculate the Monin-Obukhov length

# Arguments

- `u✶`: m.s-1, friction velocity. `u✶` for light, moderate and strong wind conditions
can be selected as 0.1, 0.3 and 0.6 m/s respectively

- `H`: W.m-2, sensible heat flux. For mostly sunny conditions during the daytime,
`H` may be set 200 W.m-2; for cloudy conditions, day or night, `H` ~ 0 W.m-2;
for clear skies at night `H` ~ -4 W.m-2 for light wind conditions, and
  `H` ~ -20 W.m-2 for moderate or strong wind.

- `T`: ᵒC, air temperature, default `T` = 20 ᵒC

# Example:
```julia-repl
julia> Monin_Obukhov_length(0.3,-20)
```
# Reference
[1]. Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
"""
function Monin_Obukhov_length(u✶,H,T=20)
  T += Kelvin

  L = (-ϱ0*cP*T*u✶^3) / (κ*g0*H)

  return L
end

# correction under non-neutral conditions
"""

  𝛹m(z/L)

correct deviations from the wind profile under non-neutral conditions

# Argument

- `z/L`: where `L` is the Monin-Obukhov length
"""
function 𝛹m(x)
      if x>0 # stable conditions
        return -5*x
      else # unstable conditions
        x2 = ( 1 - 15*(x) )^(1/4)
        value = 2*log( (1+x2)/2 ) + log( (1+x2^2)/2 ) - 2*atan(x2) + π/2
        return value
    end
end

"""

    𝛷m(z/L)

correct deviations from the wind profile under non-neutral conditions

# Argument

- `zg/L`: where `L` is the Monin-Obukhov length

"""
function 𝛷m(x)
  if x>0
    return 1+5*x
  else
    return (1 - 15*x)^(-1/4)
  end
end

"""

  Coriolis(lat)

calculate Coriolis frequency, f

# Arguments

- `lat`: deg, latitude

# Example
```julia-repl
julia> Coriolis(45)
```
"""
function Coriolis(lat)
  ω = 7292115*10^-11 # rad/s, rotation rate of the Earth
  f = 2*ω*sin(lat*π/180) # Coriolis parameter
  return f
end

"""

    windprofile_logLL(h,z0,vref,href)

    windprofile_logLL(h,z0,vref,href, u✶, H)

calculate wind speed at a given height using Log-linear law

# Arguments

- `h`: m, height need to extrapolate wind speed

- `z0`: m, roughness length

- `vref`: m/s, wind speed at reference height

- `href`: m, reference height

- `u✶`: m.s-1, friction velocity, default `u✶` = 0.3 for moderate wind condition

- `H`: W.m-2, sensible heat flux, default `H` = -20 for clear sky at night with
moderate wind condition

# Example
```julia-repl
julia> windprofile_logLL(10,0.05,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.

"""
function windprofile_logLL(h,z0,vref,href, u✶=0.3, H=-20)

    L = Monin_Obukhov_length(u✶,H)

    vh = vref*( log(h/z0) - 𝛹m(h/L) ) / ( log(href/z0) - 𝛹m(href/L) )

    return vh

end

"""

    windprofile_logL(h,z0,vref,href)

calculate wind speed at a given height using Log-linear law (neutral conditions)

# Arguments

- `h`: m, height need to extrapolate wind speed

- `z0`: m, roughness length

- `vref`: m/s, wind speed at reference height

- `href`: m, reference height

# Example
```julia-repl
julia> windprofile_logL(10,0.05,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.

"""
function windprofile_logL(h,z0,vref,href)
  vh = vref* log(h/z0) / log(href/z0)
  return vh
end

# Eaves and Harris model
"""

    windprofile_DH(h,z0,vref,href)

    windprofile_DH(h,z0,vref,href,u✶,lat)

calculate wind speed at a given height using Deaves and Harris model
(near neutral strong wind conditions)

# Arguments

- `h`: m, height need to extrapolate wind speed

- `z0`: m, roughness length

- `vref`: m/s, wind speed at reference height

- `href`: m, reference height

- `u✶`: m.s-1, friction velocity, default `u✶` = 0.3 for moderate wind condition

- `lat`: deg, latitude in degree, default `lat` = 30 for South Australia

# Example
```julia-repl
julia> windprofile_DH(10,0.05,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.
"""
function windprofile_DH(h,z0,vref,href, u✶=0.3,lat=30)

  f = Coriolis(lat)   # Coriolis parameter
  B =  6.0            # empirical constant

  hf = u✶/(B*f) # neutral ABL height

  vh = vref* ( log(h/z0) + 5.75*(h/hf) - 1.88*(h/hf)^2 - 1.33*(h/hf)^3 + 0.25*(h/hf)^4 ) /
     ( log(href/z0) + 5.75*(href/hf) - 1.88*(href/hf)^2 - 1.33*(href/hf)^3 + 0.25*(href/hf)^4 )

  return vh

end

# Power law
"""

    windprofile_PL_α(h,α,vref,href)

calculate wind speed at a given height using power law model

# Arguments

- `h`: m, height need to extrapolate wind speed

- `α`: dimensionless exponent

- `vref`: m/s, wind speed at reference height

- `href`: m, reference height

# Example
```julia-repl
julia> windprofile_PL_α(10,1/7,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.

"""
function windprofile_PL_α(h,α,vref,href)
  vh = vref*(h/href)^(α)

  return vh
end

"""

    windprofile_PL_α_P(h,α,vref,href)

    windprofile_PL_α_P(h,z0,vref,href, u✶ H)

calculate wind speed at a given height using MOST

# Arguments

- `h`: m, height need to extrapolate wind speed

- `α`: dimensionless exponent

- `vref`: m/s, wind speed at reference height

- `href`: m, reference height

# Example
```julia-repl
julia> windprofile_PL_α_P(10,1/7,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.
"""
function windprofile_PL_α_P(h,z0,vref,href, u✶=0.3, H=-20)
  zg = (h*href)^0.5

  L = Monin_Obukhov_length(u✶,H)

  α = 𝛷m(zg/L) / ( log(zg/z0) - 𝛹m(zg/L) )

  vh = vref*(h/href)^α

  return vh

end

"""

    windprofile_PL_n_JM(z2,z1,c1,k1)

    windprofile_PL_n_JM(z2,z1,c1,k1)

calculate wind speed at a given height using Weibull probability function

# Arguments

- `z2`: m, height need to extrapolate wind speed

- `z1`: height at available data

- `c1`: Weibull parameter

- `k1`: Weibull parameter

# Example
```julia-repl
julia> windprofile_PL_α_P(10,1/7,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.

"""
function windprofile_PL_n_JM(z2,z1,c1,k1)

  zr = 10 # reference height
  k2 = k1*( 1 - 0.0881*log(z1/zr) ) / ( 1 - 0.0881*log(z2/zr) )

  n = (0.37 - 0.0881*log(c1)) / ( 1 - 0.0881*log(z1/zr) )

  c2 = c1*(z2/z1)^n

  return c2, k2

end

"""

    windprofile_PL_n_SR(z2,z1,c1,k1)

    windprofile_PL_n_SR(z2,z1,c1,k1)

calculate wind speed at a given height using Weibull probability function

# Arguments

- `z2`: m, height need to extrapolate wind speed

- `z1`: height at available data

- `c1`: Weibull parameter

- `k1`: Weibull parameter

# Example
```julia-repl
julia> windprofile_PL_N_SR(10,1/7,3,1.5)
```

# Reference
[1]. Gualtieri, Giovanni. "A comprehensive review on wind resource extrapolation models applied in wind energy."
Renewable and Sustainable Energy Reviews 102 (2019): 215-233.

"""
function windprofile_PL_n_SR(z2,z1,c1,k1,z0)

  zr = 10 # reference height
  vh = 67 # m/s
  α0 = (z0/zr)^0.2

  k2 = k1*( 1 - 0.0881*log(z1/zr) ) / ( 1 - 0.0881*log(z2/zr) )

  n = α0*( 1 - log(c1)/log(vh) ) / ( 1 - α0*log(z1/zr)/log(vh) )

  c2 = c1*(z2/z1)^n

  return c2, k2

end


text:
import Base: size, ∘

"""
    Nodes{Dual/Primal}

`Nodes` is a wrapper for scalar-valued data that lie at the centers of either dual cells or
primary cells. A `Nodes` type can be accessed by indexing like any other array,
and allows the use of [size].

# Constructors
- `Nodes(C,dims)` creates a field of zeros in cells of type `C` (where `C` is
  either `Dual` or `Primal`), on a grid of dimensions `dims`. Note that `dims`
  represent the number of dual cells on the grid, even if `C` is `Primal`.
- `Nodes(C,w)` performs the same construction, but uses existing field data `w`
  of `Nodes` type to determine the size of the grid.
"""
struct Nodes{C <: CellType, NX, NY, T <: Number} <: ScalarGridData{NX,NY,T}
    data::Matrix{T}
end

# Number of indices
# Based on number of dual nodes, return the number of nodes
node_inds(::Type{Dual},   dualnodedims) = (dualnodedims[1], dualnodedims[2])
node_inds(::Type{Primal}, dualnodedims) = (dualnodedims[1]-1, dualnodedims[2]-1)

# Constructors

function Nodes(T::Type{C}, dualnodedims::Tuple{Int, Int};dtype=Float64) where {C <: CellType}
    dims = node_inds(T, dualnodedims)
    Nodes{T, dualnodedims...,dtype}(zeros(dtype,dims))
end

# This allows easy construction of nodes of either type from existing nodes of either
# type on the same grid.
Nodes(C, ::GridData{NX,NY,T};dtype=T) where {NX, NY,T <: Number} = Nodes(C, (NX, NY),dtype=dtype )

Nodes(C, nx::Int, ny::Int;dtype=Float64) = Nodes(C,(nx,ny),dtype=dtype)
(::Type{Nodes{C,NX,NY,T}})() where {C,NX,NY,T} = Nodes(C, (NX, NY),dtype=T)

Base.similar(::Nodes{C,NX,NY,T};element_type=T) where {C,NX,NY,T} = Nodes(C, (NX, NY),dtype=element_type)

function Base.show(io::IO, nodes::Nodes{C, NX, NY,T}) where {C, NX, NY, T}
    nodedims = "(nx = $NX, ny = $NY)"
    dims = "(nx = $(size(nodes,1)), ny = $(size(nodes,2)))"
    println(io, "$C nodes in a $nodedims cell grid of type $T data")
    print(io, "  Number of $C nodes: $dims")
end


text:
# NbodyGradient integration to compute points about each transit for series expansion.

abstract type AbstractTransitTimes end
struct ComputedTimes <: AbstractTransitTimes end
struct ProvidedTimes <: AbstractTransitTimes end

struct TransitSeries{T<:Real, TT<:AbstractTransitTimes}
    times::Vector{T}      # Transit times, sequentially
    bodies::Vector{Int64} # Transiting body at the ith transit time
    points::Array{T, 4}   # [body, transit, 7, 2]
    dpoints::Array{T, 6}  # [body, transit, 7, 2, body, parameter]

    # Internal values/arrays
    h::T                  # Stepsize for series points
    ntt::Int64            # Total number of transits
    intr_times::Vector{T} # Times to integrate to anc compute points
    count::Vector{Int64}  # Number of transits for ith body
    s_prior::State{T}
end

"""Pass a set of transit times"""
function TransitSeries(times::Matrix{T}, ic::InitialConditions{T}; h::T=T(2e-2)) where T<:Real
    ntt = sum(times .> ic.t0) # Total transits, assumes no transits <= ic.t0

    # Flatten transit times and assign a body index
    count = zeros(Int64, ic.nbody)
    tmp = zeros(T,ntt,2)
    for (i,t) in enumerate(eachrow(times))
        for j in eachindex(t)
            tc = t[j]
            if tc > ic.t0
                count[i] += 1
                ind = sum(count)
                tmp[ind,1] = tc
                tmp[ind,2] = i
            end
        end
    end

    # Now sort by transit time
    tmp .= sortslices(tmp, dims=1)

    times = copy(tmp[:,1])
    bodies = copy(Int64.(tmp[:,2]))

    # Compute the time of the first expansion point for each transit
    intr_times = zeros(T,length(times))
    for (i,t) in enumerate(times)
        intr_times[i] =  t - 3*h
    end

    points = zeros(T,ic.nbody,ntt,7,2)
    dpoints = zeros(T,ic.nbody,ntt,7,2,ic.nbody,7)
    s_prior = State(ic)
    TransitSeries{T, ProvidedTimes}(times, bodies, points, dpoints, h, ntt, intr_times, count, s_prior)
end

"""Need to compute transit times"""
function TransitSeries(tmax::T, ic::InitialConditions{T}; h::T=T(2e-2)) where T<:Real
    ntt::Int64 = 0  # Total expected transits, assumes no transits <= ic.t0
    for P in ic.elements[:,2]
        ntt += round(Int64, tmax/P)
    end

    times = zeros(T, ntt)
    bodies = zeros(Int64, ntt)
    intr_times = zeros(T,length(times))
    points = zeros(T,ic.nbody,ntt,7,2)
    dpoints = zeros(T,ic.nbody,ntt,7,2,ic.nbody,7)
    s_prior = State(ic)
    TransitSeries{T, ComputedTimes}(times, bodies, points, dpoints, h, ntt, intr_times, count, s_prior)
end

"""Compute 7 points about each transit time.

Integrate to right before the first transit expansion point. Save state and
integrate for 7 steps at h=ts.h. Revert to pre-transit state and continue to
next transit.
"""
function (intr::Integrator)(s::State{T},ts::TransitSeries{T, ProvidedTimes},d::Union{Derivatives,Nothing}=nothing; grad::Bool=false) where T<:Real
    if grad && d == nothing; d = Derivatives(T, s.n); end
    if ~grad && d !== nothing; grad = true; end

    # Run integrator and record sky positions for list of integration times
    nstep::Int64 = 0; t0 = s.t[1]
    for (itime,t) in enumerate(ts.intr_times)
        ibody = ts.bodies[itime]

        # Integrate to right before the first point in the transit
        while s.t[1] < t
            # Break if we will pass a transit on this step
            if t-s.t[1] < intr.h; break; end

            if grad
                intr.scheme(s,d,intr.h)
            else
                intr.scheme(s,intr.h)
            end
            nstep += 1
            s.t[1] = t0 + nstep*intr.h
        end

        # Save the pre-transit state
        set_state!(ts.s_prior, s)

        # Now integrate to each expansion point with step ts.h
        for k in 1:7

            # Step to the first point if needed
            k == 1 ? h_points = t-ts.s_prior.t[1] : h_points = ts.h

            if grad
                intr.scheme(s,d,h_points)
            else
                intr.scheme(s,h_points)
            end

            # record points for transiting body
            ts.points[ibody,itime,k,1] = s.x[1,ibody]
            ts.points[ibody,itime,k,2] = s.x[2,ibody]

            if grad
                # Get gradients of points wrt initial orbital elements and masses
                # p body, q element
                for p in 1:s.n, q in 1:7
                    ts.dpoints[ibody,itime,k,1,p,q] = s.jac_step[(ibody-1)*7+1,(p-1)*7+q]
                    ts.dpoints[ibody,itime,k,2,p,q] = s.jac_step[(ibody-1)*7+2,(p-1)*7+q]
                end
            end
        end

        # Return to state before transit
        set_state!(s, ts.s_prior)
    end
end


text:
# # Logistic regression

import MaximumLikelihoodProblems

import Distributions
import ForwardDiff
import LogDensityProblems
import StatsFuns
import TransformVariables

struct LogisticRegression{Ty, TX}
    y::Ty
    X::TX
end

function (problem::LogisticRegression)(θ)
    y = problem.y
    X = problem.X

    β = θ.β

    η = X * β
    μ = StatsFuns.logistic.(η)
    log_likelihood = sum(Distributions.logpdf.(Distributions.Bernoulli.(μ), y))
    return log_likelihood
end

N = 10_000

## the first column (the column of all ones) is the intercept
X = hcat(ones(N), randn(N))

size_β = (2,)
β_true = [1.0, 2.0]
η_true = X * β_true
μ_true = StatsFuns.logistic.(η_true)
y = rand.(Distributions.Bernoulli.(μ_true))

problem = LogisticRegression(y, X)
transformation = TransformVariables.as((β = TransformVariables.as(Array, size_β), ))
transformed_problem = LogDensityProblems.TransformedLogDensity(transformation,
                                                               problem)
transformed_gradient_problem = LogDensityProblems.ADgradient(:ForwardDiff,
                                                             transformed_problem)

β_hat_initial_guess = zeros(size_β)
θ_hat_initial = (; β = β_hat_initial_guess)

# θ_hat:

θ_hat = MaximumLikelihoodProblems.fit(transformed_gradient_problem,
                                      θ_hat_initial)

# β_hat:

β_hat = θ_hat[:β]

# Value of the log likelihood function evaluated at θ_hat:

MaximumLikelihoodProblems.loglikelihood(transformed_gradient_problem, θ_hat)


text:
# ================
# Homogeneous case
# ================

function reach_homog!(R::Vector{ReachSet{<:Zonotope{Float64}}},
                      Ω0::Zonotope,
                      Φ::AbstractMatrix,
                      N::Int,
                      δ::Float64)
    # initial reach set
    t0, t1 = zero(δ), δ
    R[1] = ReachSet(Ω0, t0, t1)

    k = 2
    while k <= N
        Rₖ = linear_map(Φ, set(R[k-1]))
        t0 = t1
        t1 += δ
        R[k] = ReachSet(Rₖ, t0, t1)

        k += 1
    end
    return R
end

# ==================
# Inhomogeneous case
# ==================

function reach_inhomog!(R::Vector{ReachSet{<:Zonotope{Float64}}},
                        Ω0::Zonotope,
                        U::ConstantInput,
                        Φ::AbstractMatrix,
                        N::Int,
                        δ::Float64,
                        max_order::Int)
    # initial reach set
    t0, t1 = zero(δ), δ
    R[1] = ReachSet(Ω0, t0, t1)

    V = next_set(U)
    Wk₊ = V
    Φ_power_k = copy(Φ)
    Φ_power_k_cache = similar(Φ)

    k = 2
    while k <= N
        Rₖ = minkowski_sum(linear_map(Φ_power_k, Ω0), Wk₊)
        Rₖ = reduce_order(Rₖ, max_order)
        t0 = t1
        t1 += δ
        R[k] = ReachSet(Rₖ, t0, t1)

        Wk₊ = minkowski_sum(Wk₊, linear_map(Φ_power_k, V))
        Wk₊ = reduce_order(Wk₊, max_order)

        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
        k += 1
    end
    return R
end


text:
module Geometry

export centroid, squaredradius
"""
Computes the centroid of a simplex defined by vertices in 'X',
which is a matrix of dimension nx(n+1) and each column is a vertex.

Input arguments
---------------
simplex::Array{Float64, 2}  Simplex represented by an n-by-(n+1) matrix, where each column is a vertex.

Returns
-------
centroid::Array{Float64, 1}   Row vector. The centroid of the simplex.
"""
function centroid(simplex::AbstractArray{Float64, 2})
    n = size(simplex)[1]
    centroid = simplex * ones(n + 1, 1)/(n + 1)
    return(centroid)
end

"""
Computes the radius of a simplex defined by vertices in 'X',
which is a matrix of dimension nx(n+1) and each column is a vertex.

Input arguments
---------------
simplex::Array{Float64, 2}  Simplex represented by an n-by-(n+1) matrix, where each column is a vertex.

Returns
-------
radius::Float64  The radius of the simplex.

"""
function squaredradius(simplex)
    n = size(simplex)[1]
    C1 = repeat(centroid(simplex), 1, (n + 1))
    radius = maximum(ones(1, n) * ((simplex - C1)))^2
    return(radius::Float64)
end

end


text:
using TimeseriesSurrogates

export bbnue

# Keep this for when a common interface for optimized variable selection methods has been established
# """
#export BBNUE
#     BBNUE(est) <: TransferEntropyEstimator

# The bootstrap-based non-uniform embedding estimator (BB-NUE) for conditional transfer entropy (Montalto et al., 2014).
# Uses the estimator `est` to compute relevant marginal entropies. (e.g. `VisitationFrequency(RectangularBinning(3))`)

# [^Montalto2014]: Montalto, A.; Faes, L.; Marinazzo, D. MuTE: A MATLAB toolbox to compare established and novel estimators of the multivariate transfer entropy. PLoS ONE 2014, 9, e109462.
# """
# struct BBNUE{E} <: TransferEntropyEstimator
#     est::E
# end

"""
    bbnue(source, target, [cond], est::BBNUE; q = 0.95, η = 1,
        nsurr = 100, uq = 0.95,
        include_instantaneous = true,
        method_delay = "ac_min",
        maxlag::Union{Int, Float64} = 0.05
        ) → te, js, τs, idxs_source, idxs_target, idxs_cond

Estimate transfer entropy using the bootstrap-based non-uniform embedding (BBNUE) estimator,
which uses a bootstrap-basedcriterion to identify the most relevant and minimally redundant
variables from the present/past of `source`, present/past `cond` (if given) and the past of
`target` that contribute most to `target`'s future. `η` is the forward prediction lag.
Multivariate `source`, `target` and `cond` (if given) are all possible.

For significance testing of a variable, `nsurr` circular shift surrogates are generated,
and if transfer entropy for the original variables exceeds the `uq`-quantile of that of the
surrogate ensemble, then the variable is included.

If `instantaneous` is `true`, then instantaneous interactions are also considered, i.e. effects like
`source(t) → target(t)` are allowed.

In this implementation, the maximum lag for each embedding variable is determined using `estimate_delay`
from `DelayEmbeddings`. The keywords `method_delay` (default is "ac_min") controls the method
for estimating the delay, and `maxlag` is the maximum allowed delay (if `maxlag ∈ [0, 1]` is a fraction,
then the maximum lag is that fraction of the input time series length, and if `maxlag` is an integer,
then the maximum lag is `maxlag`).

## Implementation details

Currently, only this implementation is optimized for the bin-estimator approach from
Montalto et al. (2014)[^Montalto2014], which uses a conditional entropy minimization criterion for
select variables.  Their bin-estimator approach corresponds to using
the `VisitationFrequency` estimator with bins whose sides are equal-length, e.g.
`VisitationFrequency(RectangularBinning(0.5))`. Here, you can use any desired rectangular binning.

It is also possible to use other entropy estimators than `VisitationFrequency` to compute entropies,
but this implementation will use conditional entropy minimization regardless of the choice of
estimator.

## Returns

A 6-tuple is returned, consisting of:
- `te`: The computed transfer entropy value. If no relevant variables were selected, then `te = 0.0`.
- `js`: The indices of the selected variables. `js[i]` is the `i`-th entry in the array `[idxs_source..., idxs_target..., idxs_cond...,]`.
- `τs`: The embedding lags of the selected variables. `τs[i]` corresponds to `js[i]`.
- `idxs_source`: The indices of the source variables.
- `idxs_target`: The indices of the target variables.
- `idxs_cond`: The indices of the conditional variables (empty if `cond` is not given).

## Example

```julia
using CausalityTools, DynamicalSystems
sys = ExampleSystems.logistic2_unidir(c_xy = 1.5)
orbit = trajectory(sys, 10000, Ttr = 10000)
x, y = columns(orbit)

# Use a coarse-grained rectangular binning with subdivisions in each dimension,
# to keep computation costs low and to ensure the probability distributions
# over the bins don't approach the uniform distribution (need enough points
# to fill bins).
est = VisitationFrequency(RectangularBinning(3))
te_xy, params_xy = bbnue(x, y, BBNUE(est))
te_yx, params_yx = bbnue(y, x, BBNUE(est))

te_xy, te_yx
```

[^Montalto2014]: Montalto, A.; Faes, L.; Marinazzo, D. MuTE: A MATLAB toolbox to compare established and novel estimators of the multivariate transfer entropy. PLoS ONE 2014, 9, e109462.
"""
function bbnue(source, target, cond, est; q = 0.95, η = 1, nsurr = 100, uq = 0.95,
        include_instantaneous = true, method_delay = "ac_min", maxlag::Union{Int, Float64} = 0.05)

    Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond =
        embed_candidate_variables(
            process_input(source),
            process_input(target),
            process_input(cond),
            η = η)

    return optim_te(Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond, est, q = q, nsurr = 19, uq = uq)
end

function bbnue(source, target, est; q = 0.95, η = 1, nsurr = 100, uq = 0.95,
        include_instantaneous = true, method_delay = "ac_min", maxlag::Union{Int, Float64} = 0.05)

    Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond =
        embed_candidate_variables(process_input(source), process_input(target), η = η)

    return optim_te(Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond, est, q = q, nsurr = 19, uq = uq)
end