text:
module TestModels

include("util.jl")

using QuantEcon, QuantEcon.Models
using Distributions
using Base.Test
using FactCheck

function set_up_data(cp::ConsumerProblem)
    f = get_data_file()

    # read in or compute v_star
    if exists(f, "ifp/v_star")
        v_star = read(f, "ifp/v_star")
    else
        v_init, c_init = init_values(cp)
        v_star = compute_fixed_point(x-> bellman_operator(cp, x), v_init,
                                     max_iter=5000, err_tol=1e-10,
                                     verbose=true, print_skip=10)
        write(f, "ifp/v_star", v_star)
    end

    # read in or compute c_star_vfi
    if exists(f, "ifp/c_star_vfi")
        c_star_vfi = read(f, "ifp/c_star_vfi")
    else
        c_star_vfi = get_greedy(cp, v_star)
        write(f, "ifp/c_star_vfi", c_star_vfi)
    end

    # read in or compute c_star
    if exists(f, "ifp/c_star")
        c_star = read(f, "ifp/c_star")
    else
        v_init, c_init = init_values(cp)
        c_star = compute_fixed_point(x-> coleman_operator(cp, x), c_init,
                                     max_iter=5000, err_tol=1e-10,
                                     verbose=true, print_skip=10)
        write(f, "ifp/c_star", c_star)
    end

    close(f)

    return v_star, c_star_vfi, c_star
end

function set_up_data(jv::JvWorker)
    f = get_data_file()

    # read in or compute v_star
    if exists(f, "jv/v_star")
        v_star = read(f, "jv/v_star")
    else
        v_init = [jv.x_grid .* 0.5]
        v_star = compute_fixed_point(x-> bellman_operator(jv, x), v_init,
                                     max_iter=5000, err_tol=1e-10,
                                     verbose=true, print_skip=10)
        write(f, "jv/v_star", v_star)
    end

    # read in or compute policies
    if exists(f, "jv/s_star")
        s_star = read(f, "jv/s_star")
    else
        v_init = [jv.x_grid .* 0.5]
        s_star, phi_star = get_greedy(jv, v_star)
        write(f, "jv/s_star", s_star)
    end

    if exists(f, "jv/phi_star")
        phi_star = read(f, "jv/phi_star")
    else
        if !isdefined(:phi_star)  # otherwise computed above
            _, phi_star = get_greedy(jv, v_star)
        end
        write(f, "jv/phi_star", phi_star)
    end

    close(f)

    return v_star, s_star, phi_star
end

function set_up_data(lt::LucasTree)
    f = get_data_file()

    # read in or compute prices
    if exists(f, "lt/prices")
        prices = read(f, "lt/prices")
    else
        prices = compute_lt_price(lt, err_tol=1e-10, max_iter=5000)
        write(f, "lt/prices", prices)
    end

    close(f)

    return prices
end

function set_up_data(sp::SearchProblem)
    f = get_data_file()

    # read in or compute v_star
    if exists(f, "odu/v_star")
        v_star = read(f, "odu/v_star")
    else
        v_init = fill(sp.c / (1 - sp.bet), sp.n_w, sp.n_pi)
        v_star = compute_fixed_point(x-> bellman_operator(sp, x), v_init,
                                     print_skip=1, err_tol=1e-10, max_iter=5000)
        write(f, "odu/v_star", v_star)
    end

    # read in or compute phi_vfi
    if exists(f, "odu/phi_vfi")
        phi_vfi = read(f, "odu/phi_vfi")
    else
        phi_vfi = get_greedy(sp, v_star)
        write(f, "odu/phi_vfi", phi_vfi)
    end

    # read in or compute phi_pfi
    if exists(f, "odu/phi_pfi")
        phi_pfi = read(f, "odu/phi_pfi")
    else
        phi_init = ones(sp.n_pi)
        phi_pfi = compute_fixed_point(x-> res_wage_operator(sp, x), phi_init,
                                      print_skip=1, err_tol=1e-10, max_iter=5000)
        write(f, "odu/phi_pfi", phi_pfi)
    end

    close(f)

    return v_star, phi_vfi, phi_pfi
end

function set_up_data(gm::GrowthModel)
    f = get_data_file()

    # read in or compute v_star
    if exists(f, "gm/v_star")
        v_star = read(f, "gm/v_star")
    else
        v_init = 5 .* gm.u([gm.grid]) .- 25
        v_star = compute_fixed_point(x->bellman_operator(gm, x), v_init,
                                     err_tol=1e-10, max_iter=5000, print_skip=1)
        write(f, "gm/v_star", v_star)
    end

    close(f)

    return v_star
end

facts("Testing asset_pricing.jl") do
    n = 5
    P = 0.0125 .* ones(n, n)
    P .+= diagm(0.95 .- 0.0125 .* ones(5))
    s = [1.05, 1.025, 1.0, 0.975, 0.95]
    γ = 2.0
    β = 0.94
    ζ = 1.0
    p_s = 150.0
    ap = AssetPrices(β, P, s, γ)

    tree_p = tree_price(ap)
    consol_p = consol_price(ap, ζ)
    call_option_p = call_option(ap, ζ, p_s)

    context("test shapes") do
        @fact size(ap.P, 1) --> size(ap.P, 2)
        @fact size(ap.P, 1) --> ap.n
        @fact length(tree_p) --> ap.n
        @fact length(consol_p) --> ap.n
        @fact length(call_option_p[1]) --> ap.n
    end

    context("test P_tilde") do
        P_tilde = ap.P_tilde
        P_tilde2 = similar(P_tilde)
        for i=1:ap.n, j=1:ap.n
            P_tilde2[i, j] = ap.P[i, j] * ap.s[j]^(1.0-ap.gamm)
        end
        @fact P_tilde --> roughly(P_tilde2)
    end

    context("test P_check") do
        P_check = ap.P_check
        P_check2 = similar(P_check)
        for i=1:ap.n, j=1:ap.n
            P_check2[i, j] = ap.P[i, j] * ap.s[j]^(- ap.gamm)
        end
        @fact P_check --> roughly(P_check2)
    end

    context("test multiple periods call option") do
        w_bars = call_option(ap, ζ, p_s, [5, 7])[2]
        @fact length(w_bars) --> 2
    end
end  # facts

facts("Testing career.jl") do
    cp = CareerWorkerProblem()
    v_init = rand(cp.N, cp.N)
    v_prime = bellman_operator(cp, v_init)
    greedy = get_greedy(cp, v_init)

    context("test shapes") do
        @fact size(v_init) --> size(v_prime)
        @fact size(v_init) --> size(greedy)
    end

    context("test model intuition") do
        # new life with worst job and career
        if any(greedy .== 3)
            @fact greedy[1, 1] --> 3
        end

        # new job with best career and worst job
        @fact greedy[end, 1] --> 2

        # say put with best job and career
        if any(greedy .== 1)
            @fact greedy[end, end] --> 1
        end
    end

    context("test mutating functions") do
        # new array to update inplace in bellman_operator! and get_greedy!
        out = copy(v_init)
        get_greedy!(cp, v_init, out)
        @fact out --> roughly(greedy, atol=1e-12)

        # now test bellman_operator!
        bellman_operator!(cp, v_init, out)
        @fact out --> roughly(v_prime)
    end
end  # facts

facts("Testing ifp.jl") do
    cp = ConsumerProblem()
    v_star, c_star_vfi, c_star = set_up_data(cp)

    context("test bellman and coleman policies agree") do
        @fact Base.maxabs(c_star_vfi - c_star) <= 0.2 --> true
    end

    context("test bellman solution is fixed point") do
        @fact v_star --> roughly(bellman_operator(cp, v_star); atol=1e-6)
    end

    context("test coleman solution is fixed point") do
        @fact c_star --> roughly(coleman_operator(cp, c_star); atol=1e-6)
    end

    context("test shape of init_values") do
        v_init, c_init = init_values(cp)
        shapes = (length(cp.asset_grid), length(cp.z_vals))
        @fact size(v_init) --> shapes
        @fact size(c_init) --> shapes
    end

    context("test mutating functions") do
        # new array to update inplace in bellman_operator! and get_greedy!
        out = copy(v_star)
        get_greedy!(cp, v_star, out)
        @fact out --> roughly(c_star_vfi, atol=1e-12)

        # now test bellman_operator!
        bellman_operator!(cp, v_star, out)
        @fact out --> roughly(v_star, atol=1e-12)

        # now test coleman_operator!
        coleman_operator!(cp, c_star, out)
        @fact out --> roughly(c_star; atol=1e-6)
    end
end  # facts

facts("Testing jv.jl") do
    A = 1.4
    α = 0.6
    β = 0.96
    grid_size = 50

    jv = JvWorker(A=A, alpha=α, bet=β, grid_size=grid_size)
    v_star, s_star, phi_star = set_up_data(jv)
    n = length(jv.x_grid)

    context("test model intuition") do
        # s preferred to phi with low x?
        # low x is an early index
        @fact s_star[1] > phi_star[1] --> true

        # phi preferred to s with high x?
        # high x is a late index
        @fact phi_star[end] >  s_star[end] --> true
    end

    context("policies correct size") do
        @fact length(s_star) --> n
        @fact length(phi_star) --> n
    end

    context("solution to bellman is fixed point") do
        @fact v_star --> roughly(bellman_operator(jv, v_star); atol=1e-6)
    end

    context("test mutating functions") do
        # new array to update inplace in bellman_operator! and get_greedy!
        out1 = copy(v_star)
        out2 = copy(v_star)
        get_greedy!(jv, v_star, (out1, out2))
        @fact out1 --> roughly(s_star, atol=1e-12)
        @fact out2 --> roughly(phi_star, atol=1e-12)

        # now test bellman_operator!
        bellman_operator!(jv, v_star, out1)
        @fact out1 --> roughly(v_star, atol=1e-12)
    end
end  # facts

facts("Testing lucastree.jl") do
    # model parameters
    γ = 2.0
    β  = 0.95
    α = 0.90
    σ = 0.1

    lt = LucasTree(γ, β, α, σ)
    prices = set_up_data(lt)

    @fact size(prices) --> size(lt.grid)

    context("test integrate") do
        g(x) = x .* 0.0 + 1.0
        est = Models.integrate(lt, g)

        # Same as bounds in lucastree.jl
        _int_min, _int_max = exp(-4 * σ), exp(4 * σ)
        exact = cdf(lt.phi, _int_max) - cdf(lt.phi, _int_min)

        @fact abs(est - exact) <= 0.1 --> true
    end

    context("test lucas_operator fp") do
        old_f = prices ./ (lt.grid.^γ)
        new_f = lucas_operator(lt, old_f)
        new_p = new_f .* lt.grid.^γ

        @fact prices --> roughly(new_p; atol=1e-3)
    end

    context("test prices increasing in y") do
        @fact prices --> sort(prices)
    end
end  # facts

facts("Testing odu.jl") do
    β = 0.95
    c = 0.6
    F_a = 1
    F_b = 1
    G_a = 3
    G_b = 1.2
    w_max = 2
    w_grid_size = 40
    pi_grid_size = 40

    sp = SearchProblem(β, c, F_a, F_b, G_a, G_b, w_max, w_grid_size,
                       pi_grid_size)

    v_star, phi_vfi, phi_pfi = set_up_data(sp)

    context("tests shapes of vfi outputs") do
        @fact size(v_star) --> size(phi_vfi)
    end

    context("phi_vfi increasing?") do
        phi_vfi_sorted = true
        for col=1:sp.n_pi
            if !(issorted(phi_vfi[:, col]))
                phi_vfi_sorted = false
            end
        end
        @fact phi_vfi_sorted --> true
    end

    context("v_star increasing?") do
        v_star_sorted = true
        for col=1:sp.n_pi
            if !(issorted(v_star[:, col]))
                v_star_sorted = false
            end
        end
        @fact v_star_sorted --> true
    end

    context("test model intuition") do
        # phi_pfi increasing?
        @fact issorted(phi_pfi, rev=true) --> true

        # v_star fixed_point?
        @fact v_star --> roughly(bellman_operator(sp, v_star); atol=1e-5)
    end

    context("phi_pfi fixed_point?") do
        @fact phi_pfi --> roughly(res_wage_operator(sp, phi_pfi); atol=1e-5)
    end

    context("test mutating functions") do
        # new array to update inplace in bellman_operator! and get_greedy!
        out = copy(v_star)
        get_greedy!(sp, v_star, out)
        @fact out --> roughly(phi_vfi, atol=1e-12)

        # now test bellman_operator!
        bellman_operator!(sp, v_star, out)
        @fact out --> roughly(v_star, atol=1e-12)

        # now test res_wage_operator!
        out2 = copy(phi_pfi)
        res_wage_operator!(sp, phi_pfi, out2)
        @fact out2 --> roughly(phi_pfi, atol=1e-12)
    end
end  # facts

facts("Testing optgrowth.jl") do
    # model parameters
    α = 0.65
    func(k) = k.^α
    β = 0.95
    u = log
    grid_max = 2
    grid_size = 150
    gm = GrowthModel(func, β, u, grid_max, grid_size)

    # get value/policy functions we solved for
    v_star = set_up_data(gm)
    sigma = get_greedy(gm, v_star)

    # compute true policy rule
    true_sigma = (1 - α*β) .* gm.grid.^α

    # compute true value function
    ab = α * β
    c1 = (log(1 - ab) + log(ab) * ab / (1 - ab)) / (1 - β)
    c2 = α / (1 - ab)
    true_v_star(k) = c1 + c2.*log(k)

    context("test computed/analytical policies are close") do
        @fact sigma --> roughly(true_sigma; atol=1e-2)
    end

    context("test computed/analytical values are close") do
        # Note first point is garbage b/c of interp
        @fact v_star[2:end] --> roughly(true_v_star(gm.grid)[2:end]; atol=5e-2)
    end

    context("test v_star fixed point.") do
        @fact v_star[2:end] --> roughly(bellman_operator(gm, v_star)[2:end];
                                       atol=5e-2)
    end

    context("test mutating functions") do
        # new array to update inplace in bellman_operator! and get_greedy!
        out = copy(v_star)
        get_greedy!(gm, v_star, out)
        @fact out --> roughly(sigma, atol=1e-12)

        # now test bellman_operator!
        bellman_operator!(gm, v_star, out)
        @fact out --> roughly(v_star, atol=1e-12)
    end
end  # facts

end  # module


text:
using Test

include("dfa.jl")

M = DFA(Dict(
         ("q0", 'a') => "q1",
         ("q0", 'b') => "q0",
         ("q1", 'a') => "q2",
         ("q1", 'b') => "q1",
         ("q2", 'a') => "q3",
         ("q2", 'b') => "q2",
         ("q3", 'a') => "q4",
         ("q3", 'b') => "q3",
         ("q4", 'a') => "q4",
         ("q4", 'b') => "q4"
       ), "q0", Set(["q0", "q1", "q2", "q3"]))

@test evaluate!(M, "")
@test evaluate!(M, "aaab")
@test evaluate!(M, "bbbb")
@test !evaluate!(M, "aaaab")

ext: jl

max_issues_repo_name: MadhavJivrajani/DFA-Sim

max_forks_repo_name: MadhavJivrajani/DFA-Sim


text:
using CVGSurfaceFile

files = "file\\surface.dat"
surface_struct = cvgseek(files)
save_cvgfile("file\\out.dat", surface_struct)

cvgseek2(files)
open(files)

max_stars_repo_name: powerzf/CVGSurfaceFile

max_issues_repo_name: powerzf/CVGSurfaceFile

max_forks_repo_name: powerzf/CVGSurfaceFile


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_draggableandresizable

"""
    ''_draggableandresizable(;kwargs...)
    ''_draggableandresizable(children::Any;kwargs...)
    ''_draggableandresizable(children_maker::Function;kwargs...)

A DraggableAndResizable component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `width` (Real; optional)
"""
function ''_draggableandresizable(; kwargs...)
        available_props = Symbol[:children, :id, :width]
        wild_props = Symbol[]
        return Component("''_draggableandresizable", "DraggableAndResizable", "python_version", available_props, wild_props; kwargs...)
end

''_draggableandresizable(children::Any; kwargs...) = ''_draggableandresizable(;kwargs..., children = children)
''_draggableandresizable(children_maker::Function; kwargs...) = ''_draggableandresizable(children_maker(); kwargs...)

max_stars_repo_licenses: CNRI-Python Xnet X11

max_issues_repo_licenses: CNRI-Python Xnet X11

max_forks_repo_licenses: CNRI-Python Xnet X11


text:
@testset "FockState in StateVector" begin
    dim = DIM
    T = ComplexF64

    vacuum_state_v = zeros(T, dim)
    vacuum_state_v[0 + 1] = 1

    vacuum_state = StateVector{T}(vacuum_state_v, dim)

    @test FockState(0, dim=dim) == vacuum_state
    @test NumberState(0, dim=dim) == vacuum_state
    @test VacuumState(dim=dim) == vacuum_state
    @test SinglePhotonState(dim=dim) == FockState(1, dim=dim)

    @test FockState(T, 0, dim=dim) == vacuum_state
    @test NumberState(T, 0, dim=dim) == vacuum_state
    @test VacuumState(T, dim=dim) == vacuum_state
    @test SinglePhotonState(T, dim=dim) == FockState(1, dim=dim)
end

@testset "FockState in StateMatrix" begin
    dim = DIM
    T = ComplexF64

    vacuum_state_𝛒 = zeros(T, dim, dim)
    vacuum_state_𝛒[0+1, 0+1] = 1

    vacuum_state = StateMatrix{T}(vacuum_state_𝛒, dim)

    @test FockState(0, dim=dim, rep=StateMatrix) == vacuum_state
    @test NumberState(0, dim=dim, rep=StateMatrix) == vacuum_state
    @test VacuumState(dim=dim, rep=StateMatrix) == vacuum_state
    @test SinglePhotonState(dim=dim, rep=StateMatrix) == FockState(1, dim=dim, rep=StateMatrix)

    @test FockState(T, 0, dim=dim, rep=StateMatrix) == vacuum_state
    @test NumberState(T, 0, dim=dim, rep=StateMatrix) == vacuum_state
    @test VacuumState(T, dim=dim, rep=StateMatrix) == vacuum_state
    @test SinglePhotonState(T, dim=dim, rep=StateMatrix) == FockState(1, dim=dim, rep=StateMatrix)
end


text:
@testset "Density" begin
	# Corners of the unit square
	x = [1.0 ; 0.0 ; 0.0 ; 1.0]
	y = [1.0 ; 1.0 ; 0.0 ; 0.0]

	@test density(x, y) ≈ sqrt(2)/2
end


text:
#To deal with DESI hardware sepcific functionality
"""Generate the bounding grid of focal plane
"""
struct focal_plane_Bound
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
    dx_zone::Float64
    dy_zone::Float64
    nx::Int32
    ny::Int32
end

""" Adapted from desimodel.io
https://github.com/desihub/desimodel/blob/d5f7f871873c547892a40ece30890b8610f53768/py/desimodel/io.py#L336
Loads platescale.txt.
    Returns
    -------
    :class:`~numpy.recarray`\n
        The data table read from the file.
    Notes
    -----
    The returned object has these columns:
    radius
        Radius from center of focal plane [mm].
    theta
        Radial angle that has a centroid at this radius [deg].
    radial_platescale
        Meridional (radial) plate scale [um/arcsec].
    az_platescale:
        Sagittal (azimuthal) plate scale [um/arcsec].
    arclength:
        Unknown description.

"""
function load_platescale(focalplane_dir)

    plate_scale_file="$(focalplane_dir)platescale.txt"
    plate_scale_columns=Dict("radius"=>1 ,"theta"=>2 ,"radial_platescale"=>7,
                   "az_platescale"=>8)

    rzs_file="$(focalplane_dir)rzsn.txt"
    rzs_column=Dict("R"=>1 ,"Z"=>2, "S"=>3, "N"=>4)

    #load and setup the plate scale disctionary
    platescale_dic=Dict()
    tab1=readdlm(plate_scale_file, Float64,comments=true)
    for col in keys(plate_scale_columns)
        platescale_dic[col]=tab1[:,plate_scale_columns[col]]
    end

    # One of the header column is not commented hence skipstart
    tab1=readdlm(rzs_file, Float64,comments=true,skipstart=8)

    #This function take yas first input and then x as second!!
    #https://htmlpreview.github.io/?https://github.com/PumasAI/DataInterpolations.jl/blob/v2.0.0/example/DataInterpolations.html
    rz_itp = QuadraticInterpolation(tab1[:,rzs_column["Z"]],tab1[:,rzs_column["R"]])

    platescale_dic["arclength"]=rz_itp.(platescale_dic["radius"])

    return platescale_dic
end

"""Adapted from desimodel.io
Returns an array of radii in mm given an array of radii in degrees using the platescale data
    relative to the center of the focal plane as (0,0). Supports scalar and vector inputs.\n
    Parameters
    ----------
    theta : :class:`float` or array-like
        An array that represents the angle from the center of the focal plane.
    Returns
    -------
    :class:`float` or array-like
        Radii in mm.
"""
function get_radius_mm(platescale::Dict{Any,Any},theta::Vector)

    #platescale = load_platescale(focalplane_dir)

    # Uses a quadratic one-dimensional interpolation to approximate the radius in degrees versus radius in mm
    #This function take yas first input and then x as second!!
    #https://htmlpreview.github.io/?https://github.com/PumasAI/DataInterpolations.jl/blob/v2.0.0/example/DataInterpolations.html
    rz_itp = QuadraticInterpolation(platescale["radius"],platescale["theta"])

    radius = rz_itp.(theta)

    #This to plot the interpolation for testing
    #scatter(platescale["theta"], platescale["radius"], label="input data")
    #plot!(rz_itp)
    #plot!(theta,radius) |> display

    return radius
end

function get_radius_mm(focalplane_dir::String,theta)
    platescale = load_platescale(focalplane_dir)
    radius=get_radius_mm(platescale,theta)
    return radius
end

""" Adapted from desimodel.focalplane.geometry \n
Returns arrays of the x, y positions of given celestial objects
    on the focal plane given an arbitrary telescope pointing in RA and Dec and
    arrays of the `ra` and `dec` of celestial objects in the sky.\n
    Parameters
    ----------
    telra : :class:`float`
        The telescope's RA pointing in degrees.
    teldec : :class:`float`
        The telescope's Dec pointing in degrees.
    ra : array-like
        An array of RA values for locations in the sky.
    dec : array-like
        An array of Dec values for locations in the sky.
    Returns
    -------
    tuple
        The x, y positions corrsponding to `ra`, `dec`.
    Notes
    -----
    Implements the Haversine formula.
"""
function xyz_unit_2_xy_focalplane(telra, teldec, x_unit,y_unit,z_unit,platescale::Dict{Any,Any})
    # Inclination is 90 degrees minus the declination in degrees
    nobj=size(x_unit,1)

    telra_rad=deg2rad(telra)
    teldec_rad=deg2rad(teldec)

    # Clockwise rotation around the z-axis by the right ascension of the tile center
    #Note that this uses column major order which mean
    #the first three values are in first column and not first row
    rarotate = SMatrix{3,3}(cos(telra_rad),-sin(telra_rad),0,
                            sin(telra_rad),cos(telra_rad),0,
                            0,0,1)

    # Counter-Clockwise rotation around y axis by declination of the tile center
    decrotate = SMatrix{3,3}(cos(teldec_rad),0,-sin(teldec_rad),
                            0,1,0,
                            sin(teldec_rad),0,cos(teldec_rad))

    dec_ra_rotate=decrotate*rarotate

    x_new = dec_ra_rotate[1,1].*x_unit .+ dec_ra_rotate[1,2].*y_unit .+ dec_ra_rotate[1,3].*z_unit
    y_new = dec_ra_rotate[2,1].*x_unit .+ dec_ra_rotate[2,2].*y_unit .+ dec_ra_rotate[2,3].*z_unit
    z_new = dec_ra_rotate[3,1].*x_unit .+ dec_ra_rotate[3,2].*y_unit .+ dec_ra_rotate[3,3].*z_unit

    newteldec = 0
    newtelra = 0

    ra_rad = atan.(y_new, x_new)

    dec_rad = (pi / 2) .- acos.(z_new) #assuming x,y,z on unit sphere so norm is 1

    radius_rad= (sin.(0.5 .* dec_rad).^2) .+ ((sin.(0.5 .* ra_rad).^2).*cos.(dec_rad))
    @. radius_rad = 2.0 .* asin.(sqrt.(radius_rad))

    #Finally convert this in degree but keep the variable name same for memory
    radius_rad .= radius_rad .* (180/pi)

    q_rad = atan.(z_new, -y_new)

    radius_mm = get_radius_mm(platescale,radius_rad)

    x_focalplane =@. radius_mm * cos(q_rad)
    y_focalplane =@. radius_mm * sin(q_rad)

    return x_focalplane, y_focalplane
end

function radec2xy_focalplane(telra, teldec, ra, dec,platescale::Dict{Any,Any})

    #Convert dec to inclination angle in radian
    x_unit,y_unit,z_unit = embed_sphere(ra,dec)

    x_focalplane, y_focalplane= xyz_unit_2_xy_focalplane(telra, teldec, x_unit,y_unit,z_unit,platescale)
    return x_focalplane, y_focalplane
end

"""
Estimates the focal plane boundaries for the efficient neighbour finding zone definition
"""
function get_bound_focal_plane(xfocal,yfocal,grid_scale;pad=0.2)
    xlim=[minimum(xfocal)-pad,maximum(xfocal)+pad]
    ylim=[minimum(yfocal)-pad,maximum(yfocal)+pad]

    nx=floor(Int,(xlim[2]-xlim[1])/grid_scale)
    ny=floor(Int,(ylim[2]-ylim[1])/grid_scale)

    dx_zone=(xlim[2]-xlim[1])/(nx-1)
    dy_zone=(ylim[2]-ylim[1])/(ny-1)

    fp_bound=focal_plane_Bound(xlim[1],xlim[2],ylim[1],ylim[2],dx_zone,dy_zone,nx,ny)

    return fp_bound
end

function Assign_xy_zone(xfocal,yfocal,fp_bound)
    nobj=size(xfocal,1)

    zone_x=zeros(Int32,nobj)
    zone_y=zeros(Int32,nobj)

    @. zone_x = floor(Int32,(xfocal-fp_bound.xmin)/fp_bound.dx_zone)
    @. zone_y = floor(Int32,(yfocal-fp_bound.ymin)/fp_bound.dy_zone)

    if(false)#To make some plots for visualization
        plot([0],[0],seriestype=:scatter,color=:green)
        for ii in 0:fp_bound.nx-1
            for jj in 0:fp_bound.ny-1
                indsel= findall((zone_x .== ii) .& (zone_y .== jj))
                if(size(indsel,1)>0)
                    println(ii,' ',jj,' ',size(indsel,1))
                    plot!(xfocal[indsel],yfocal[indsel],seriestype=:scatter,color=:auto)
                end
            end
        end

        plot!([0],[0],seriestype=:scatter,color=:blue) |> display
    end

    return zone_x,zone_y
end

"""
Loads the status of focal plane based on given date
look for following file with appropriate dates
    desi-exclusion_2019-09-16T00:00:00.yaml
    desi-focalplane_2019-09-16T00:00:00.ecsv
    desi-state_2019-09-16T00:00:00.ecsv
"""
function load_focal_plane_hardware(config;date="2019-09-16T00:00:00")

    harware_fp=Dict()

    file_fp="$(focalplane_dir)desi-focalplane_$(date).ecsv"
    df=CSV.read(file_fp, DataFrame;delim=" ",comment="#",header=1)

    return df

end

"""Load the focal plane file from jld2 format
if jld2 format doesn't exists then load the ascii version
and convert it to jld2 seperating into two group
one called base consisting of column we should need for fba
and other called extra contain everything else
It only loads the base part on any future call if jld2 file exists
"""
function load_hw_focalplane_file(config;date="2019-09-16T00:00:00")

    base_cols=["OFFSET_X","OFFSET_Y","LOCATION","MIN_T","MAX_T","MIN_P","MAX_P",
                    "OFFSET_T","OFFSET_P","LENGTH_R1","LENGTH_R2","DEVICE_TYPE"]

    fp_jld2="$(config["focalplane_dir_jld2"])desi-focalplane_$(date).jld2"

    if(isfile(fp_jld2))
        #fp_base=load(fp_jld2,"base")
        fp_base=read_columns_jld2(fp_jld2,base_cols)
    else
        file_fp="$(config["focalplane_dir"])desi-focalplane_$(date).ecsv"
        if(! isfile(file_fp))
            warn("File no Found $(file_fp)")
        else#Load the csv file convert to jld2 file
            df=CSV.read(file_fp, DataFrame;delim=" ",comment="#",header=1)

            DataFrame_TO_JLD2(df,fp_jld2,base_cols)
            fp_base=read_columns_jld2(fp_jld2,base_cols)
        end
    end
    return fp_base
end

function load_hw_state_file(config;date="2019-09-16T00:00:00")

    base_cols=["LOCATION", "STATE", "EXCLUSION"]
    st_jld2="$(config["focalplane_dir_jld2"])desi-state_$(date).jld2"

    if(isfile(st_jld2))
        st_base=read_columns_jld2(st_jld2,base_cols)
    else
        file_st="$(config["focalplane_dir"])desi-state_$(date).ecsv"
        if(! isfile(file_st))
            warn("File not Found $(file_st)")
        else#Load the csv file convert to jld2 file
            df=CSV.read(file_st, DataFrame;delim=" ",comment="#",header=1)
            #Write the JLD2 file
            DataFrame_TO_JLD2(df,st_jld2,base_cols)
            #st_base=load(st_jld2,"base")
            st_base=read_columns_jld2(st_jld2,base_cols)
        end
    end
    return st_base
end

"""transforms the nested ditionary loaded from yaml to flat dictionary
"""
function transform_dictionary_flatJLD(jldfile,dicin)
    jldopen(jldfile,"a+",compress=LZ4FrameCompressor()) do file
        for tkey in keys(dicin)
            for pkey in keys(dicin[tkey])
                for lkey in keys(dicin[tkey][pkey])
                    nel=length(dicin[tkey][pkey][lkey])
                    if(nel==0)
                        continue
                    elseif(nel>1)
                        @warn(" Found more than one element $(tkey) > $(pkey) > $(lkey) : $(nel)\n
                            Ignoring elements after the first one please chekc the raw file for errors")
                    end

                    tarr=dicin[tkey][pkey][lkey][1]
                    nseg=length(tarr)
                    if(lkey=="circles")
                        out_arr=[tarr[1][1],tarr[1][2],tarr[2][1]]
                    else
                        out_arr=zeros(nseg,2)
                        for ii in 1:nseg
                            out_arr[ii,:]=tarr[ii]
                        end
                    end
                    write(file,"$(tkey)-$(pkey)-$(lkey)",out_arr,compress=true)
                end
            end
        end
    end
end

"""load the flat dictionary of array and converts it into a StaticArray
"""
function load_exclusion(jldfile)
    tdic=load(jldfile)

    for tkey in keys(tdic)
        nsize=size(tdic[tkey])
        if(length(nsize)==2)
            tdic[tkey]=SMatrix{nsize[1],nsize[2],Float32}(tdic[tkey])
        else
            tdic[tkey]=SVector{nsize[1],Float32}(tdic[tkey])
        end
        #println(tkey,nsize[1])
    end
    return tdic
end

function load_hw_exclusion_file(config;date="2019-09-16T00:00:00")

    exc_jld2="$(config["focalplane_dir_jld2"])desi-exclusion_$(date).jld2"

    if(isfile(exc_jld2))
        exc=load_exclusion(exc_jld2)
    else
        file_exc="$(config["focalplane_dir"])desi-exclusion_$(date).yaml"
        #loads the yaml file
        desi_exc=YAML.load_file(file_exc)
        #flattens the ditionary, convert to array and compress and save
        transform_dictionary_flatJLD(exc_jld2,desi_exc)

        #finally load the new format and covert them to static arrays for speed
        exc=load_exclusion(exc_jld2)
    end
    return exc
end

function DataFrame_TO_JLD2(df,jldfile,base_cols;compress=false)
    #save the base property
    #tdic_base=Dict()
    #tdic_extra=Dict()

    #jldopen(jldfile, "a+",compress=LZ4FrameCompressor()) do file
    jldopen(jldfile, "a+") do file
        for col in names(df)
            #if(col in base_cols)
            if df[!,col] isa Array
                #tdic_base[col]=df[!,col]
                write(file,"$(col)",df[!,col])#,compress=false)
                #println(col,"  isa")
            else
                #tdic_extra[col]=df[!,col]
                write(file,"$(col)",df[!,col])#,compress=false)
                #println(col," not isa")
            end
        end

        #delete if needed or just write
        #write(file,"extra",tdic_extra,compress=compress)
    end

    println("saved base and extra properties in $(jldfile)")
end

function read_columns_jld2(jldfile,columns)
    tdic=Dict()
    jldopen(jldfile,"r") do file
        for col in columns
            try
               tdic[col]=read(file,col)
            catch e
                println(jldfile,col,e)
            end
        end
    end
    return tdic
end

"""
To load the full status of focal plane on certain date
first loads the focalplane, state and exclusion library
It then combines the focal plane and state together in one coherent structure
example:
config_fp=Dict("focalplane_dir" => "focalplane/",
            "focalplane_dir_jld2" => "focalplane_jld2/",
            "date" => "2019-09-16T00:00:00")

tdic=load_hw_focalplane_file(config_fp;date="2019-09-16T00:00:00")
"""
function load_hw_full_FocalPlane!(config;date="2019-09-16T00:00:00",
      default_prop=Dict("STATE"=> 0 ,"EXCLUSION"=>"default") )
    #load the focal plane
    fp_dic=load_hw_focalplane_file(config;date=date)
    #load the state
    state_dic=load_hw_state_file(config;date=date)
    #load the exclusion polygons
    exclusion_dic=load_hw_exclusion_file(config;date=date)

    #Now append a state column to the original fp_dic
    nobj=size(fp_dic["OFFSET_X"],1)
    fp_dic["STATE"]=zeros(Int8,nobj)
    fp_dic["STATE"] .= default_prop["STATE"]
    fp_dic["EXCLUSION"]=[default_prop["EXCLUSION"] for ii=1:nobj]#zeros(String,nobj)
    #In the future we might want to update other positioner info like R1,R2 or theta,phi limits

    isort_state=sortperm(state_dic["LOCATION"])
    isort_fp=sortperm(fp_dic["LOCATION"])

    nst=size(isort_state,1)
    ist=1
    for ii in 1:nobj
        if(ist>nst)#fill default
            break
        elseif(fp_dic["LOCATION"][isort_fp[ii]]==state_dic["LOCATION"][isort_state[ist]])
            for prop in ["STATE","EXCLUSION"]
               fp_dic[prop][isort_fp[ii]]= state_dic[prop][isort_state[ist]]
            end
        #elseif(fp_dic["LOCATION"][isort_fp[ii]]<state_dic["LOCATION"][isort_state[ist]])
        #    #fill the default
        #    for prop in ["STATE","EXCLUSION"]
        #       fp_dic[prop][isort_fp[ii]]= default_state[prop]
        #    end
        end
        ist +=1
    end

    #add date to the dictionaries
    fp_dic["tile_date"] =date
    exclusion_dic["tile_date"]=date

    return fp_dic,exclusion_dic
end

"""convert positioner x,y to theta phi for the positioner
"""
function Hardware_xy_to_thetaphi(pos,x,y,rad2)
    dx=x-pos.xy[1]
    dy=y-pos.xy[2]

    if(isapprox(rad2,pos.arms2_lim[2]))
        # We are at the maximum arm extension.  Force phi angle to zero
        # and compute theta.
        phi = 0.0;
        theta = atan(dy, dx)
    elseif (isapprox(rad2,pos.arms2_lim[1]))
        # We are at the limit of the arm folded inwards.  Force phi angle
        # to PI and compute theta.
        phi = pi;
        theta = atan(dy,dx)
    else #We are on neither limit.
        theta_r2=pos.arms[1]*pos.arms[1]
        phi_r2=pos.arms[2]*pos.arms[2]

        # Use law of cosines to compute "opening" angle at the "elbow".
        opening = acos((theta_r2 + phi_r2 - rad2)/ (2.0 * pos.arms[1]*pos.arms[2]))

        # The PHI angle is just the supplement of this.
        phi = pi - opening

        # Compute the theta angle.
        # Use law of cosines to compute angle from theta arm to the line from
        # the origin to the X/Y position.
        nrm_offset = sqrt(rad2)
        txy = acos((theta_r2 + rad2 - phi_r2)/ (2 * pos.arms[1] * nrm_offset))
        theta = atan(dy, dx) - txy
    end

    return theta,phi
end

function Hardware_xy_to_thetaphi(pos,x,y)
    rad2=(x-pos.xy[1])^2 + (y-pos.xy[2])^2
    return Hardware_xy_to_thetaphi(pos,x,y,rad2)
end

"""
Takes the difference of two angle
"""
function angle_diff(hi, low)
    twopi = 2.0 * pi
    # range reduction to [-Pi, Pi)
    if (hi >= pi)
        hi -= twopi
    elseif(hi < -pi)
        hi += twopi
    end

    if (low >= pi)
        low -= twopi
    elseif(low < -pi)
        low += twopi
    end

    diff = hi - low
    if (diff > 0.0)
        return diff
    else
        return twopi + diff
    end
end

"""
theta_lim consists of MIN_T,MAX_T,OFFSET_T
phi_lim consists of MIN_P, MAX_P, OFFSET_P
"""
function outside_theta_phi_range(theta,phi,theta_lim,phi_lim)
    #= Adapted from dsihub/fibreassign
    // The definitions of the theta / phi angles and ranges can be found in
    // DESI-0899.

    // Theta angle.  The theta offset loaded from desimodel is already in global
    // focalplane coordinates (not petal-local).  The theta min / max values are
    // relative to this offset (not the coordinate system).  Compute the angle
    // difference rotating both directions from theta_zero and see if we can reach
    // it at least one of those ways.  Note that theta_min is negative to indicate
    // a clockwise rotation from theta_zero.
    =#
    diff_hi = angle_diff(theta, theta_lim[3])
    diff_lo = angle_diff(theta_lim[3], theta)

    #@show "theta",diff_hi,diff_lo,theta_lim
    #=
    // Since calling code may add theta_offset (aka theta_zero here),
    // which we then subtract off, we want to put a little margin on
    // the min/max checks.  (In most cases where we are passing
    // theta/phi values through from the desimodel inputs and we need
    // to add offsets, we also ignore this range check, so this is
    // maybe overly cautious.)
    =#
    eps = 1e-12;

    if ((   ( diff_hi > (theta_lim[2] + eps)) || ( diff_hi < (theta_lim[1] - eps)))
        & ((-diff_lo > (theta_lim[2] + eps)) || (-diff_lo < (theta_lim[1] - eps)))
        )
        return true
    end

    #=
    // Phi angle.  The phi offset (phi_zero) is relative to the coordinate system
    // defined by the theta arm along the X axis.  The phi min / max values are
    // relative to this offset (not the theta-arm coordinate system).  Note that
    // a negative phi_min indicates a clockwise rotation from phi_zero.
    =#
    diff_hi = angle_diff(phi, phi_lim[3])
    diff_lo = angle_diff(phi_lim[3], phi)
    #@show "phi",diff_hi,diff_lo,phi_lim
    if ( (   ( diff_hi > (phi_lim[2] + eps)) || ( diff_hi < (phi_lim[1] - eps)) )
        & ((-diff_lo > (phi_lim[2] + eps)) || (-diff_lo < (phi_lim[1] - eps)))
        )
        return true
    end

    #// not outside range
    return false
end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
"""
Package implementing half space trees.

Half space trees are an online variant of isolation forests. They
work well when anomalies are spread out. However, they do not work
well if anomalies are packed together in windows.
"""
module HalfSpaceTrees

export HalfSpaceTree, learn!, score

import Base

include("./utils.jl")
include("./binarytree.jl")

const Feature = Dict{String,Float64};

struct Split
    on::String
    how::String
    at::Float64
end

mutable struct Data
    rmass::Int
    lmass::Int
    split::Split

    Data(rmass, lmass) = new(rmass, lmass)
    Data(rmass, lmass, split) = new(rmass, lmass, split)
end

struct WalkTreeFeature
    root::BinaryNode
    feature::Feature
end

Base.iterate(wtf::WalkTreeFeature) = (wtf.root, wtf.root)
function Base.iterate(wtf::WalkTreeFeature, state::BinaryNode)
    node = state
    if !isdefined(node.data, :split)
        return nothing
    end

    try
        node = next(node, wtf.feature)
    catch ex
        if isa(ex, KeyError)
            if node.right.data.lmass > node.left.data.lmass
                node = node.right
            else
                node = node.left
            end
        else
            rethrow(ex)
        end
    end
    return (node, node)
end

function next(node::BinaryNode, x::Feature)
    if node.data.split.how == "lt"
        if x[node.data.split.on] < node.data.split.at
            return node.left
        else
            return node.right
        end
    end
end

# HACK: Given recursion, stack overflow may happen here.
function maketree(
    limits::Dict{String,Tuple{Float64,Float64}},
    height::Int,
    padding::Float64,
    rmass::Int,
    lmass::Int)

    if height == 0
        return BinaryNode(Data(rmass, lmass))
    end

    on = rand(keys(limits))

    (a, b) = limits[on]
    at = random(a + padding * (b - a), b - padding * (b - a))

    tmp = limits[on]
    limits[on] = (tmp[1], at)
    left = maketree(limits, height - 1, padding, rmass, lmass)
    limits[on] = tmp

    tmp = limits[on]
    limits[on] = (at, tmp[2])
    right = maketree(limits, height - 1, padding, rmass, lmass)
    limits[on] = tmp

    res = BinaryNode(
        Data(rmass, lmass, Split(on, "lt", at)),
        left,
        right
    )

    return res
end

"""
    HalfSpaceTree(; <keyword arguments>)

Half space tree, an online variant of isolation forests.

# Arguments
- `ntrees::Int`: The number of trees, defaults to 10.
- `height::Int`: The height of each tree, defaults to 8.
- `windowsize::Int`: The window size, defaults to 250.
- `limits::Union{Dict{String,Tuple{Float64,Float64}},Nothing}`: Limits of each feature,
  defaults to nothing and hence to [0, 1] for each feature.

# Examples
```julia-repl
julia> hst = HalfSpaceTree(ntrees=10, height=3, windowsize=3);
```
"""
mutable struct HalfSpaceTree
    ntrees::Int
    height::Int
    windowsize::Int
    limits::Union{Dict{String,Tuple{Float64,Float64}},Nothing}
    trees::Any
    counter::Int
    firstwindow::Bool

    function HalfSpaceTree(;ntrees=10, height=8, windowsize=250, limits=nothing)
        if !(ntrees > 0)
            throw(ArgumentError("ntrees must be > 0"))
        end
        if !(height > 0)
            throw(ArgumentError("height must be > 0"))
        end
        if !(windowsize > 0)
            throw(ArgumentError("windowsize must be > 0"))
        end

        new(ntrees, height, windowsize, limits, nothing, 0, true)
    end
end

sizelimit(hst::HalfSpaceTree) = 0.1 * hst.windowsize
maxscore(hst::HalfSpaceTree) = hst.ntrees * hst.windowsize * (2^(hst.height + 1) - 1)

"""
    learn!(hst::HalfSpaceTree, x::Dict{String,Float64})

Learn feature `x`. Assumes all values of `x` to be between [0, 1]. If
not use `limits` while constucting the half space tree. Returns updated
`hst`.

# Examples
```julia-repl
julia> x = [Dict("x" => e, "y" => e, "z" => e) for e in [0.5, 0.45, 0.43, 0.44, 0.445, 0.45, 0.0]];

julia> for e in x[1:3]
           learn!(hst, e)
       end
```
"""
function learn!(hst::HalfSpaceTree, x::Feature)
    if hst.trees === nothing
        limit = (0.0, 1.0)
        limits = Dict(key => limit for key in keys(x))
        if hst.limits !== nothing
            limits = merge(limits, hst.limits)
        end

        hst.trees = [
            maketree(limits, hst.height, 0.15, 0, 0)
            for _ in 1:hst.ntrees
        ]
    end

    for tree in hst.trees
        for node in WalkTreeFeature(tree, x)
            node.data.lmass += 1
        end
    end

    hst.counter += 1
    if hst.counter == hst.windowsize
        for tree in hst.trees
            for node in inordertraversal(tree)
                node.data.rmass = node.data.lmass
                node.data.lmass = 0
            end
        end
        hst.firstwindow = false
        hst.counter = 0
    end

    return hst
end

"""
    score(hst::HalfSpaceTree, x::Dict{String,Float64})

Score feature `x`. Assumes all values of `x` to be between [0, 1]. If
not use `limits` while constucting the half space tree. Returns a score
between [0, 1] for `x`, with higher values more likely to be an anomaly.

# Examples
```julia-repl
# learn at least windowsize features
julia> score(hst, x[end - 1])
"""
function score(hst::HalfSpaceTree, x::Feature)
    if hst.firstwindow
        return 0.0
    end

    score = 0.0
    for tree in hst.trees
        for (depth, node) in enumerate(WalkTreeFeature(tree, x))
            score += node.data.rmass * 2^(depth - 1)
            if node.data.rmass < sizelimit(hst)
                break
            end
        end
    end
    score /= maxscore(hst)

    return 1 - score
end

end


text:
using Random: GLOBAL_RNG

struct ELBO <: VariationalObjective end

function (elbo::ELBO)(alg, q, logπ, num_samples; kwargs...)
    return elbo(GLOBAL_RNG, alg, q, logπ, num_samples; kwargs...)
end

function (elbo::ELBO)(
    rng::AbstractRNG,
    alg::VariationalInference,
    q,
    model::Model,
    num_samples;
    weight = 1.0,
    kwargs...
)
    return elbo(rng, alg, q, make_logjoint(model; weight = weight), num_samples; kwargs...)
end

const elbo = ELBO()


text:
"""
    PrescribedConditions{T}

Type which contains the forces, moments, displacements, and/or rotations prescribed at a point.

# Fields:
 - `isforce`: Flag for each degree of freedom indicating whether forces/moments or
    displacements are prescribed
 - `value`: Magnitude of the prescribed force/moment and/or displacement associated
    with each degree of freedom
 - `follower`: Magnitude of the follower force/moment associated with each degree of
    freedom (if applicable)
"""
struct PrescribedConditions{T}
    isforce::SVector{6, Bool}
    value::SVector{6, T}
    follower::SVector{6, T}
end
Base.eltype(::PrescribedConditions{T}) where T = T

function PrescribedConditions{T}(p::PrescribedConditions) where T
    PrescribedConditions{T}(p.isforce, p.value, p.follower)
end

PrescribedConditions(isforce, value, follower) = PrescribedConditions(isforce, promote(value, follower)...)

"""
    PrescribedConditions(t = 0.0; kwargs...)

Return the prescribed conditions at a point at time `t`.

Individual prescribed conditions may be assigned as either a scalar parameter or
as a function of time.

Prescribed Wiener-Milenkovic parameters must satisfy the following inequality:
sqrt(theta_x^2 + theta_y^2 + theta_z^2) <= 4.  Note that this restriction still
allows all possible rotations to be represented.

# Keyword Arguments
 - `ux`: Prescribed x-displacement (in the body frame)
 - `uy`: Prescribed y-displacement (in the body frame)
 - `uz`: Prescribed z-displacement (in the body frame)
 - `theta_x`: Prescribed first Wiener-Milenkovic parameter
 - `theta_y`: Prescribed second Wiener-Milenkovic parameter
 - `theta_z`: Prescribed third Wiener-Milenkovic parameter
 - `Fx`: Prescribed x-direction force
 - `Fy`: Prescribed y-direction force
 - `Fz`: Prescribed z-direction force
 - `Mx`: Prescribed x-direction moment
 - `My`: Prescribed y-direction moment
 - `Mz`: Prescribed z-direction moment
 - `Fx_follower`: Prescribed x-direction follower force
 - `Fy_follower`: Prescribed y-direction follower force
 - `Fz_follower`: Prescribed z-direction follower force
 - `Mx_follower`: Prescribed x-direction follower moment
 - `My_follower`: Prescribed y-direction follower moment
 - `Mz_follower`: Prescribed z-direction follower moment
"""
function PrescribedConditions(t = 0.0;
    ux = nothing,
    uy = nothing,
    uz = nothing,
    theta_x = nothing,
    theta_y = nothing,
    theta_z = nothing,
    Fx = nothing,
    Fy = nothing,
    Fz = nothing,
    Mx = nothing,
    My = nothing,
    Mz = nothing,
    Fx_follower = nothing,
    Fy_follower = nothing,
    Fz_follower = nothing,
    Mx_follower = nothing,
    My_follower = nothing,
    Mz_follower = nothing)

    # Set function for first slot
    if isnothing(ux)
        f1 = true
        v1 = isnothing(Fx) ? 0.0 : typeof(Fx) <: Number ? Fx : Fx(t)
        v1_f = isnothing(Fx_follower) ? 0.0 : typeof(Fx_follower) <: Number ? Fx_follower : Fx_follower(t)
    elseif !isnothing(ux) && (isnothing(Fx) && isnothing(Fx_follower))
        f1 = false
        v1 = isnothing(ux) ? 0.0 : typeof(ux) <: Number ? ux : ux(t)
        v1_f = 0.0
    else
        error("Both `ux` and `Fx` or `Fx_follower` cannot be specified at the same time")
    end

    # Set function for second slot
    if isnothing(uy)
        f2 = true
        v2 = isnothing(Fy) ? 0.0 : typeof(Fy) <: Number ? Fy : Fy(t)
        v2_f = isnothing(Fy_follower) ? 0.0 : typeof(Fy_follower) <: Number ? Fy_follower : Fy_follower(t)
    elseif !isnothing(uy) && (isnothing(Fy) && isnothing(Fy_follower))
        f2 = false
        v2 = isnothing(uy) ? 0.0 : typeof(uy) <: Number ? uy : uy(t)
        v2_f = 0.0
    else
        error("Both `uy` and `Fy` or `Fy_follower` cannot be specified at the same time")
    end

    # Set function for third slot
    if isnothing(uz)
        f3 = true
        v3 = isnothing(Fz) ? 0.0 : typeof(Fz) <: Number ? Fz : Fz(t)
        v3_f = isnothing(Fz_follower) ? 0.0 : typeof(Fz_follower) <: Number ? Fz_follower : Fz_follower(t)
    elseif !isnothing(uz) && (isnothing(Fz) && isnothing(Fz_follower))
        f3 = false
        v3 = isnothing(uz) ? 0.0 : typeof(uz) <: Number ? uz : uz(t)
        v3_f = 0.0
    else
        error("Both `uz` and `Fz` or `Fz_follower` cannot be specified at the same time")
    end

    # Set function for fourth slot
    if isnothing(theta_x)
        f4 = true
        v4 = isnothing(Mx) ? 0.0 : typeof(Mx) <: Number ? Mx : Mx(t)
        v4_f = isnothing(Mx_follower) ? 0.0 : typeof(Mx_follower) <: Number ? Mx_follower : Mx_follower(t)
    elseif !isnothing(theta_x) && (isnothing(Mx) && isnothing(Mx_follower))
        f4 = false
        v4 = isnothing(theta_x) ? 0.0 : typeof(theta_x) <: Number ? theta_x : theta_x(t)
        v4_f = 0.0
    else
        error("Both `theta_x` and `Mx` or `Mx_follower` cannot be specified at the same time")
    end

    # Set function for fifth slot
    if isnothing(theta_y)
        f5 = true
        v5 = isnothing(My) ? 0.0 : typeof(My) <: Number ? My : My(t)
        v5_f = isnothing(My_follower) ? 0.0 : typeof(My_follower) <: Number ? My_follower : My_follower(t)
    elseif !isnothing(theta_y) && (isnothing(My) && isnothing(My_follower))
        f5 = false
        v5 = isnothing(theta_y) ? 0.0 : typeof(theta_y) <: Number ? theta_y : theta_y(t)
        v5_f = 0.0
    else
        error("Both `theta_y` and `My` or `My_follower` cannot be specified at the same time")
    end

    # Set function for sixth slot
    if isnothing(theta_z)
        f6 = true
        v6 = isnothing(Mz) ? 0.0 : typeof(Mz) <: Number ? Mz : Mz(t)
        v6_f = isnothing(Mz_follower) ? 0.0 : typeof(Mz_follower) <: Number ? Mz_follower : Mz_follower(t)
    elseif !isnothing(theta_z) && (isnothing(Mz) && isnothing(Mz_follower))
        f6 = false
        v6 = isnothing(theta_z) ? 0.0 : typeof(theta_z) <: Number ? theta_z : theta_z(t)
        v6_f = 0.0
    else
        error("Both `theta_z` and `Mz` or `Mz_follower` cannot be specified at the same time")
    end

    force = SVector(f1, f2, f3, f4, f5, f6)
    value = SVector(v1, v2, v3, v4, v5, v6)
    follower = SVector(v1_f, v2_f, v3_f, v4_f, v5_f, v6_f)

    return PrescribedConditions(force, value, follower)
end

"""
    DistributedLoads{T}

Type which contains pre-integrated distributed forces and moments applied to a beam element.

# Fields
 - f1: Integrated non-follower distributed force corresponding to the start of the beam element.
 - f2: Integrated non-follower distributed force corresponding to the end of the beam element.
 - m1: Integrated non-follower distributed moment corresponding to the start of the beam element.
 - m2: Integrated non-follower distributed moment corresponding to the end of the beam element.
 - f1_follower: Integrated follower distributed force corresponding to the start of the beam element.
 - f2_follower: Integrated follower distributed force corresponding to the end of the beam element.
 - m1_follower: Integrated follower distributed moment corresponding to the start of the beam element.
 - m2_follower: Integrated follower distributed moment corresponding to the end of the beam element.
"""
struct DistributedLoads{T}
    f1::SVector{3, T}
    f2::SVector{3, T}
    m1::SVector{3, T}
    m2::SVector{3, T}
    f1_follower::SVector{3, T}
    f2_follower::SVector{3, T}
    m1_follower::SVector{3, T}
    m2_follower::SVector{3, T}
end
Base.eltype(::DistributedLoads{T}) where T = T

function DistributedLoads{T}(d::DistributedLoads) where T
    return DistributedLoads{T}(d.f1, d.f2, d.m1, d.m2, d.f1_follower, d.f2_follower,
        d.m1_follower, d.m2_follower)
end

DistributedLoads(f1, f2, m1, m2, f1_follower, f2_follower, m1_follower, m2_follower) =
    DistributedLoads(promote(f1, f2, m1, m2, f1_follower, f2_follower, m1_follower, m2_follower)...)

"""
    DistributedLoads(assembly, ielem; kwargs...)

Pre-integrate distributed loads on a beam element for use in an analysis.

# Arguments
 - `assembly`: Beam element assembly (of type [`Assembly`](@ref))
 - `ielem`: Beam element index

# Keyword Arguments
 - `s1 = 0.0`: Start of the beam element (used solely for integrating the distributed loads)
 - `s2 = 1.0`: End of the beam element (used solely for integrating the distributed loads)
 - `method = (f, s1, s2) -> gauss_quadrature(f, s1, s2)`: Method which integrates function
    `f` from `s1` to `s2`. Defaults to the Gauss-Legendre quadrature with 4 points
    on each element.
 - `fx = (s) -> 0.0`: Distributed x-direction force
 - `fy = (s) -> 0.0`: Distributed y-direction force
 - `fz = (s) -> 0.0`: Distributed z-direction force
 - `mx = (s) -> 0.0`: Distributed x-direction moment
 - `my = (s) -> 0.0`: Distributed y-direction moment
 - `mz = (s) -> 0.0`: Distributed z-direction moment
 - `fx_follower = (s) -> 0.0`: Distributed x-direction follower force
 - `fy_follower = (s) -> 0.0`: Distributed y-direction follower force
 - `fz_follower = (s) -> 0.0`: Distributed z-direction follower force
 - `mx_follower = (s) -> 0.0`: Distributed x-direction follower moment
 - `my_follower = (s) -> 0.0`: Distributed y-direction follower moment
 - `mz_follower = (s) -> 0.0`: Distributed z-direction follower moment
"""
function DistributedLoads(assembly, ielem;
    s1 = 0.0,
    s2 = 1.0,
    method = gauss_quadrature,
    fx = (s) -> 0.0,
    fy = (s) -> 0.0,
    fz = (s) -> 0.0,
    mx = (s) -> 0.0,
    my = (s) -> 0.0,
    mz = (s) -> 0.0,
    fx_follower = (s) -> 0.0,
    fy_follower = (s) -> 0.0,
    fz_follower = (s) -> 0.0,
    mx_follower = (s) -> 0.0,
    my_follower = (s) -> 0.0,
    mz_follower = (s) -> 0.0,
    )

    return DistributedLoads(assembly, ielem, 0.0;
        s1 = s1,
        s2 = s2,
        method = method,
        fx = (s, t) -> fx(s),
        fy = (s, t) -> fy(s),
        fz = (s, t) -> fz(s),
        mx = (s, t) -> mx(s),
        my = (s, t) -> my(s),
        mz = (s, t) -> mz(s),
        fx_follower = (s, t) -> fx_follower(s),
        fy_follower = (s, t) -> fy_follower(s),
        fz_follower = (s, t) -> fz_follower(s),
        mx_follower = (s, t) -> mx_follower(s),
        my_follower = (s, t) -> my_follower(s),
        mz_follower = (s, t) -> mz_follower(s),
        )
end

"""
    DistrubutedLoads(assembly, ielem, t; kwargs...)

Pre-integrate distributed loads on a beam element for use in an analysis.

# Arguments
 - `assembly`: Beam element assembly (of type [`Assembly`](@ref))
 - `ielem`: Beam element index
 - `t`: Current time

# Keyword Arguments
 - `s1 = 0.0`: Start of the beam element (used solely for integrating the distributed loads)
 - `s2 = 1.0`: End of the beam element (used solely for integrating the distributed loads)
 - `method = (f, s1, s2) -> gauss_quadrature(f, s1, s2)`: Method which integrates function
    `f` from `s1` to `s2`. Defaults to the Gauss-Legendre quadrature with 4 points on each element.
 - `fx = (s, t) -> 0.0`: Distributed x-direction force
 - `fy = (s, t) -> 0.0`: Distributed y-direction force
 - `fz = (s, t) -> 0.0`: Distributed z-direction force
 - `mx = (s, t) -> 0.0`: Distributed x-direction moment
 - `my = (s, t) -> 0.0`: Distributed y-direction moment
 - `mz = (s, t) -> 0.0`: Distributed z-direction moment
 - `fx_follower = (s, t) -> 0.0`: Distributed x-direction follower force
 - `fy_follower = (s, t) -> 0.0`: Distributed y-direction follower force
 - `fz_follower = (s, t) -> 0.0`: Distributed z-direction follower force
 - `mx_follower = (s, t) -> 0.0`: Distributed x-direction follower moment
 - `my_follower = (s, t) -> 0.0`: Distributed y-direction follower moment
 - `mz_follower = (s, t) -> 0.0`: Distributed z-direction follower moment
"""
function DistributedLoads(assembly, ielem, t;
    s1 = 0.0,
    s2 = 1.0,
    method = gauss_quadrature,
    fx = (s, t) -> 0.0,
    fy = (s, t) -> 0.0,
    fz = (s, t) -> 0.0,
    mx = (s, t) -> 0.0,
    my = (s, t) -> 0.0,
    mz = (s, t) -> 0.0,
    fx_follower = (s, t) -> 0.0,
    fy_follower = (s, t) -> 0.0,
    fz_follower = (s, t) -> 0.0,
    mx_follower = (s, t) -> 0.0,
    my_follower = (s, t) -> 0.0,
    mz_follower = (s, t) -> 0.0,
    )

    # get element length
    ΔL = assembly.elements[ielem].L

    # create function for general coordinate
    ξ = (s) -> (s-s1)/(s2-s1)

    # element length scaling
    scaling = ΔL/(s2-s1)

    # integrate to get f1 at this point in time
    f1 = SVector(
        method((s)->(1-ξ(s))*fx(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*fy(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*fz(s, t)*scaling, s1, s2)
    )

    # integrate to get f2 at this point in time
    f2 = SVector(
        method((s)->ξ(s)*fx(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*fy(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*fz(s, t)*scaling, s1, s2)
    )

    # integrate to get m1 at this point in time
    m1 = SVector(
        method((s)->(1-ξ(s))*mx(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*my(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*mz(s, t)*scaling, s1, s2)
    )

    # integrate to get m2 at this point in time
    m2 = SVector(
        method((s)->ξ(s)*mx(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*my(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*mz(s, t)*scaling, s1, s2)
    )

    # integrate to get f1_follower at this point in time
    f1_follower =SVector(
        method((s)->(1-ξ(s))*fx_follower(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*fy_follower(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*fz_follower(s, t)*scaling, s1, s2)
    )

    # integrate to get f2_follower at this point in time
    f2_follower = SVector(
        method((s)->ξ(s)*fx_follower(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*fy_follower(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*fz_follower(s, t)*scaling, s1, s2)
    )

    # integrate to get m1_follower at this point in time
    m1_follower = SVector(
        method((s)->(1-ξ(s))*mx_follower(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*my_follower(s, t)*scaling, s1, s2),
        method((s)->(1-ξ(s))*mz_follower(s, t)*scaling, s1, s2)
    )

    # integrate to get m2_follower at this point in time
    m2_follower = SVector(
        method((s)->ξ(s)*mx_follower(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*my_follower(s, t)*scaling, s1, s2),
        method((s)->ξ(s)*mz_follower(s, t)*scaling, s1, s2)
    )

    return DistributedLoads(f1, f2, m1, m2, f1_follower, f2_follower, m1_follower, m2_follower)
end

"""
    combine_loads(loads)

Combine the distributed loads in the iterable collection `loads`
"""
function combine_loads(loads)
    f1 = @SVector zeros(3)
    f2 = @SVector zeros(3)
    m1 = @SVector zeros(3)
    m2 = @SVector zeros(3)
    f1_follower = @SVector zeros(3)
    f2_follower = @SVector zeros(3)
    m1_follower = @SVector zeros(3)
    m2_follower = @SVector zeros(3)
    for load in loads
        f1 += load.f1
        f2 += load.f2
        m1 += load.m1
        m2 += load.m2
        f1_follower += load.f1_follower
        f2_follower += load.f2_follower
        m1_follower += load.m1_follower
        m2_follower += load.m2_follower
    end
    return DistributedLoads(f1, f2, m1, m2, f1_follower, f2_follower, m1_follower, m2_follower)
end

"""
    PointMass{T}

Type which contains the aggregated inertial properties of one or more point masses which
are rigidly attached to the center of an element.

# Fields:
 - `mass`: Mass matrix corresponding to the point masses.
"""
struct PointMass{T}
    mass::SMatrix{6,6,T,36}
end
Base.eltype(::PointMass{T}) where T = T

function PointMass{T}(p::PointMass) where T
    PointMass{T}(p.mass)
end

"""
    PointMass(mass)

Define a point mass given its mass matrix
"""
PointMass(mass::AbstractMatrix{T}) where T = PointMass{T}(SMatrix{6,6,T,36}(mass))

"""
    PointMass(m, p, I)

Define a point mass given its mass `m`, offset `p`, and inertia matrix `I`
"""
PointMass(m, p, I) = PointMass([m*I3 -m*tilde(p); m*tilde(p) I-m*tilde(p)*tilde(p)])

"""
    combine_masses(masses)

Combine the point masses in the iterable collection `masses`
"""
function combine_masses(masses)
    M = @SMatrix zeros(6,6)
    for mass in masses
        M += mass.M
    end
    return PointMass(M)
end

"""
    acceleration_loads(mass11, mass12, mass21, mass22, CtCab, u, a, α)

Calculate the integrated distributed loads on an element caused by acceleration.
"""
@inline function acceleration_loads(mass11, mass12, mass21, mass22, CtCab, a, α)

    # force and moment per unit length due to accelerating reference frame
    f = -CtCab*(mass11*CtCab'*a + mass12*CtCab'*α)
    m = -CtCab*(mass21*CtCab'*a + mass22*CtCab'*α)

    # integrate force and moment per unit length
    f1 = f2 = f/2
    m1 = m2 = m/2

    # return result
    return f1, f2, m1, m2
end

@inline function acceleration_loads_jacobian(mass11, mass12, mass21, mass22, a, α,
    Cab, CtCab, C_θ1, C_θ2, C_θ3, Ct_θ1, Ct_θ2, Ct_θ3)

    # force and moment per unit length due to accelerating reference frame
    f_u = -CtCab*mass11*CtCab'*tilde(α)
    m_u = -CtCab*mass21*CtCab'*tilde(α)

    f_θ = -mul3(Ct_θ1, Ct_θ2, Ct_θ3, Cab*(mass11*CtCab'*a + mass12*CtCab'*α)) -
        CtCab*mass11*Cab'*mul3(C_θ1, C_θ2, C_θ3, a) -
        CtCab*mass12*Cab'*mul3(C_θ1, C_θ2, C_θ3, α)
    m_θ = -mul3(Ct_θ1, Ct_θ2, Ct_θ3, Cab*(mass21*CtCab'*a + mass22*CtCab'*α)) -
        CtCab*mass21*Cab'*mul3(C_θ1, C_θ2, C_θ3, a) -
        CtCab*mass22*Cab'*mul3(C_θ1, C_θ2, C_θ3, α)

    # calculate integrated force and moment per unit length
    f1_u = f2_u = f_u/2
    m1_u = m2_u = m_u/2
    f1_θ = f2_θ = f_θ/2
    m1_θ = m2_θ = m_θ/2

    return f1_u, f2_u, m1_u, m2_u, f1_θ, f2_θ, m1_θ, m2_θ
end


text:
__precompile__()
module GeoArrayOps
using ProgressMeter

include("utils.jl")
include("pmf.jl")
include("smf.jl")
include("psf.jl")
include("plot.jl")
include("spread.jl")
include("terrain.jl")

export pmf, smf, psf
export pssm
export pitremoval
export spread, spread2
export roughness, TRI, TPI

# function __init__()
#     A = rand(Float64, 3, 3)
#     pssm(A)
#     spread(A, A, A)
#     spread2(A, A, A)
#     spread(A, A, 1.)
#     spread(A, 1., 1.)
#     roughness(A)
#     TRI(A)
#     TPI(A)
# end

precompile(pmf, (Matrix{Float64},))
precompile(psf, (Matrix{Float64},))
precompile(pitremoval, (Matrix{Float64},))
precompile(smf, (Matrix{Float64},))
precompile(pssm, (Matrix{Float64},))
precompile(spread, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}))
precompile(spread2, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}))
precompile(spread, (Matrix{Float64}, Matrix{Float64}, Float64))
precompile(spread, (Matrix{Float64}, Float64, Float64))
precompile(roughness, (Matrix{Float64},))
precompile(TRI, (Matrix{Float64},))
precompile(TPI, (Matrix{Float64},))

end # module