text:
"""
    met2cart(angle_met)

Convert from meteorological polar system (CW, 0 rad.=N, wind from) to cartesian polar system
(CCW, 0 rad.=E, wind to).

# Arguments
- `angle_met::Number`: an angle in radians in a meteorological coordinate system
"""
function met2cart(angle_met)
    angle_cartesian = (3.0*pi/2.0 - angle_met)

    if angle_cartesian < 0.0
        angle_cartesian += 2.0*pi
    end
    return angle_cartesian
end

"""
    rotate_to_wind_direction(xlocs, ylocs, wind_direction_met)

Rotates wind farm coordinates to be in wind direction reference where wind direction is to
the positive x.

# Arguments
- `xlocs::Array`: contains turbine east-west locations in the global reference frame
- `ylocs::Array`: contains turbine north-south locations in the global reference frame
- `wind_direction_met::Array`: contains wind direction in radians in meteorological standard
    system (N=0 rad, proceeds CW, wind from direction given)
"""
function rotate_to_wind_direction(xlocs, ylocs, wind_direction_met; center=[0.0,0.0])
    # use radians

    # convert from meteorological polar system (CW, 0 rad.=N) to standard polar system (CCW, 0 rad.=E)
    wind_direction_cart = met2cart(wind_direction_met)

    cos_wdr = cos(-wind_direction_cart)
    sin_wdr = sin(-wind_direction_cart)

    # convert to cartesian coordinates with wind to positive x
    x_cart = (xlocs.-center[1])*cos_wdr - (ylocs.-center[2])*sin_wdr
    y_cart = (xlocs.-center[1])*sin_wdr + (ylocs.-center[2])*cos_wdr

    return x_cart.+center[1], y_cart.+center[2]
end

"""
    latlong_to_xy(latitude, longitude, utm_zone; isnorth=true, units="m")

Converts arrays of points from latitude and longitude to x and y in meters in a local
coordinate frame based on the point with the lowest magnitude latitude,

# Arguments
- `latitude::Array{Float,N}`
- `longitude::Array{Float,N}`
- `isnorth::Float`: specifies if the point is in the northern hemisphere (defaul: true)
"""
function latlong_to_xy(latitude, longitude, utm_zone; isnorth=true)

    # get number of points
    npoints = length(latitude)

    # set up utm region
    utmregion = gd.UTMfromLLA(utm_zone, isnorth, wgs84)

    # get zero point index
    zp = argmin(latitude)

    # get zero point lat long
    minlla = gd.LLA(latitude[zp], longitude[zp])

    # get zero point utm
    minutm = utmregion(minlla)

    # get x and y for all points
    x = zeros(npoints)
    y = zeros(npoints)
    for i in 1:npoints
        # lat and long for current point
        lla = gd.LLA(latitude[i], longitude[i])

        # convert to utm coordinates
        utm = utmregion(lla)

        # get x location
        x[i] = utm.x - minutm.x

        # get y location
        y[i] = utm.y - minutm.y

    end

    return x, y
end

"""
    hermite_spline(x, x0, x1, y0, dy0, y1, dy1)

Produces the y and (optionally) dy values for a hermite cubic spline
interpolating between two end points with known slopes

# Arguments
- `x::Float`: x position of output y
- `x0::Float`: x position of upwind endpoint of spline
- `x1::Float`: x position of downwind endpoint of spline
- `y0::Float`: y position of upwind endpoint of spline
- `dy0::Float`: slope at upwind endpoint of spline
- `y1::Float`: y position of downwind endpoint of spline
- `dy1::Float`: slope at downwind endpoint of spline
"""
function hermite_spline(x, x0, x1, y0, dy0, y1, dy1; return_deriv=false)

    # initialize coefficients for parametric Hermite cubic spline
    c3 = (2.0*(y1))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) -
        (2.0*(y0))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) +
        (dy0)/(x0^2 - 2.0*x0*x1 + x1^2) +
        (dy1)/(x0^2 - 2.0*x0*x1 + x1^2)

    c2 = (3.0*(y0)*(x0 + x1))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) -
        ((dy1)*(2.0*x0 + x1))/(x0^2 - 2.0*x0*x1 + x1^2) - ((dy0)*(x0 +
        2.0*x1))/(x0^2 - 2.0*x0*x1 + x1^2) - (3.0*(y1)*(x0 + x1))/(x0^3 -
        3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3)

    c1 = ((dy0)*(x1^2 + 2.0*x0*x1))/(x0^2 - 2.0*x0*x1 + x1^2) + ((dy1)*(x0^2 +
        2.0*x1*x0))/(x0^2 - 2.0*x0*x1 + x1^2) - (6.0*x0*x1*(y0))/(x0^3 -
        3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) + (6.0*x0*x1*(y1))/(x0^3 -
        3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3)

    c0 = ((y0)*(- x1^3 + 3.0*x0*x1^2))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 -
        x1^3) - ((y1)*(- x0^3 + 3.0*x1*x0^2))/(x0^3 - 3.0*x0^2*x1 +
        3.0*x0*x1^2 - x1^3) - (x0*x1^2*(dy0))/(x0^2 - 2.0*x0*x1 + x1^2) -
        (x0^2*x1*(dy1))/(x0^2 - 2.0*x0*x1 + x1^2)

    # Solve for y and dy values at the given point
    y = c3*x^3 + c2*x^2 + c1*x + c0
    dy_dx = c3*3*x^2 + c2*2*x + c1

    if return_deriv
        return y, dy_dx
    else
        return y
    end
end

"""
    overlap_area_func(turbine_y, turbine_z, rotor_diameter, wake_center_y,
    wake_center_z, wake_diameter; tol=1E-6)

Produces the y and (optionally) dy values for a hermite cubic spline
interpolating between two end points with known slopes

All dimensions should be in meters

# Arguments
- `turbine_y::Float`: cross wind location of turbine hub
- `turbine_z::Float`: vertical location of turbine hub
- `rotor_diameter::Float`
- `wake_center_y::Float`: cross wind location of wake center
- `wake_center_z::Float`: vertical location of wake center
- `wake_diameter::Float`
- `tol::Float`: extra distance for comparisons. Default 1E-6
"""
# calculates the overlap area between a given wake and a rotor area
function overlap_area_func(turbine_y, turbine_z, rotor_diameter, wake_center_y,
    wake_center_z, wake_diameter; tol=1E-6)

    # distance between wake center and rotor center
    if (wake_center_z > (turbine_z + tol)) || (wake_center_z < (turbine_z - tol))
        OVdYd = norm([wake_center_y-turbine_y, wake_center_z - turbine_z])
    elseif (wake_center_y > (turbine_y + tol))
        OVdYd = wake_center_y - turbine_y
    elseif (turbine_y > (wake_center_y + tol))
        OVdYd = turbine_y - wake_center_y
    else
        OVdYd = 0.0
    end

    # find rotor radius
    OVr = rotor_diameter/2.0

    # find wake radius
    OVRR = wake_diameter/2.0

    # make sure the distance from wake center to turbine hub is positive
    OVdYd = abs_smooth(OVdYd, eps())

    # determine if there is overlap
    if (OVdYd < (OVr+OVRR)) # if the rotor overlaps the wake zone

        # check that turbine and wake centers are not perfectly aligned
        if (OVdYd > (0.0 + tol))

            # check if the rotor is wholly contained in the wake
            if ((OVdYd + OVr) < OVRR + tol)
                # wake_overlap = pi*OVr*OVr
                wake_overlap = 3.1415926535897*OVr*OVr
            elseif ((OVdYd + OVRR) < OVr + tol)
                # wake_overlap = pi*OVRR*OVRR
                wake_overlap = 3.1415926535897*OVRR*OVRR
            else
                # calculate the distance from the wake center to the chord connecting the lens cusps
                OVL = (-OVr*OVr+OVRR*OVRR+OVdYd*OVdYd)/(2.0*OVdYd)

                OVz = sqrt(OVRR*OVRR-OVL*OVL)
                OVz2 = sqrt(OVr*OVr-(OVdYd-OVL)*(OVdYd-OVL))

                wake_overlap = OVRR*OVRR*acos(OVL/OVRR) + OVr*OVr*acos((OVdYd-OVL)/OVr) - OVL*OVz - (OVdYd-OVL)*OVz2
            end

        # perfect overlap case where the wake is larger than the rotor
        elseif (OVRR > OVr)
            wake_overlap = pi*OVr*OVr
        # perfect overlap case where the rotor is larger than the wake
        else
            wake_overlap = pi*OVRR*OVRR
        end

    # case with no overlap
    else
        wake_overlap = 0.0
    end

end

"""
    smooth_max_ndim(x; s=100.0)

Calculate the smooth-max (a.k.a. softmax or LogSumExponential) of the elements in x.

Based on John D. Cook's writings at
(1) https://www.johndcook.com/blog/2010/01/13/soft-maximum/
and
(2) https://www.johndcook.com/blog/2010/01/20/how-to-compute-the-soft-maximum/

# Arguments
- `x::Float`: first value for comparison
- `y::Float`: second value for comparison
- `s::Float` : controls the level of smoothing used in the smooth max
"""
function smooth_max(x, y; s=10.0)

    # LogSumExponential Method - used this in the past
    # g = (x*exp(s*x)+y*exp(s*y))/(exp(s*x)+exp(s*y))

    # non-overflowing version of Smooth Max function (see ref 2 above)
    max_val = max(x, y)
    min_val = min(x, y)
    r = (log(1.0 + exp(s*(min_val - max_val))) + s*max_val)/s

    return r

end

"""
    smooth_max(x; s=10.0)

Calculate the smoothmax (a.k.a. softmax or LogSumExponential) of the elements in x.

Based on John D. Cook's writings at
(1) https://www.johndcook.com/blog/2010/01/13/soft-maximum/
and
(2) https://www.johndcook.com/blog/2010/01/20/how-to-compute-the-soft-maximum/

And based on article in FeedlyBlog
(3) https://blog.feedly.com/tricks-of-the-trade-logsumexp/

# Arguments
- `x::Array{Float,1}` : vector with all the input values
- `s::Float` : controls the level of smoothing used in the smooth max
"""
function smooth_max(x; s=10.0)

    # non-overflowing version of Smooth Max function (see ref 2 and 3 above)

    # get the maximum value and the index of maximum value
    max_val, max_ind = findmax(real(x))

    # get the indices of x
    indices = collect(1:length(x))

    # remove the index of the maximum value
    splice!(indices, max_ind)

    # LogSumExp with smoothing factor s
    r = (log(sum([1.0; exp.(s*(x[indices].-max_val))])) + s*max_val)/s

    return r

end

"""
    closeBndryLists(bndryPts_x, bndryPts_y)

Appends the 1st element to the end of each array for a closed boundary.
Note, this will not function properly if there is only one region.
For only one region, use `closeBndryList(bndryPts_x, bndryPts_y)` (note the
singular, not plural 'List' in the function title)

# Arguments
- `bndryPts_x::Array{Float,1}` : N-D array of x-coordinates for the vertices
        around N-many closed boundaries
- `bndryPts_y::Array{Float,1}` : N-D array of y-coordinates for the vertices
        around N-many closed boundaries
"""
function closeBndryLists(region_bndry_x, region_bndry_y)
    # Determine how many regions and points per region were passed
    nRegions = length(region_bndry_x)
    # For every region
    for i in 1:nRegions
        # Append the initial points to the end of that row (if needed)
        region_bndry_x[i], region_bndry_y[i] = closeBndryList(region_bndry_x[i], region_bndry_y[i])
    end

    return region_bndry_x, region_bndry_y
end

"""
    closeBndryList(bndryPts_x, bndryPts_y)

Appends the 1st element to the end of the coordinate arrays if it is not already
repeated. Note, this will only work on 1-D arrays. For an array of 1-D arrays,
use `closeBndryLists(bndryPts_x, bndryPts_y)` (note the
plural, not singular 'Lists' in the function title)

# Arguments
- `bndryPts_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndryPts_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
"""
function closeBndryList(bndry_x, bndry_y)
    # If the end point isn't already a repeat of the staring point
    if !((bndry_x[1] == bndry_x[end]) && (bndry_y[1] == bndry_y[end]))
        # Append the first point to the end of the array
        bndryPts_x_clsd = push!(bndry_x, bndry_x[1])
        bndryPts_y_clsd = push!(bndry_y, bndry_y[1])
    end

    return bndryPts_x_clsd, bndryPts_y_clsd
end

"""
    PointsOnCircum(center_x, center_y, r, n = 100)

Given a circle center, radius, and number of discrete points, returns an array
of discrete points along the circle's circumference

# Arguments
- `center_x::Float64` : cartesian x-coordinate for the center of the circle
- `center_y::Float64` : cartesian y-coordinate for the center of the circle
- `r::Float64` : distance from circle's center to the circumference points
- `n::Float64` : defaults to 100, is the number of discrete evenly-spaced points
        that will be returned along the circle's circumference
"""
function DiscreteCircum(center_x, center_y, r, n = 100)
    bndry_x = zeros(n+1)
    bndry_y = zeros(n+1)
    for i in 0:n
        bndry_x[i+1] = center_x + cos(2*pi/n*i)*r
        bndry_y[i+1] = center_y + sin(2*pi/n*i)*r
    end
    return bndry_x, bndry_y
end

"""
    calcMinorAngle(bndry_x, bndry_y, bndry_z=[0,0,0])

Given three points in space, calculates the magnitude of the non-reflex angle
formed at the center point. Created to be used in VR_bounary_startup()

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
- `bndry_z::Array{Float,1}` : 1-D array of z-coordinates for the vertices
        around a singlar closed boundary. Default to [0,0,0] for X-Y plane
"""
function calcMinorAngle(bndry_x, bndry_y, bndry_z=[0,0,0])
    # Calculates the magnitude of the non-reflex angle formed at the center point
    ABx = bndry_x[2]-bndry_x[1]
    ABy = bndry_y[2]-bndry_y[1]
    ABz = bndry_z[2]-bndry_z[1]

    BCx = bndry_x[2]-bndry_x[3]
    BCy = bndry_y[2]-bndry_y[3]
    BCz = bndry_z[2]-bndry_z[3]

    Num = (ABx*BCx) + (ABy*BCy) + (ABz*BCz) # Dot Product

    Denom = norm([ABx, ABy, ABz]) * norm([BCx, BCy, BCz]) # Multiplication of magnitudes
    Theta = acosd(Num/Denom) # Get the angle formed

    # If it's greater than 180, get the
    if (Theta > 180.0)
        Theta = 360.0 - Theta    # Get the explementary angle
    end

    return Theta
end

"""
    calcSmallestAngle(bndry_x_clsd, bndry_y_clsd)

Given a 1-D closed array of boundary verticies (with first point repeated at the
end) it determines the smallest non-reflex angle created by any three
consecutive verticies along the boundary. Created to be used in
VR_bounary_startup()

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
"""
function calcSmallestAngle(bndry_x_clsd, bndry_y_clsd)
    num_angles = length(bndry_x_clsd)-1
    # Loop the second point back on so we get the angle at the beginning
    if !(((bndry_x_clsd[1] == bndry_x_clsd[end-1]) && (bndry_y_clsd[1] == bndry_y_clsd[end-1])) && ((bndry_x_clsd[2] == bndry_x_clsd[end]) && (bndry_y_clsd[2] == bndry_y_clsd[end])))
        bndryPts_x_loopd = vcat(bndry_x_clsd, bndry_x_clsd[2])
        bndryPts_y_loopd = vcat(bndry_y_clsd, bndry_y_clsd[2])
    end

    #- Calculate the smallest angle -#
    smallest_angle = 360
    for i in 1:num_angles
        temp_angle = calcMinorAngle(bndryPts_x_loopd[i:i+2], bndryPts_y_loopd[i:i+2])

        if (temp_angle < smallest_angle)
            smallest_angle = temp_angle
        end
    end

    return smallest_angle
end

"""
    getPerimeterLength(bndry_x_clsd, bndry_y_clsd)

Given a 1-D closed array of boundary verticies (with first point repeated at the
end) returns the length along the perimeter. Created to be used in
VR_bounary_startup()

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
"""
function getPerimeterLength(bndry_x_clsd, bndry_y_clsd)
    num_bndry_pts = length(bndry_x_clsd)-1
    nLength = zeros(num_bndry_pts)
    for i in 1:num_bndry_pts
        nLength[i] = norm([bndry_x_clsd[i+1]-bndry_x_clsd[i], bndry_y_clsd[i+1]-bndry_y_clsd[i]])
    end

    return nLength
end

"""
    coordDist(x1, y1, x2, y2)

Given a two points (x1, y1) and (x2, y2), returns the euclidean distance between
them

# Arguments
- `x1::Float64` : x-coord of the first point
- `y1::Float64` : y-coord of the first point
- `x2::Float64` : x-coord of the second point
- `y2::Float64` : y-coord of the second point
"""
function coordDist(x1, y1, x2, y2)
    xDiff = x1 - x2
    yDiff = y1 - y2
    return norm([xDiff, yDiff])
end

"""

    single_boundary_normals_calculator(boundary_vertices)

Outputs the unit vectors perpendicular to each edge of a polygon, given the Cartesian
coordinates for the polygon's vertices.

# Arguments
- `boundary_vertices::Array{Float,1}` : m-by-2 array containing all the boundary vertices, counterclockwise

"""

function single_boundary_normals_calculator(boundary_vertices)

    # get number of vertices in shape
    nvertices = length(boundary_vertices[:, 1])

    # add the first vertex to the end of the array to form a closed-loop
    boundary_vertices = [boundary_vertices; boundary_vertices[1,1] boundary_vertices[1,2]]

    # initialize array to hold boundary normals
    boundary_normals = zeros(nvertices, 2)

    # iterate over each boundary
    for i = 1:nvertices

        # create a vector normal to the boundary
        boundary_normals[i, :] = [ -(boundary_vertices[i+1, 2] - boundary_vertices[i, 2]) ; boundary_vertices[i+1, 1] - boundary_vertices[i, 1] ]

        # normalize the vector
        boundary_normals[i, :] = boundary_normals[i, :] / norm(boundary_normals[i, :])

    end

    return boundary_normals

end

"""

    boundary_normals_calculator(boundary_vertices; nboundaries=1)

Outputs the unit vectors perpendicular to each edge of each polygon in a set of polygons,
given the Cartesian coordinates for each polygon's vertices.

# Arguments
- `boundary_vertices::Array{Float,1}` : ragged array of arrays containing all the boundary vertices of each polygon, counterclockwise
- `nboundaries::Int` : the number of boundaries in the set
"""
function boundary_normals_calculator(boundary_vertices; nboundaries=1)

    if nboundaries == 1
        boundary_normals = single_boundary_normals_calculator(boundary_vertices)
    else
        boundary_normals = deepcopy(boundary_vertices)
        for i = 1:nboundaries
            boundary_normals[i] = single_boundary_normals_calculator(boundary_vertices[i])
        end
    end
    return boundary_normals
end

"""

    _remove_perimeter_points!(n; alpha=0.0)

Internal function. Removes points outside or outside and on the border of the rotor-swept
    area

# Arguments
- `y::AbstractArray`: horizontal point locations
- `z::AbstractArray`: vertical point locations
- `use_perimeter_points::Bool`: flag that determines whether or not to include points on the
    boundary of the rotor-swept area
"""
function _remove_out_of_bounds_points(y, z, use_perimeter_points)

    # get x and y values separately and crop to only include points in the swept area
    if use_perimeter_points
        yy = y[hypot.(y,z) .<= 1]
        zz = z[hypot.(y,z) .<= 1]
    else
        yy = y[hypot.(y,z) .< 1]
        zz = z[hypot.(y,z) .< 1]
    end

    # return new arrays with points outside boundary removed
    return yy, zz

end

"""

    sunflower_points(n; alpha=0.0)

Generates points in a circle of radius=1 using the sunflower packing algorithm.

# Arguments
- `n::Float`: number of points to generate
- `alpha::Float`: Controls the smoothness of the boundary. alpha=0 is the standard "jagged edge" sunflower algoirthm and
    alpha=1 results in a smooth boundary.
"""
function sunflower_points(n; alpha=0.0)
    # this function generates n points within a circle in a sunflower seed pattern
    # the code is based on the example found at
    # https://stackoverflow.com/questions/28567166/uniformly-distribute-x-points-inside-a-circle

    function radius(k, n, b)
        if k > n - b
            r = 1 # put on the boundary
        else
            r = sqrt(k - 1.0/2.0)/sqrt(n - (b + 1.0)/2.0)  # apply squareroot
        end
        return r
    end

    x = zeros(n)
    y = zeros(n)

    b = round(alpha*sqrt(n)) # number of boundary points

    phi = (sqrt(5.0) + 1.0)/2.0  # golden ratio

    for k in 1:n
        r = radius(k, n, b)

        theta = 2.0*pi*k/phi^2

        x[k] = r*cos(theta)
        y[k] = r*sin(theta)
    end

    return x, y
end

"""

    grid_points(n)

Generates points in a grid. If n is not a perfect square, then the nearest square root will
be used for the side length of the grid.

# Arguments
- `n::Float`: number of points to generate
"""
function grid_points(n)
    # this function generates n points within a circle using a grid pattern
    # the code is based on NREL's floris approach

    # determine length of x and y in grid, round if not perfect square
    sidepoints = Int(round(sqrt(n), digits=0))

    # generate horizontal points
    y = -1.0:2.0/(sidepoints-1.0):1.0

    # generate vertical points
    z = -1.0:2.0/(sidepoints-1.0):1.0

    # generate grid
    grid = repeat(y',length(z),1),repeat(z,1,length(y))

    # get new number of points
    nnew = length(grid[1])

    # return 1D arrays for x and y
    return reshape(grid[1],nnew), reshape(grid[2], nnew)
end

"""
    rotor_sample_points(nsamplepoints=1)

Initializes the sampling locations in the rotor-swept-area. Returns values such that
zero is at the turbine hub location and 1 is at the tip of the blades. If a single
sample is requested, it will be at the hub location. Otherwise, the points will be located
using the sunflower packcing algorithm.

# Arguments
- `nsamplepoints::Int`: controlls how many sample points to generate
- `alpha::Float`: Controls smoothness of the sunflower algorithm boundary. alpha=0 is the standard "jagged edge" sunflower algoirthm and
    alpha=1 results in a smooth boundary.
- `pradius::Float`: the percent of the rotor radius to use in generating initial point grid
- `use_perimeter_points`: whether or not to include point exactly on the perimeter of the
    rotor swept area
"""
function rotor_sample_points(nsamplepoints=1; method="sunflower", alpha=0.0, use_perimeter_points=true, pradius=1.0)

    if nsamplepoints > 1
        if method == "sunflower"
            rotor_sample_points_y, rotor_sample_points_z = ff.sunflower_points(nsamplepoints, alpha=alpha)
        elseif method == "grid"
            rotor_sample_points_y, rotor_sample_points_z = ff.grid_points(nsamplepoints)
        end

        # adjust to desired radius
        rotor_sample_points_y .*= pradius
        rotor_sample_points_z .*= pradius

        # remove any points outside the swept area
        rotor_sample_points_y, rotor_sample_points_z = _remove_out_of_bounds_points(rotor_sample_points_y, rotor_sample_points_z, use_perimeter_points)

    else
        rotor_sample_points_y = rotor_sample_points_z = [0.0]
    end

    return rotor_sample_points_y, rotor_sample_points_z
end

"""
    print_state_layouts_in_cartesian_frame(turbinex, turbiney, winddirections)

Given a wind farm layout in the global reference frame, print the layout rotated to the
cartesian frame with wind to the positive x axis (right) for all wind directions.

# Arguments
- `turbinex::Array{T,1}`: x locations of turbines in global reference frame
- `turbiney::Array{T,1}`: y locations of turbines in global reference frame
- `winddirections::Array{T,1}`: all wind directions in radians in meteorological coordinates (0 rad. = from North)
"""
function print_layout_in_cartesian_frame_excel(turbinex, turbiney, rotordiameter, winddirections, outfile; center=[0.0,0.0], plot_layouts=false, round_locations=false)
    # get number of directions
    ndirections = length(winddirections)

    # initialize excel file
    XLSX.openxlsx(outfile, mode="w") do xf

        sheet = xf[1]
        sheet["A1"] = "wind direction"
        sheet["B1"] = "turbine x"
        sheet["C1"] = "turbine y"

        # loop through directions
        for i in 1:ndirections

            # rotate turbinex and turbiney to current direction
            rotx, roty = rotate_to_wind_direction(turbinex, turbiney, winddirections[i], center=center)

            # round if desired
            if round_locations
                rotx = round.(rotx)
                roty = round.(roty)
            end

            # plot layouts if desired, pause for 5 seconds on each
            if plot_layouts
                fig, ax = plt.subplots()
                ax.set(aspect="equal", xlim=[minimum(turbinex)-200,maximum(turbinex)+200], ylim=[minimum(turbiney)-200,maximum(turbiney)+200])
                ff.plotlayout!(ax, rotx, roty, rotordiameter)
                ff.plotlayout!(ax, [rotx[1], rotx[20]], [roty[1], roty[20]], rotordiameter, fill=true, color="r")
                plt.title("$(winddirections[i]*180.0/pi)")
                plt.show()

                # sleep(1)

            end

            XLSX.rename!(sheet, "rotated layouts")

            # print rotated coordinates to specified output
            sheet["A$(i+1)"] = winddirections[i]*180.0/pi
            sheet["B$(i+1)"] = join(rotx, ",")
            sheet["C$(i+1)"] = join(roty, ",")

        end
    end
end

function wrap_180(x)
    """
    Wrap an angle to between -180 and 180
    adapted from NREL's floris
    """

    x[x .<= -180.0] .+= 360.0
    x[x .> 180.0] .-= 360.0

    return(x)
end

# the following is from floris for calcculating wake counts
"""
    wake_count_iec(turbinex, turbiney, winddirection, diameter; return_turbines=true)

    Adapted from NREL's floris

    Finds the number of turbines waking each turbine for the given
    wind direction. Waked directions are determined using the formula
    in Figure A.1 in Annex A of the IEC 61400-12-1:2017 standard.

# Arguments
- `turbinex::Array{T,1}`: x locations of turbines in global reference frame
- `turbiney::Array{T,1}`: y locations of turbines in global reference frame
- `winddirection::Float`: wind direction in radians in meteorological coordinates (0 rad. = from North)
- `diameter::Array{T,1}`: diameters of all wind turbines
"""
function wake_count_iec(turbinex, turbiney, wd::Real, diameter; return_turbines=true)

    # convert wind direction to degrees
    wd *= 180.0/pi

    # get number of turbines
    nturbines = length(turbinex)

    # set indices of all turbines
    turbines = collect(1:nturbines)

    # initialize wake count list
    wake_list = []

    # loop through turbines
    for turbi = 1:nturbines

        # get list of all other turbines
        other_turbines = turbines[turbines .!= turbi]

        # calculate distance in diameters from turbi to all other turbines
        dists = hypot.(turbinex[other_turbines] .- turbinex[turbi], turbiney[other_turbines] .- turbiney[turbi])./diameter[other_turbines]

        # calculate angles in degrees from other turbines to turbi
        angles = rad2deg.(atan.(turbinex[other_turbines] .- turbinex[turbi], turbiney[other_turbines] .- turbiney[turbi]))

        waked = dists .<= 2.0
        waked = waked .| (&).(
            (dists .<= 20.0),
            (
                abs.(wrap_180(wd .- angles)) .<= 0.5 .* (1.3 .* rad2deg.(atan.(2.5 ./ dists .+ 0.15)) .+ 10)
            )
        )

        push!(wake_list, sum(waked))

    end

    return wake_list
end

function wake_count_iec(turbinex, turbiney, wd::AbstractArray, diameter; return_turbines=true)

    # initialize wake list
    wake_list = zeros(Int, (length(wd), length(turbinex)))

    # call wake_count_iec for each direction
    for i = 1:length(wd)
        wake_list[i,:] = wake_count_iec(turbinex, turbiney, wd[i], diameter)
    end

    return wake_list
end

"""
    find_upstream_turbines(turbinex, turbiney, winddirection, diameter; inverse=false)

A convenience function to quickly find either which turbines are waked, or those that are
not.

# Arguments
- `turbinex::Array{T,1}`: x locations of turbines in global reference frame
- `turbiney::Array{T,1}`: y locations of turbines in global reference frame
- `winddirection::Real` or `winddirection::AbstractArray`: wind direction in radians in meteorological coordinates (0 rad. = from North)
- `diameter::Array{T,1}`: diameters of all wind turbines
"""
function find_upstream_turbines(turbinex, turbiney, winddirection::AbstractArray, diameter; inverse=false)

    # find wake count for all turbines in given wind direction
    wake_count = []

    for wd in winddirection
        push!(wake_count, wake_count_iec(turbinex, turbiney, wd, diameter))
    end

    if inverse
        # return waked turbines
        returnarray = []
        for i = 1:length(winddirection)
            push!(returnarray, collect(1:length(turbinex))[wake_count[i] .!= 0])
        end
        return returnarray
    else
        # return unwaked turbines
        returnarray = []
        for i = 1:length(winddirection)
            push!(returnarray, collect(1:length(turbinex))[wake_count[i] .== 0])
        end
        return returnarray
    end

end

function find_upstream_turbines(turbinex, turbiney, winddirection::Real, diameter; inverse=false)

    # find wake count for all turbines in given wind direction
    wake_count = wake_count_iec(turbinex, turbiney, winddirection, diameter)

    if inverse
        # return waked turbines
        return collect(1:length(turbinex))[wake_count .!= 0]
    else
        # return unwaked turbines
        return collect(1:length(turbinex))[wake_count .== 0]
    end

end

function round_farm_concentric_start(rotor_diameter, center, radius; min_spacing=2.)

    # normalize inputs
    radius /= rotor_diameter
    center /= rotor_diameter

    # calculate how many circles can be fit in the wind farm area
    nCircles = floor(radius/min_spacing)
    radii = range((radius/nCircles), radius, length = Int(nCircles))
    alpha_mins = 2.0*asin.(min_spacing./(2.0.*radii))
    nTurbines_circles = floor.(2.0 * pi ./ alpha_mins)
    # println(radius)
    # println(radii)
    nTurbines = Int(sum(nTurbines_circles))+1
    # println(nCircles)
    # println(nTurbines_circles)

    alphas = 2.0*pi./nTurbines_circles

    turbineX = zeros(nTurbines)
    turbineY = zeros(nTurbines)

    index = 1
    turbineX[index] = center[1]
    turbineY[index] = center[2]
    index += 1
    for circle in 1:Int(nCircles)
        for turb in 1:Int(nTurbines_circles[circle])
            angle = alphas[circle]*(turb-1)
            w = radii[circle]*cos(angle)
            h = radii[circle]*sin(angle)
            x = center[1] + w
            y = center[2] + h
            turbineX[index] = x
            turbineY[index] = y
            index += 1
        end
    end

    return turbineX*rotor_diameter, turbineY*rotor_diameter
end

"""
    round_farm_random_start(rotor_diameter, center, radius; min_spacing=2., min_spacing_random=3., method="individual")

Generates starting locations for multi-start optimization approaches when the farm boundary is round.

# Arguments
- `rotor_diameter::Number`: wind turbine diameter
- `center::Number`: wind farm center
- `radius::Number`: wind farm radius
- `diameter::Array{T,1}`: diameters of all wind turbines
"""
function round_farm_random_start(rotor_diameter, center, radius; nturbines=nothing, min_spacing=2., min_spacing_random=3., method="individual", alpha=0.0)
    # normalize inputs
    radius /= rotor_diameter
    center /= rotor_diameter

    if nturbines === nothing
        # calculate how many circles can be fit in the wind farm area
        ncircles = floor(radius / min_spacing)
        radii = range((radius/ncircles), radius, length = Int(ncircles))
        alpha_mins = 2.0 * asin.(min_spacing ./ (2.0 * radii))
        nturbines_circles = floor.(2.0 * pi ./ alpha_mins)

        nturbines = Int(sum(nturbines_circles)) + 1
    end

    turbinex = zeros(nturbines)
    turbiney = zeros(nturbines)

    if method == "individual"
        # generate random points within the wind farm boundary
        for i = 1:nturbines

            good_point = false

            while !good_point

                # generate random point in containing rectangle
                turbinex[i], turbiney[i] = rand(1,2).*2.0 .- 1.0

                turbinex[i] *= radius
                turbiney[i] *= radius

                turbinex[i] += center[1]
                turbiney[i] += center[2]

                # calculate signed distance from the point to each boundary facet
                distance = radius - sqrt((turbinex[i] - center[1])^2 + (turbiney[i] - center[2])^2)

                # determine if the point is inside the wind farm boundary
                if distance > 0.0
                    n_bad_spacings = 0
                    for turb = 1:nturbines
                        if turb >= i
                            continue
                        end
                        spacing = sqrt((turbinex[turb] - turbinex[i])^2 + (turbiney[turb] - turbiney[i])^2)
                        if spacing < min_spacing_random
                            n_bad_spacings += 1
                        end
                    end
                    if n_bad_spacings == 0
                        good_point = true
                    end
                end
            end
            print(".")
        end
    elseif method == "angle"
        turbinex, turbiney = round_farm_concentric_start(copy(rotor_diameter), copy(center*rotor_diameter), copy(radius*rotor_diameter), min_spacing=min_spacing_random)

        turbinex /= rotor_diameter
        turbiney /= rotor_diameter

        turbinex .-= center[1]
        turbiney .-= center[2]

        # get rotation angle
        step = 0.001
        rotationangle = rand(0:step:(2*pi-step)) # in radians

        # rotate
        turbinex[:], turbiney[:] = rotate_to_wind_direction(turbinex, turbiney, rotationangle)

        # translate
        turbinex .+= center[1]
        turbiney .+= center[2]

    elseif method == "angle-each-circle"

        turbinex, turbiney = round_farm_concentric_start(copy(rotor_diameter), copy(center*rotor_diameter), copy(radius*rotor_diameter), min_spacing=min_spacing_random)

        turbinex /= rotor_diameter
        turbiney /= rotor_diameter

        turbinex .-= center[1]
        turbiney .-= center[2]

        # get rotation angle
        step = 0.001
        circleidx = [2:7, 8:19, 20:38]
        for i in 1:3
            rotationangle = rand(0:step:(2*pi-step)) # in radians

            # rotate circle i
            turbinex[circleidx[i]], turbiney[circleidx[i]] = rotate_to_wind_direction(turbinex[circleidx[i]], turbiney[circleidx[i]], rotationangle)
        end
        # translate
        turbinex .+= center[1]
        turbiney .+= center[2]

    elseif method == "concentric"
        # calculate how many circles can be fit in the wind farm area
        maxcircles = floor(radius / min_spacing_random)

        # get max number of turbines that will fit on the boundary
        alpha_min = 2.0 * asin.(min_spacing_random / (2.0 * radius))
        max_boundary_turbines = floor.(2.0 * pi / alpha_min)

        if max_boundary_turbines >= nturbines - 1
            mincircles = 1
        else
            mincircles = 2
        end
        mincircles = 3

        # choose how many circles (random)
        ncircles = Int(rand(mincircles:maxcircles))

        # initialize circles
        radii = range((radius/ncircles), radius, length = Int(ncircles))

        remaining_turbines = nturbines - 1
        circle_turbines = zeros(ncircles)
        for i = 1:length(radii)-1
            # get max number of turbines that will fit on the circle
            alpha_min = 2.0 * asin.(min_spacing_random / (2.0 * radii[i]))
            max_circle_turbines = minimum([floor.(2.0 * pi / alpha_min), remaining_turbines*0.5])

            min_circle_turbines = 4*i

            # select how many turbines should be placed on the circle
            circle_turbines[i] = rand(min_circle_turbines:max_circle_turbines)

            # update remaining turbines
            remaining_turbines -= circle_turbines[i]
        end

        # set n boundary turbines
        circle_turbines[end] = remaining_turbines

        alphas = 2.0*pi./circle_turbines

        turbinex = zeros(nturbines)
        turbiney = zeros(nturbines)

        index = 1 # start at first turbine
        index += 1 # leave the first turbine at the (0,0)
        for circle in 1:Int(ncircles)
            for turb in 1:Int(circle_turbines[circle])
                angle = alphas[circle]*(turb-1)
                w = radii[circle]*cos(angle)
                h = radii[circle]*sin(angle)
                turbinex[index] = w
                turbiney[index] = h
                index += 1
            end
        end

        # get rotation angle
        step = 0.001
        rotationangle = rand(0:step:(2*pi-step)) # in radians

        # rotate
        turbinex[:], turbiney[:] = rotate_to_wind_direction(turbinex, turbiney, rotationangle)

        # translate
        turbinex .+= center[1]
        turbiney .+= center[2]

    elseif method == "grid"

    elseif method == "vrsunflower"
        # puting random number of turbines on the boundary, distribute the rest using sunflower

        # set minimum number of turbines for the boundary
        min_boundary_turbines = floor(0.3*nturbines)

        # get max number of turbines that will fit on the boundary
        alpha_min = 2.0 * asin.(min_spacing_random / (2.0 * radius))
        max_boundary_turbines = floor.(2.0 * pi / alpha_min)
        if max_boundary_turbines > floor(0.6*nturbines)
            max_boundary_turbines = 0.6*nturbines
        end

        # get number of turbines to put on the boundary
        nbturbines = Int(rand(min_boundary_turbines:max_boundary_turbines))

        # get angle for chosen number of boundary turbines
        alpha_b = 2*pi/nbturbines

        # place boundary turbines
        for turb in 1:nbturbines
            angle = alpha_b*(turb-1)
            w = radius*cos(angle)
            h = radius*sin(angle)
            turbinex[Int(turb)] = w
            turbiney[Int(turb)] = h
        end

        # get number of interior turbines
        niturbines = nturbines - nbturbines

        # place interior turbines
        xi, yi = rotor_sample_points(niturbines, alpha=1).*(radius - min_spacing_random)
        if niturbines > 0
            turbinex[nbturbines+1:end] = xi
            turbiney[nbturbines+1:end] = yi
        end

        # check spacing
        for i = 1:nturbines
            for j = i+1:nturbines
                dist = hypot(turbinex[i]-turbinex[j], turbiney[i]-turbiney[j])

                if dist < min_spacing_random
                    # println(dist)
                    throw(ErrorException("too many turbines to use $method in given space"))
                end
            end
        end

        # get rotation angle
        step = 0.001
        rotationangle = rand(0:step:(2*pi-step)) # in radians

        # rotate
        turbinex[:], turbiney[:] = rotate_to_wind_direction(turbinex, turbiney, rotationangle)

        # translate
        turbinex .+= center[1]
        turbiney .+= center[2]

    elseif method == "sunflower"
        # get locations
        x, y = rotor_sample_points(38, alpha=1.0).*radius

        # check spacing
        for i = 1:nturbines
            for j = i+1:nturbines
                dist = hypot(x[i]-x[j], y[i]-y[j])

                if dist < min_spacing_random
                    # println(dist)
                    throw(ErrorException("too many turbines to use sunflower in given space"))
                end
            end
        end

        # translate
        x .+= center[1]
        y .+= center[2]

        # get rotation angle
        step = 0.001
        rotationangle = rand(0:step:(2*pi-step)) # in radians

        # rotate
        turbinex[:], turbiney[:] = rotate_to_wind_direction(x, y, rotationangle)
    else
        throw(ErrorException("invalid layout generation method: $method"))
    end
    print("\n")

    return turbinex.*rotor_diameter, turbiney.*rotor_diameter
end

function generate_round_layouts(nlayouts, rotor_diameter; farm_center=0., farm_diameter=4000., base_spacing=5., min_spacing=1.,
                           output_directory=nothing, show=false, save_layouts=false, startingindex=1, method="individual")

    if nlayouts > 10 && show == true
        error("do you really want to see $nlayouts plots in serial?")
    end

    boundary_radius = 0.5*(farm_diameter - rotor_diameter)
    # println(boundary_radius)
    area = pi * boundary_radius^2

    turbinex, turbiney = round_farm_concentric_start(copy(rotor_diameter), copy(farm_center), copy(boundary_radius), min_spacing=copy(base_spacing))

    nturbines = size(turbinex)[1]
    # println("nturbines: $nturbines")
    effective_rows = sqrt(nturbines)
    effective_row_length = sqrt(area)
    effective_spacing = effective_row_length / (effective_rows - 1.)
    l = 1

    if save_layouts
        df = DataFrame(turbinex=turbinex./rotor_diameter, turbiney=turbiney./rotor_diameter)
        CSV.write(output_directory*"nTurbs$(nturbines)_spacing$(base_spacing)_layout_$(l+startingindex-1).txt",
                 df, header=["turbineX", "turbineY"])
    end
    if show
        fig, ax = subplots(1)
        plotlayout!(ax, turbinex, turbiney, ones(nturbines).*rotor_diameter)
        plt.show()
    end

    if nlayouts > 1
        for l in 2:nlayouts
            println("Generating Layout $l")

            turbinex, turbiney = round_farm_random_start(copy(rotor_diameter), copy(farm_center),
                                    copy(boundary_radius), min_spacing=copy(base_spacing),
                                    min_spacing_random=min_spacing, method=method, nturbines=nturbines)

            if save_layouts
                df = DataFrame(turbinex=turbinex./rotor_diameter, turbiney=turbiney./rotor_diameter)
                CSV.write(output_directory*"nTurbs$(nturbines)_spacing$(min_spacing)_layout_$(l+startingindex-1).txt",
                           df, header=["turbineX", "turbineY"])
            end
            if show
                fig, ax = plt.subplots(1)
                plotlayout!(ax, turbinex, turbiney, ones(nturbines).*rotor_diameter, xlim=[-boundary_radius-rotor_diameter, boundary_radius+rotor_diameter], ylim=[-boundary_radius-rotor_diameter, boundary_radius+rotor_diameter])
                circle = matplotlib.patches.Circle(farm_center, boundary_radius, fill=false, color="r")
                ax.add_patch(circle)
                plt.show()
            end
        end
    end
end

"""
    pointonline(p, v1, v2; tol=1E-6)

Given a line determined two points (v1 and v2) determine if the point (p) lies on the line
between those points.

Returns true if the point lies on the line (within the given tolerance), false otherwise.

# Arguments
- `p::Vector{Number}(2)`: point of interest
- `v1::Vector{Number}(2)`: first vertex of the line
- `v2::Vector{Number}(2)`: second vertex of the line
- `tol::Number`: how close the cumulative distance from v1 to p to v2 must be to the distance from v1 to v2 to count as being co-linear
"""
function pointonline(p, v1, v2; tol=1E-6)
    # based on https://stackoverflow.com/a/11912171/5128616
    d1 = norm(v1 - p) + norm(p - v2)
    d2 = norm(v1 - v2)
    if isapprox(d1, d2, atol=tol/2.0)
        return true
    else
        return false
    end
end

"""
    pointinpolygon(point, vertices, normals=nothing; s=700, method="raycasting", shift=1E-10, return_distance=true)

Given a polygon determined by a set of vertices, determine the signed distance from the point
to the polygon.

Returns the negative (-) distance if the point is inside or on the polygon, positive (+)
otherwise. If return_distance is set to false, then returns -1 if in polygon or on the
boundary, and 1 otherwise.

# Arguments
- `point::Vector{Number}(2)`: point of interest
- `vertices::Vector{Matrix{Number}(2)`: vertices of polygon
- `normals::Vector{Matrix{Number}(2)`: if not provided, they will be calculated
- `s::Number`: smoothing factor for ksmax function (smoothmax)
- `method::String`: currently only raycasting is available
- `shift::Float`: how far to shift point if it lies on an edge or vertex
- `return_distance::Bool`: if true, return distance. if false, return -1 if in polygon or on the boundary, and 1 otherwise.
"""
function pointinpolygon(point, vertices, normals=nothing; s=700, method="raycasting", shift=1E-10, return_distance=true)
    # println(point)
    if return_distance && typeof(point[1]) <: Int
        throw(ArgumentError("point coordinates may not be given as Ints, must use Floats of some kind. point used $(typeof(point[1]))"))
    end

    if normals === nothing
        normals = boundary_normals_calculator(vertices)
    end

    # number of turbines and boundary vertices
    nvertices = size(vertices)[1]

    # initialize intersection counter
    intersection_counter = 0

    # initialize array to hold distances from each turbine to closest boundary face
    turbine_to_face_distance = zeros(typeof(point[1]), nvertices)

    # get vector from turbine to the first vertex in first face
    turbine_to_first_facepoint = vertices[1, :] - point # dy/dp = -1

    # add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
    vertices = [vertices; vertices[1,1] vertices[1,2]]

    # flag for if in polygon
    inpolygon = false

    # flag for if on vertex
    onvertex = false
    onedge = false

    # initial distance value
    c = 0.0

    # make sure that the point is not exactly on a vertex
    # this check needs to be separate from the edge check because vertex check trumps edge
    # check and the edge check shifts the point slightly.

    for i = 1:nvertices
        # println(point, vertices[i,:])
        if isapprox(point, vertices[i,:], atol=shift/2.0)
            # println("ON VERTEX")
            onvertex = true
            break
            # # if the point is approximately on a vertex or face, move the point slightly
            # # this introduces some slight error, but should be well within the error
            # # for actual turbine placement.

            # # The direction moved is perpendicular to line between the previous and
            # # following vertices to avoid moving along an adjacent face
            # if i == 1
            #     pre_direction_vector = vertices[i+1,:] - vertices[nvertices, :]
            # else
            #     pre_direction_vector = vertices[i+1,:] - vertices[i-1, :]
            # end

            # # get a vector perpendicular to the pre_direction_vector
            # perpendicular_direction = [pre_direction_vector[2], -pre_direction_vector[1]]

            # # normalize perpendicular vector to make it a unit vector
            # # perpendicular_direction ./= norm(perpendicular_direction)
            # perpendicular_direction ./= nansafesqrt(sum(perpendicular_direction.^2))

            # # move the point by shift in the direction of the perpendicular vector
            # point .+= shift*perpendicular_direction

            # break

        elseif pointonline(point, vertices[i,:], vertices[i+1,:], tol=shift/2.0)
            onedge = true
            break
            # # if the point is approximately on a vertex or face, move the point slightly
            # # this introduces some slight error, but should be well within the error
            # # for actual turbine placement.

            # # The direction moved is perpendicular to line between the previous and
            # # following vertices to avoid moving along an adjacent face
            # if i == 1
            #     pre_direction_vector = vertices[i+1,:] - vertices[nvertices, :]
            # else
            #     pre_direction_vector = vertices[i+1,:] - vertices[i-1, :]
            # end

            # # get a vector perpendicular to the pre_direction_vector
            # perpendicular_direction = [pre_direction_vector[2], -pre_direction_vector[1]]

            # # normalize perpendicular vector to make it a unit vector
            # # perpendicular_direction ./= norm(perpendicular_direction)
            # perpendicular_direction ./= nansafesqrt(sum(perpendicular_direction.^2))

            # # move the point by shift in the direction of the perpendicular vector
            # point .+= shift*perpendicular_direction

            # break
            # # if return_distance
            # #     return 0.0
            # # else
            # #     return -1
            # # end

        end
    end

    # check and handle if vertex is on an edge
    # for i = 1:nvertices

    #     onedge = pointonline(point, vertices[i,:], vertices[i+1,:], tol=shift/2.0)
    #     # break
    #     if onedge || onvertex

    #     #     # get a vector for the edge
    #     #     vectoredge = vertices[i+1,:] - vertices[i,:]

    #     #     # get a vector perpendicular to the edge
    #     #     vectorperpendicular = [vectoredge[2], -vectoredge[1]]

    #     #     # get a unit vector perpendicular to the edge
    #     #     vectorperpendicularhat = vectorperpendicular./nansafenorm(vectorperpendicular)

    #     #     # get distance from edge to point
    #     #     c = -abs_smooth(dot(point, vectorperpendicularhat), eps())

    #     #     break
    #     # end
    #     #     if return_distance
    #     #         return 0.0 #-abs_smooth(d, eps())
    #     #     else
    #     #         return -1
    #     #     end
    #     # end
    #         # if the point is approximately on a vertex or face, move the point slightly
    #         # this introduces some slight error, but should be well within the error
    #         # for actual turbine placement.

    #         # The direction moved is perpendicular to line between the previous and
    #         # following vertices to avoid moving along an adjacent face
    #         if i == 1
    #             pre_direction_vector = vertices[i+1,:] - vertices[nvertices, :]
    #         else
    #             pre_direction_vector = vertices[i+1,:] - vertices[i-1, :]
    #         end

    #         # get a vector perpendicular to the pre_direction_vector
    #         perpendicular_direction = [pre_direction_vector[2], -pre_direction_vector[1]]

    #         # normalize perpendicular vector to make it a unit vector
    #         # perpendicular_direction ./= norm(perpendicular_direction)
    #         perpendicular_direction ./= nansafesqrt(sum(perpendicular_direction.^2))

    #         # move the point by shift in the direction of the perpendicular vector
    #         point .+= shift*perpendicular_direction
    #     end
    # end

    # iterate through each boundary
    for j = 1:nvertices

        # check if x-coordinate of turbine is between the x-coordinates of the two boundary vertices
        if real(vertices[j, 1]) <= real(point[1]) < real(vertices[j+1, 1]) || real(vertices[j, 1]) >= real(point[1]) > real(vertices[j+1, 1])

            # check to see if the turbine is below the boundary
            y = (vertices[j+1, 2] - vertices[j, 2]) / (vertices[j+1, 1] - vertices[j, 1]) * (point[1] - vertices[j, 1]) + vertices[j, 2]
            if real(point[2]) < real(y) #(vertices[j+1, 2] - vertices[j, 2]) / (vertices[j+1, 1] - vertices[j, 1]) * (point[1] - vertices[j, 1]) + vertices[j, 2]
                # println((point[1]-vertices[j, 1]))
                # check for case in boundary with sharp angle and point is at same x value as the vertex in the sharp angle
                if isapprox(point[1], vertices[j, 1])
                    # println("here1")
                    if j == 1
                        vprev = vertices[nvertices]
                    else
                        vprev = vertices[j-1]
                    end
                    vnext = vertices[j+1]

                    if !((vprev[1] <= vertices[j, 1] && vnext[1] <= vertices[j, 1]) || (vprev[1] >= vertices[j, 1] && vnext[1] >= vertices[j, 1]))
                        # the vertical ray intersects the boundary
                        intersection_counter += 1
                    end
                else
                    # the vertical ray intersects the boundary
                    intersection_counter += 1
                end
                # intersection_counter += 1

            end
        # elseif vertices[j, 1] == real(point[1])
        #     if j == 1
        #         vprev = vertices[nvertices]
        #     else
        #         vprev = vertices[j-1]
        #     end
        #     vnext = vertices[j+1]

        #     if (vprev[1] < vertices[j, 1] < vnext[1]) || (vprev[1] > vertices[j, 1] > vnext[1])
        #         continue
        #     else
        #         # the vertical ray intersects the boundary
        #         intersection_counter += 1
        #     end
        end

        if return_distance #&& !onedge
            # define the vector from the turbine to the second point of the face
            turbine_to_second_facepoint = vertices[j+1, :] - point # dy/dp = -1

            # find perpendicular distance from turbine to current face (vector projection)

            # get vector defining the boundary
            boundary_vector = vertices[j+1, :] - vertices[j, :]

            # check if perpendicular distance is the shortest
            if real(dot(boundary_vector, -turbine_to_first_facepoint)) > 0 && real(dot(boundary_vector,turbine_to_second_facepoint)) > 0
            # if boundary_vector <= turbine_to_first_facepoint && boundary_vector <= turbine_to_second_facepoint

                # perpendicular distance from turbine to face
                d = dot(turbine_to_first_facepoint, normals[j,:])
                if !onedge && !onvertex
                    turbine_to_face_distance[j] = abs(d)
                else
                    turbine_to_face_distance[j] = abs(d+eps())
                end

            # check if distance to first facepoint is shortest
            elseif real(dot(boundary_vector, -turbine_to_first_facepoint)) < 0

                # distance from turbine to first facepoint
                # turbine_to_face_distance[j] = norm(turbine_to_first_facepoint)
                turbine_to_face_distance[j] = nansafenorm(turbine_to_first_facepoint)

            # distance to second facepoint is shortest
            else

                # distance from turbine to second facepoint
                # turbine_to_face_distance[j] = norm(turbine_to_second_facepoint)
                turbine_to_face_distance[j] = nansafenorm(turbine_to_second_facepoint)

            end

            # reset for next face iteration
            turbine_to_first_facepoint = turbine_to_second_facepoint # dy/dx = 1       # (for efficiency, so we don't have to recalculate for the same vertex twice)
        end
    end

    if return_distance
        # magnitude of the constraint value
        # if !onedge
        #     c = -ff.smooth_max(-turbine_to_face_distance, s=s)
        #     # c = -ksmax(-turbine_to_face_distance, s)
        # end
        c = -ff.smooth_max(-turbine_to_face_distance, s=s)
            # c = -ksmax(-turbine_to_face_distance, s)

        # sign of the constraint value (- is inside, + is outside)
        if mod(intersection_counter, 2) == 1 || onvertex || onedge
            c = -c
        end
    else
        if mod(intersection_counter, 2) == 1 || onedge || onvertex
            c = -1
        else
            c = 1
        end
    end

    return c
end

"""
    nansafesqrt(a)

Calculate the square root of a number, but if the number is less than the given tolerance
then use the line y = a(sqrt(eps())/eps()) so that the derivative is well defined.

# Arguments
- `a::Number`: takes the square root of this value, or approximates it with a line for a < eps()
"""
function nansafesqrt(a::Number)
    tol = eps()
    if real(a) < tol
        return a*sqrt(tol)/tol
    else
        return sqrt(a)
    end

end

"""
    nansafenorm(v)

Calculate the norm of a vector, but if the sum of the squares is less than the given tolerance
then use the line y = a(sqrt(eps())/eps()) so that the derivative is well defined.

# Arguments
- `v::Vector{}`: takes the norm of this vector, but avoids NaN by using a linear
    approximation of sqrt near 0.
"""
function nansafenorm(v::Vector)

    return nansafesqrt(sum(v.^2))

end

"""
    star_boundary(n)

Generate the points for a star with n points

# Arguments
- `n::Int`: The number of points the star should have
- `ri::Float`: inner radius
- `ro::Float`: outer radius
- `rotation::Float`: rotation angle in radians
"""

function star_boundary(n, ri, ro, rotation=0.0)

    # outer angles
    ao = collect(range(2.0*pi/n, 2.0*pi, length=n) .+ rotation)
    # ao .- 2.0*pi/n

    # inner angles
    ai = ao .+ (ao[2] - ao[1])/2.0

    # enforce angles to be between 0 and 2pi
    for i = 1:n
        if ao[i] > 2.0*pi
            ao[i] -= 2.0*pi
        end
        if ai[i] > 2.0*pi
            ai[i] -= 2.0*pi
        end
    end

    # outer x values
    xo = ro.*cos.(ao)

    # outer y values
    yo = ro.*sin.(ao)

    # inner x values
    xi = ri.*cos.(ai)

    # inner y values
    yi = ri.*sin.(ai)

    # generate point set
    vertices = Matrix{Float64}(undef, n*2,2)
    k = 1
    for i = 1:n
        for j = 1:2
            if j == 1
                point = [xo[i] yo[i]]
            else
                point = [xi[i] yi[i]]
            end
            vertices[k,:] = point
            k += 1
        end
    end

    # return
    return round.(vertices, digits=6)
end


text:
@testset "Node composition report" begin

@testset "function nodeCompositionReport" begin
    @test isempty( nodeCompositionReport( mpSim, [-12.0, -8.0], "A junior",
        "Foo", "B senior", "Career", "Branch A", "Bar" ) )
    @test isempty( nodeCompositionReport( mpSim, [12.0, 8.0], "Foo", "Bar" ) )
    @test isempty( nodeCompositionReport( mpSim, -12, "A junior", "Foo",
        "B senior", "Career", "Branch A", "Bar" ) )
    report = nodeCompositionReport( mpSim, 12, "A junior", "Empty", "Foo",
        "B senior", "Career", "Branch A", "Bar" )
    @test haskey( report, "Base nodes" ) &&
        ( ["A junior", "B senior"] ⊆ names( report["Base nodes"] ) )
    @test all( [haskey( report, "Career" ), haskey( report, "Branch A" ),
        !haskey( report, "Foo" ), !haskey( report, "Bar" ),
        haskey( report, "Empty" )] )
    @test ( size( report["Empty"], 2 ) == 2 ) &&
        all( report["Empty"][:, "Empty"] .== 0 )
end  # @testset "function nodeCompositionReport"

end  # @testset "Node composition report"

ext: jl


text:
#=
fld-write.jl
write AVD .fld to file
2019-05-12, Jeff Fessler, University of Michigan
=#

export fld_write
export ir_test_dir, ir_test_dir!

"""
`dir = ir_test_dir!(path)`

set default test directory
"""
function ir_test_dir!(path::String)
	ENV["MIRT_TEST_DIR"] = path
end

"""
`dir = ir_test_dir`

return default test directory `ENV["MIRT_TEST_DIR"]`; default is `/tmp`
"""
function ir_test_dir()
#	get(ENV, "MIRT_TEST_DIR", "/tmp")
	try
		return ENV["MIRT_TEST_DIR"]
	catch
		return "/tmp"
	end
end

"""
`fld_write(file, data, varargin)`

write data into AVS format `.fld` file

See ASPIRE user's guide (or AVS online documentation) for file format.

in
- `file` name of file typically ending in `.fld`
- `data` real data array

option
- `check::Bool`			report error if file exists; default `true`
- `dir::String`			directory name to prepend file name; default `""`
- `endian::`Symbol`		`:le` little endian (default), `:be` big endian
- `head::Array{String}`	comment information for file header
- `raw::Bool`			put raw data in `name.raw`, header in `name.fld`
						where `file` = `name.fld`; default `false`
"""

# * `how::Symbol`			file data type; default `:data` infer from data

function fld_write(file::String, data::AbstractArray{<:Real};
	check::Bool = true,
	dir::String = "",
	endian::Symbol = :le,
	head::Array{String} = empty([""]),
#	how::Symbol = :data,
	raw::Bool = false,
	)

	data = fld_write_data_fix(data) # make data suitable for writing in .fld

	endian != :le && endian != :be && throw("endian '$endian' unknown")

	typedict = Dict([
		(Float32, endian == :le ? "float_le" : "xdr_float"),
		(Float64, endian == :le ? "double_le" : "xdr_double"),
		(UInt8, "byte"),
		(Int16, endian == :le ? "short_le" : "short_be"),
		(Int32, endian == :le ? "int_le" : "xdr_int"),
		])

	dtype = eltype(data)
	datatype = typedict[dtype] # throws error if not there

	file = joinpath(dir, file)
#	@show file

	check && isfile(file) && throw("file '$file' exists")

	if raw # if writing data to separate ".raw" file, ensure it does not exist
		fileraw = file[1:(end-3)] * "raw"
		check && isfile(fileraw) && throw("file $fileraw exists")
	end

	# open output avs file for writing
	fid = open(file, "w")
	fraw = fid # default
	if raw
        fraw = open(fileraw, "w")
	end

	# write header
	ndim = ndims(data)

	println(fid, "# AVS field file ($(basename(@__FILE__)))")
	for ii=1:length(head)
		println(fid, "# $(head[ii])")
	end

	println(fid, "ndim=$ndim")
	for ii=1:ndim
		println(fid, "dim$ii=$(size(data,ii))")
	end
	println(fid, "nspace=$ndim")
	println(fid, "veclen=1")
	println(fid, "data=$datatype")
	println(fid, "field=uniform")

	if raw
		println(fid, "variable 1 file=$fileraw filetype=binary")
	else
		write(fid, "\f\f") # two form feeds: char(12)
	end

	# finally, write the binary data
	host_is_le = ENDIAN_BOM == 0x04030201
	if host_is_le == (endian == :le) # host/file same endian
		write(fraw, data) # finally, write the binary data
	elseif host_is_le && (endian == :be)
		write(fraw, hton.(data))
	elseif !host_is_le && (endian == :le)
		write(fraw, htol.(data))
	end

	close(fid)

	if raw
		close(fraw)
	end

end

"""
`data = fld_write_data_fix(data)`
convert data to format suitable for writing to .fld file
"""
function fld_write_data_fix(data::AbstractArray{<:Real})

	dtype = eltype(data)

	if dtype == BigFloat
		@warn "BigFloat downgraded to Float64"
		return Float64.(data)
	end

	if dtype == Float16
		@warn "Float16 upgraded to Float32"
		return Float32.(data)
	end

	if dtype == BigInt
		@warn "BigInt downgraded to Int32"
		return Int32.(data)
	end

	if dtype == Int128
		@warn "Int128 downgraded to Int64"
		return Int64.(data)
	end

	if dtype == Int64
		@warn "Int64 downgraded to Int32"
		return Int32.(data)
	end

	if eltype(data) == Bool
		@warn "Bool upgraded to UInt8"
		return UInt8.(data)
	end

	!in(dtype, (Float32, Float64, UInt8, Int16, Int32)) &&
		throw("unsupported type $dtype")

	return data
end

"""
`fld_write_test1(file, data, ...)`
"""
function fld_write_test1(file, data; raw::Bool=false, chat::Bool=false, kwarg...)
	fld_write(file, data; raw=raw, kwarg...)
	tmp = fld_read(file, chat=chat)
	tmp != data && throw("test failed for file = $file")

	rm(file)
	if raw
		tmp = file[1:(end-4)] * ".raw"
		rm(tmp)
	end
	true
end

"""
`fld_write(:test ; chat::Bool)`

tests writing and reading AVS `.fld` files of all key types
"""
function fld_write(test::Symbol; chat::Bool=false)
	test != :test && throw("bug")

	file = joinpath(ir_test_dir(), "fld-write-test.fld")

	for dtype in (UInt8, Int16, Int32, Float32, Float64)
		for endian in (:le, :be)
			for raw in (false, true)
				chat && @info "dtype=$dtype endian=$endian raw=$raw"
				data = convert.(dtype, [5:8; 1:4])
				fld_write_test1(file, data, endian=endian, raw=raw,
					check=true, chat=chat)
			end
		end
	end

	true
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectivityDestination <: SwaggerModel
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    address::Any # spec type: Union{ Nothing, String } # spec name: address
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port

    function ConnectivityDestination(;resourceId=nothing, address=nothing, port=nothing)
        o = new()
        validate_property(ConnectivityDestination, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(ConnectivityDestination, Symbol("address"), address)
        setfield!(o, Symbol("address"), address)
        validate_property(ConnectivityDestination, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        o
    end
end # type ConnectivityDestination

const _property_map_ConnectivityDestination = Dict{Symbol,Symbol}(Symbol("resourceId")=>Symbol("resourceId"), Symbol("address")=>Symbol("address"), Symbol("port")=>Symbol("port"))
const _property_types_ConnectivityDestination = Dict{Symbol,String}(Symbol("resourceId")=>"String", Symbol("address")=>"String", Symbol("port")=>"Int32")
Base.propertynames(::Type{ ConnectivityDestination }) = collect(keys(_property_map_ConnectivityDestination))
Swagger.property_type(::Type{ ConnectivityDestination }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectivityDestination[name]))}
Swagger.field_name(::Type{ ConnectivityDestination }, property_name::Symbol) =  _property_map_ConnectivityDestination[property_name]

function check_required(o::ConnectivityDestination)
    true
end

function validate_property(::Type{ ConnectivityDestination }, name::Symbol, val)
end


text:
function f2(x)
  x[1]^2 + (2.0 - x[2])^2
end
function g2(x, storage)
  storage[1] = 2.0 * x[1]
  storage[2] = -2.0 * (2.0 - x[2])
end
d2 = DifferentiableFunction(f2, g2)
initial_x = [100.0, 100.0]

results = Optim.bfgs(d2, initial_x)
@assert length(results.trace.states) == 0
@assert results.converged
@assert norm(results.minimum - [0.0, 2.0]) < 0.01


text:
# Deprecate the old name AbstractMCMCParallel in favor of AbstractMCMCEnsemble
Base.@deprecate_binding AbstractMCMCParallel AbstractMCMCEnsemble false


text:
module GeoInterface

using Extents
using Base.Iterators: flatten

export testgeometry, isgeometry, trait, geomtrait, ncoord, getcoord, ngeom, getgeom

# traits
export AbstractGeometryTrait,
    AbstractGeometryCollectionTrait,
    GeometryCollectionTrait,
    AbstractPointTrait,
    PointTrait,
    AbstractCurveTrait,
    AbstractLineStringTrait,
    LineStringTrait,
    LineTrait,
    LinearRingTrait,
    CircularStringTrait,
    CompoundCurveTrait,
    AbstractSurfaceTrait,
    AbstractCurvePolygonTrait,
    CurvePolygonTrait,
    AbstractPolygonTrait,
    PolygonTrait,
    TriangleTrait,
    RectangleTrait,
    QuadTrait,
    PentagonTrait,
    HexagonTrait,
    AbstractPolyhedralSurfaceTrait,
    PolyhedralSurfaceTrait,
    TINTrait,
    AbstractMultiPointTrait,
    MultiPointTrait,
    AbstractMultiCurveTrait,
    MultiCurveTrait,
    AbstractMultiLineStringTrait,
    MultiLineStringTrait,
    AbstractMultiSurfaceTrait,
    MultiSurfaceTrait,
    AbstractMultiPolygonTrait,
    MultiPolygonTrait,
    AbstractFeatureTrait,
    FeatureTrait,
    AbstractCollectionFeatureTrait,
    FeatureCollectionTrait

include("types.jl")
include("interface.jl")
include("fallbacks.jl")
include("utils.jl")
include("base.jl")

end # module


text:
using Pkg
Pkg.add("MultivariateStats")
using Optim
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using DataFramesMeta
using CSV
using MultivariateStats
using HTTP
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://github.com/OU-PhD-Econometrics/fall-2021/blob/master/ProblemSets/PS8-factor/nlsy.csv"

using DataFrames;
df = CSV.read(HTTP.get(url).body, DataFrame)
ols = glm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr), df)
print(ols)
#I tried to make this code run but the variables are string and I read we need to use Dict to make the actual variable names appear so it wouldn't give an error but I don't know how to apply Dict to a code#

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
#I think we would make a matrix of the asvab variables and then run a correlation of the new matrix which would be a submatrix. I had a bit of difficulty writing the code here#

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::
olsasvab= glm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabAR + asvabCS + asvabMK + asvabNO + asvabPC + asvabWK), df)
print(olsasvab)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::
M=fit(PCA, asvabMat; maxoutdim=1)
asvabPCA=MultivariateStats.transform(M,asvabMat)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::
M=fit(FactorAnalysis, asvabMat; maxoutdim=1)
asvabFactorAnalysis=MultivariateStats.transform(M,asvabMat)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 6
#:::::::::::::::::::::::::::::::::::::::::::::::::::


text:
module GMRF

using LinearAlgebra, SparseArrays, StatsBase, Distributions

import Distributions.rand, Distributions.logpdf

include("graphstructure.jl")
include("gridstructure.jl")
include("igmrf.jl")

export iGMRF, rand, logpdf, fullconditionals, fullcondlogpdf, getconditional

end # module


text:
module MicrobiomeTests

using Microbiome
using ReTest
using Microbiome.SparseArrays
using Microbiome.Tables
using Microbiome.Dictionaries
import Microbiome.MultivariateStats: MDS

@testset "Samples and Features" begin
    @testset "MicriobiomeSamples and metadata" begin
        ms = MicrobiomeSample("sample")
        @test name(ms) == "sample"
        @test isempty(metadata(ms))
        @test metadata(ms) isa Dictionary
        @test_throws Dictionaries.IndexError ms.thing = "metadata"
        @test_throws Dictionaries.IndexError ms[:thing] = "metadata"

        @test insert!(ms, :thing, "metadata") isa MicrobiomeSample
        @test ms.thing == ms[:thing] == "metadata"
        @test let
            ms.thing = "other metadata"
            @test ms.thing == ms[:thing] == "other metadata"
            ms[:thing] = "metadata"
            ms.thing == ms[:thing] == "metadata"
        end
        @test_throws MethodError ms["thing"] = "metadata"
        @test_throws Dictionaries.IndexError insert!(ms, :thing, "still other metadata")
        @test_throws Dictionaries.IndexError delete!(ms, :thing2)
        @test let
            delete!(ms, :thing)
            unset!(ms, :thing2)
            true
        end
        @test let
            set!(ms, :thing2, "metadata2")
            ms.thing2 == "metadata2"
        end

        # kwargs constructor
        ms2 = MicrobiomeSample("sample2"; age=10, birthtype="vaginal", allergies=true)
        # Dict constructor
        ms3 = MicrobiomeSample("sample3", Dict(:age=>10, :birthtype=>"vaginal", :allergies=>true))
        # NamedTuple constructor
        ms4 = MicrobiomeSample("sample4", (;age=10, birthtype="vaginal", allergies=true))

        for ms in [ms2, ms3, ms4]
            @test ms.age == 10
            @test ms.birthtype == "vaginal"
            @test Bool(ms.allergies)

            @test_throws ArgumentError insert!(ms, (; birthtype="cesarean"))
            insert!(ms, (; foo=10))
            @test ms.foo == 10
            set!(ms, (; birthtype="cesarean"))
            @test ms.birthtype == "cesarean"
        end

    end

    @testset "Taxa" begin
        txm = Taxon("taxon", missing)
        @test txm === Taxon("taxon")
        @test ismissing(taxrank(txm))
        @test !hasrank(txm)

        for (r, i) in pairs(Microbiome._ranks)
            tx = Taxon("taxon", r)
            @test taxrank(tx) == r
            @test tx === Taxon("taxon", i)
            @test tx !== txm
        end

        @test_throws ErrorException Taxon("taxon", :invalid)
        @test_throws ErrorException Taxon("taxon", 10)
        @test let tx = Taxon("taxon", :kingdom)
            hasrank(tx)
            String(tx) == string(tx) == "k__taxon"
            taxon("g__Somegenus") == Taxon("Somegenus", :genus)
        end

    end

    @testset "Gene Functions" begin
        gfm = GeneFunction("gene", missing)
        @test name(gfm) == "gene"
        @test gfm === GeneFunction("gene")
        @test ismissing(taxon(gfm))
        @test !hastaxon(gfm)
        @test !hasrank(gfm)

        gf1 = GeneFunction("gene", Taxon("sp1", :species))
        gf2 = GeneFunction("gene", Taxon("sp1"))
        @test name(gf1) == "gene"
        @test String(gf1) == string(gf1) == "gene|s__sp1"
        @test String(gf2) == string(gf2) == "gene|u__sp1"
        @test genefunction("gene|s__sp1") == gf1
        @test genefunction("gene|u__sp1") == gf2

        @test gf1 != gfm
        @test gf1 === GeneFunction("gene", Taxon("sp1", :species))
        @test gf2 === GeneFunction("gene", Taxon("sp1"))
        @test gf1 !== gf2
        @test hastaxon(gf1)
        @test !ismissing(taxon(gf1))
        @test taxon(gf1) == Taxon("sp1", :species)
        @test taxon(gf1) != taxon(gf2)
        @test hasrank(gf1)
        @test taxrank(gf1) == :species
    end
end

@testset "Profiles" begin
    _ranks = Tuple(keys(Microbiome._ranks))[1:9]
    mss = [MicrobiomeSample("sample$i") for i in 1:5]
    txs = [Taxon("taxon$i", _ranks[i]) for i in 1:9]
    push!(txs, Taxon("taxon10", missing))

    mat = spzeros(10,5)
    dmat = zeros(10,5)
    for i in 1:5;
        mat[i,i] = 1.
        dmat[i,i] = 1.
    end
    for i in 1:5
        mat[i+5,i] = 0.6
        dmat[i+5,i] = 0.6
    end

    comm = CommunityProfile(mat, txs, mss)

    @testset "Profile operations" begin
        @test CommunityProfile(mat, txs, mss) isa CommunityProfile
        @test comm == CommunityProfile(dmat, txs, mss)

        @test nsamples(comm) == 5
        @test nfeatures(comm) == 10
        @test size(comm) == (10, 5)
        @test profiletype(comm) == Taxon
        @test ranks(comm)[1:9] == [_ranks...]
        @test sampletotals(comm) == [1.6 1.6 1.6 1.6 1.6]
        @test featuretotals(comm) == reshape([1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 0.6, 0.6, 0.6, 0.6], 10, 1)
        @test features(comm) == txs
        @test samples(comm) == mss

        @test size(rankfilter(comm, :species), 1) == 1
        @test size(rankfilter(comm, :genus; keepempty=true), 1) == 2

        @test size(rankfilter(comm, 5), 1) == 1
        @test size(rankfilter(comm, 6; keepempty=true), 1) == 2
        @test_throws ErrorException rankfilter(comm, :foo)
        @test_throws ErrorException rankfilter(comm, 10)
        @test_throws ErrorException rankfilter(rankfilter(comm, :species), :genus) # will be empty

        @test featurenames(filter(f-> hasrank(f) && taxrank(f) == :species, comm)) == featurenames(rankfilter(comm, :species))

        @test present(0.1)
        @test ismissing(present(missing))
        @test present(comm) == sparse(reshape([i==j || i+5==j for i in 1:5 for j in 1:10], 10,5))
        @test !present(0.1, 0.2)
        @test_throws DomainError present(-0.1)
        @test_throws DomainError present(0., -0.1)

        @test prevalence([0.0, 0.1, 0.2, 0.3]) ≈ 0.75
        @test prevalence([0.0, 0.1, 0.2, 0.3], 0.15) ≈ 0.5
        @test prevalence((0.0, 0.1, 0.2, 0.3)) ≈ 0.75
        @test prevalence((0.0, 0.1, 0.2, 0.3), 0.15) ≈ 0.5

        @test all(≈(0.2), prevalence(comm))
        @test all(prevalence(comm, 0.7) .≈ [0.2, 0.2, 0.2, 0.2, 0.2, 0, 0, 0, 0, 0])

        @test let c2 = deepcopy(comm)
            relativeabundance!(c2)
            @test all(≈(0.625), featuretotals(c2)[1:5])
            all(≈(0.375), featuretotals(c2)[6:end])
        end
        @test let c2 = deepcopy(comm)
            relativeabundance!(c2, kind=:percent)
            @test all(≈(62.5), featuretotals(c2)[1:5])
            all(≈(37.5), featuretotals(c2)[6:end])
        end

        @test_throws ArgumentError relativeabundance!(comm, kind=:invalid)

        @test_throws ErrorException commjoin(comm, comm)
        let c3 = commjoin(comm[:,1:2], comm[:, 3:4], comm[:, 5])
            @test all(abundances(c3) .== abundances(comm))
            @test all(samples(c3) .== samples(comm))
            @test all(features(c3) .== features(comm))
        end

        filtertest = CommunityProfile(sparse(Float64[3 2 1 # 0.66, assuming minabundance 2
                                              2 2 2 # 1.0
                                              0 0 1 # 0.0
                                              2 0 0 # 0.33
                                       ]),
                               [Taxon(string(i)) for i in 1:4],
                               [MicrobiomeSample(string(i)) for i in 1:3]);
        @test size(prevalence_filter(filtertest)) == (4,3)
        @test size(prevalence_filter(filtertest, minabundance=2)) == (3,3)
        @test size(prevalence_filter(filtertest, minabundance=2, minprevalence=0.4)) == (2,3)
        @test all(<=(1.0), abundances(prevalence_filter(filtertest, renorm=true)))
        @test all(x-> isapprox(x, 1.0, atol=1e-8), sum(abundances(prevalence_filter(filtertest, renorm=true)), dims=1))
    end

    @testset "Stratified gene functions operations" begin
        strat = CommunityProfile(sparse(Float64[3 2 1
                                                2 2 2
                                                0 0 1
                                                0 1 0
                                        ]),
                                        [GeneFunction("gene1"),
                                         GeneFunction("gene1", "species1"),
                                         GeneFunction("gene1", "species2"),
                                         GeneFunction("gene2")
                                         ],
                                        [MicrobiomeSample(string(i)) for i in 1:3]);

        @test filter(hastaxon, strat)         |> nfeatures == 2
        @test filter(!hastaxon, strat)        |> nfeatures == 2
        @test strat["gene1", :]               |> nfeatures == 1
        @test strat[["gene1", "gene2"], :]    |> nfeatures == 2
        @test strat[r"gene1", :]               |> nfeatures == 3
        @test strat[r"gene[12]", :]    |> nfeatures == 4
        @test strat[GeneFunction("gene1"), :] |> nfeatures == 1
    end

    @testset "Profile Metadata" begin
        s1 = MicrobiomeSample("sample1", Dictionary(Dict(:age=> 37, :name=>"kevin", :something=>1.0)))
        s2 = MicrobiomeSample("sample2", Dictionary(Dict(:age=> 37, :name=>"kevin", :something_else=>2.0, :still_other=>"boo")))

        @testset "Single sample" begin
            c4 = CommunityProfile(sparse([1 1; 2 2; 3 3]), [Taxon(string(i)) for i in 1:3], [s1, s2])
            md1, md2 = metadata(c4)

            @test_throws Dictionaries.IndexError insert!(c4, "sample1", :something, 3.0)
            @test_throws Dictionaries.IndexError delete!(c4, "sample1", :something_else)

            @test delete!(c4, "sample1", :something) isa MicrobiomeSample
            @test !haskey(c4, "sample1", :something)
            @test get(c4, "sample1", :something_else, 42) == 42
            @test ismissing(get(c4, "sample1", :something_else))
            @test insert!(c4, "sample1", :something, 3.0) isa MicrobiomeSample
            @test get(c4, "sample1", :something, 42) == 3.0
            set!(c4, "sample1", :something, 1.0)
            @test first(metadata(c4))[:something] == 1.0
            @test haskey(c4, "sample1", :something)
            @test unset!(c4, "sample1", :something) isa MicrobiomeSample
            @test !haskey(c4, "sample1", :something)
            set!(c4, "sample1", :something, 1.0)

            @test collect(keys(c4, "sample1")) == [:age, :name, :something]
        end

        @testset "Whole community" begin
            c5 = CommunityProfile(sparse([1 1; 2 2; 3 3]), [Taxon(string(i)) for i in 1:3], [s1, s2])
            md1, md2 = metadata(c5)

            @test all(row-> row[:age] == 37, [md1, md2])
            @test all(row-> row[:name] == "kevin", [md1, md2])
            @test md1[:something] == 1.0
            @test ismissing(md2[:something])
            @test md2[:something_else] == 2.0
            @test ismissing(md1[:something_else])

            @testset "insert!" begin
                insert!(c5, "sample1", :foo, "bar")
                @test samples(c5, "sample1").foo == "bar"
                @test_throws IndexError insert!(c5, "sample1", :foo, "baz")

                insert!(c5, "sample2", (; foo="baz", greeting="hello"))
                @test samples(c5, "sample2").foo == "baz"
                @test_throws IndexError insert!(c5, "sample2", Dict(:baz=> "test", :foo=>"bar"))
                @test !haskey(c5, "sample2", :baz)

                insert!(c5, "sample2", (; graw="gnaw", biff="boof"))
                @test_throws IndexError insert!(c5, [(;sample="sample1", still_other="yes"), (;sample="sample2", still_other="no")])
                @test !haskey(c5, "sample1", :still_other)
                @test all(get(c5, :something) .=== [1.0, missing])
                @test all(get(c5, :something, 42.0) .== [1.0, 42.0])
            end

            @testset "set!" begin
                set!(c5, "sample1", :foo, "barre")
                @test samples(c5, "sample1").foo == "barre"
                set!(c5, "sample2", (; foo="bazze", greeting="hello world!"))
                @test samples(c5, "sample2").foo == "bazze"

                set!(c5, [(;sample="sample1", still_other="yes"), (;sample="sample2", still_other="no")])
                @test haskey(c5, "sample1", :still_other)
                @test samples(c5, "sample1").still_other == "yes"
                @test samples(c5, "sample2").still_other == "no"
            end
        end
    end

    @testset "Indexing and Tables integration" begin
        @test Tables.istable(comm)
        @test Tables.columnaccess(comm)
        @test Tables.rowaccess(comm)
        @test Tables.schema(comm) isa Tables.Schema

        for i in 1:5
            @test abundances(comm[:, "sample$i"]) == mat[:, [i]]
            @test abundances(comm["$(keys(Microbiome._shortranks)[i])__taxon$i", :]) == mat[[i], :]
        end

        @test abundances(comm[r"taxon1", :]) == abundances(comm[["d__taxon1", "u__taxon10"], :]) == abundances(comm[[1,10], :])
        @test abundances(comm[:, r"sample[13]"]) == abundances(comm[:,["sample1", "sample3"]]) == abundances(comm[:, [1,3]])
        @test abundances(comm[r"taxon1", r"sample[13]"]) ==
              abundances(comm[["d__taxon1", "u__taxon10"],["sample1", "sample3"]]) ==
              abundances(comm[["d__taxon1", "u__taxon10"],[r"sample1", r"sample3"]]) ==
              abundances(comm[[1,10], [1,3]])

        for (i, col) in enumerate(Tables.columns(comm))
            if i == 1
                @test col == txs
            else
                @test col ==  mat[:, [i-1]]
            end
        end

        for (i, row) in enumerate(Tables.rows(comm))
            @test row == (; :features => txs[i], (Symbol("sample$(j)") => mat[i, j] for j in 1:5)...)
        end

        tbl = Tables.columntable(comm)
        @test tbl.features == features(comm)
        @test keys(tbl) == (:features, Symbol.(samplenames(comm))...)
        @test tbl.sample1 == abundances(comm[:, "sample1"])
    end

    @testset "Diversity" begin
        R = 10
        s1 = collect(0:R-1) # high diversity
        s2 = [i % 2 == 0 ? s1[i] : 0 for i in eachindex(s1)] # low diversity
        s3 = ones(R) # uniform
        s4 = [10, zeros(R-1)...] # no diversity

        @test shannon(s1) > shannon(s2)
        @test shannon(s3) ≈ log(R)
        @test shannon(s4) ≈ 0.0
        for s in (s1, s2, s3, s4)
            shannon(s ./ sum(s)) ≈ shannon(s)
        end

        @test ginisimpson(s1) > ginisimpson(s2)
        @test ginisimpson(s3) ≈ 1.0 - 1/R
        @test ginisimpson(s4) ≈ 0.0
        for s in (s1, s2, s3, s4)
            ginisimpson(s ./ sum(s)) ≈ ginisimpson(s)
        end

        @test size(shannon(comm)) == (1, 5)
        @test size(ginisimpson(comm)) == (1, 5)

        @test let c2 = deepcopy(comm)
            shannon!(c2)
            @test all(s-> haskey(s, :shannon), samples(c2))
            @test_throws Dictionaries.IndexError shannon!(c2)
            shannon!(c2, overwrite=true)
            true
        end
        @test let c2 = deepcopy(comm)
            ginisimpson!(c2)
            @test all(s-> haskey(s, :ginisimpson), samples(c2))
            @test_throws Dictionaries.IndexError ginisimpson!(c2)
            ginisimpson!(c2, overwrite=true)
            true
        end

        @test all(braycurtis(comm) .== [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0])
        @test all(jaccard(comm) .== [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0])
        @test all(hellinger(comm) .== [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0])
        @test pcoa(comm) isa MDS
    end
end

end # module


text:
module SigmoidNumbers

  include("SigmoidTypedef.jl")
  include("macros.jl")
  include("inttools.jl")
  include("sigtools.jl")
  include("h-layer.jl")

  include("Math/iterators.jl")
  include("Math/cmp.jl")
  include("Math/constants.jl")
  include("Math/addsub.jl")
  include("Math/muldiv.jl")

  include("convert.jl")

  include("Posits/posits.jl")
  include("Posits/convert.jl")
  include("Posits/mlfunctions.jl")

  include("Valids/valids.jl")

  include("Tools/fft.jl")
  include("Tools/ubox.jl")

  include("quire.jl")
  include("fdp-matrixsolve.jl")
  include("quire_solve.jl")

  #patching julia's native lu factorization.
  include("LUpatches.jl")

  #including blas functions
  include("BLAS/blas.jl")
  include("BLAS/Level1.jl")
  include("BLAS/Level2-er.jl")
  include("BLAS/Level2-mv.jl")
  include("BLAS/Level2-sv.jl")

end # module

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
export InDims, OutDims,AsArray,AsDataFrame,AsAxisArray
const AxisDescriptorAll = Union{AxisDescriptor,String,Type{T},CubeAxis,Function} where T<:CubeAxis
import ..Cubes.Axes: get_descriptor, ByFunction
import ..Cubes: getsavefolder
import Zarr: Compressor, NoCompressor
import ...ESDL: workdir, ESDLDefaults

abstract type ArTypeRepr end
struct AsArray <: ArTypeRepr end
#struct AsAxisArray <: ArTypeRepr end
struct AsDataFrame <: ArTypeRepr
  dimcol::Bool
end
AsDataFrame()=AsDataFrame(false)
wrapWorkArray(::AsArray,a,axes) = a
# function wrapWorkArray(::AsAxisArray,a,cablabaxes)
#   newaxes = map(cablabaxes) do ax
#     AxisArrays.Axis{Symbol(axname(ax))}(ax.values)
#   end
#   AxisArrays.AxisArray(a,newaxes...)
# end
import DataFrames
function wrapWorkArray(t::AsDataFrame,a,cablabaxes)
  colnames = map(Symbol,cablabaxes[2].values)
  df = DataFrames.DataFrame(a,colnames)
  if t.dimcol
    df[Symbol(axname(cablabaxes[1]))]=collect(cablabaxes[1].values)
  end
  df
end

abstract type ProcFilter end
struct AllMissing <: ProcFilter end
struct NValid     <: ProcFilter
  n::Int
end
struct AnyMissing <: ProcFilter end
struct AnyOcean   <: ProcFilter end
struct NoFilter   <: ProcFilter end
struct StdZero    <: ProcFilter end
struct UserFilter{F} <: ProcFilter
  f::F
end

checkskip(::NoFilter,x)         = false
checkskip(::AllMissing,x::AbstractArray)   = all(ismissing,x)
checkskip(::AllMissing,df::DataFrame)  = any(map(i->all(ismissing,getindex(df,i)),names(df)))
checkskip(::AnyMissing,x::AbstractArray)   = any(ismissing,x)
checkskip(::AnyMissing,df::DataFrame)  = any(map(i->any(ismissing,getindex(df,i)),names(df)))
checkskip(nv::NValid,x::AbstractArray)     = count(!ismissing,x) <= nv.n
checkskip(uf::UserFilter,x) = uf.f(x)
checkskip(::StdZero,x)      = all(i->i==x[1],x)
docheck(pf::ProcFilter,x)::Bool = checkskip(pf,x)
docheck(pf::Tuple,x)        = reduce(|,map(i->docheck(i,x),pf))

getprocfilter(f::Function) = (UserFilter(f),)
getprocfilter(pf::ProcFilter) = (pf,)
getprocfilter(pf::NTuple{N,<:ProcFilter}) where N = pf

"""
    InDims(axisdesc;...)

Creates a description of an Input Data Cube for cube operations. Takes a single
  or a Vector/Tuple of axes as first argument. Axes can be specified by their
  name (String), through an Axis type, or by passing a concrete axis.

### Keyword arguments

* `artype` how shall the array be represented in the inner function. Defaults to `AsArray`, alternatives are `AsDataFrame` or `AsAxisArray`
* `filter` define some filter to skip the computation, e.g. when all values are missing. Defaults to
    `AllMissing()`, possible values are `AnyMissing()`, `AnyOcean()`, `StdZero()`, `NValid(n)`
    (for at least n non-missing elements). It is also possible to provide a custom one-argument function
    that takes the array and returns `true` if the compuation shall be skipped and `false` otherwise.
"""
mutable struct InDims
  axisdesc::Tuple
  artype::ArTypeRepr
  procfilter::Tuple
end
function InDims(axisdesc::AxisDescriptorAll...; artype::ArTypeRepr=AsArray(), filter = AllMissing())
  descs = map(get_descriptor,axisdesc)
  any(i->isa(i,ByFunction),descs) && error("Input cubes can not be specified through a function")
  isa(artype,AsDataFrame) && length(descs)!=2 && error("DataFrame representation only possible if for 2D inner arrays")
  InDims(descs,artype,getprocfilter(filter))
end

"""
    OutDims(axisdesc;...)

Creates a description of an Output Data Cube for cube operations. Takes a single
  or a Vector/Tuple of axes as first argument. Axes can be specified by their
  name (String), through an Axis type, or by passing a concrete axis.

- `axisdesc`: List of input axis names
- `genOut`: function to initialize the values of the output cube given its element type. Defaults to `zero`
- `finalizeOut`: function to finalize the values of an output cube, defaults to identity.
- `retCubeType`: sepcifies the type of the return cube, can be `CubeMem` to force in-memory, `TempCube` to force disk storage, or `"auto"` to let the system decide.
- `chunksize`: Chunk size for the inner dimensions, a tuple of the same length as `axisdesc`, or `:input` to copy chunksizes from input cube axes or `:max` to not chunk the inner dimensions
- `compressor`: A Zarr compressor for the specified output cube
- `retCubeType`: sepcifies the type of the return cube, can be `CubeMem` to force in-memory, `ZArrayCube` to force disk storage, or `"auto"` to let the system decide.
- `outtype`: force the output type to a specific type, defaults to `Any` which means that the element type of the first input cube is used
"""
struct OutDims
  axisdesc
  bcaxisdesc
  genOut::Any
  finalizeOut::Any
  retCubeType::Any
  update::Bool
  artype::ArTypeRepr
  chunksize::Any
  compressor::Any
  path::String
  persist::Bool
  outtype::Union{Int,DataType}
end
function OutDims(axisdesc...;
           bcaxisdesc=(),
           genOut=zero,
           finalizeOut=identity,
           retcubetype=:auto,
           update=false,
           artype::ArTypeRepr=AsArray(),
           chunksize=ESDLDefaults.chunksize[],
           compressor=ESDLDefaults.compressor[],
           path="",
           outtype=1)
  descs = map(get_descriptor,axisdesc)
  bcdescs = (map(get_descriptor,bcaxisdesc)...,)
  isa(artype,AsDataFrame) && length(descs)!=2 && error("DataFrame representation only possible if for 2D inner arrays")
  if !in(chunksize,(:input, :max)) && length(chunksize)!=length(axisdesc)
    error("Length of chunk sizes must equal number of inner Axes")
  end
  if path == ""
    persist = false
  else
    persist = true
  end
  OutDims(descs,bcdescs,genOut,finalizeOut,retcubetype,update,artype,chunksize,compressor,path,persist,outtype)
end

registerDATFunction(a...;kwargs...)=@warn("Registration does not exist anymore, ignoring....")


text:
#all OpenGL Types
typealias uint_least64_t                                              Uint64
typealias uint_least16_t                                              Uint16
typealias uint_fast32_t                                               Uint64
typealias int_fast8_t                                                 Uint8
typealias int_fast16_t                                                Int64
typealias int_least64_t                                               Int64
typealias GLubyte                                                     Cuchar
typealias int16_t                                                     Int16
typealias GLfloat                                                     Cfloat
typealias int8_t                                                      Uint8
typealias GLenum                                                      Cuint
typealias GLboolean                                                   Cuchar
typealias uint16_t                                                    Uint16
typealias GLclampd                                                    Cdouble
typealias GLvoid                                                      Void
typealias uint_fast8_t                                                Uint8
typealias int64_t                                                     Int64
typealias int_fast64_t                                                Int64
typealias uint32_t                                                    Uint32
typealias GLint                                                       Cint
typealias uint_least32_t                                              Uint32
typealias GLsizeiptr                                                  Cint
typealias ptrdiff_t                                                   Int64
typealias int_least32_t                                               Int32
typealias GLcharARB                                                   Uint8
typealias GLuint                                                      Cuint
typealias GLbyte                                                      Cuchar
typealias size_t                                                      Int64
typealias GLchar                                                      Cchar
typealias uintptr_t                                                   Uint64
typealias GLclampf                                                    Cfloat
typealias GLbitfield                                                  Cuint
typealias uint_least8_t                                               Uint8
typealias GLhandleARB                                                 Uint32
typealias GLeglImageOES                                               Ptr
typealias Pointer                                                     Ptr{Void}
typealias GLdouble                                                    Cdouble
typealias int_least8_t                                                Uint8
typealias uint64_t                                                    Uint64
typealias GLushort                                                    Cushort
typealias uint8_t                                                     Uint8
typealias int_fast32_t                                                Int64
typealias intmax_t                                                    Int64
typealias wchar_t                                                     Int32
typealias GLhalfARB                                                   Cushort
typealias GLhalf                                                      Cushort
typealias intptr_t                                                    Int64
typealias GLsizei                                                     Csize_t
typealias GLuint64                                                    Culonglong
typealias uintmax_t                                                   Uint64
typealias uint_fast16_t                                               Uint64
typealias int32_t                                                     Int32
typealias __gwchar_t                                                  Int32
typealias GLintptrARB                                                 ptrdiff_t
typealias GLhalfNV                                                    Cushort
typealias GLshort                                                     Cshort
typealias GLint64                                                     Clonglong
typealias int_least16_t                                               Int16
typealias GLintptr                                                    Cptrdiff_t
typealias uint_fast64_t                                               Uint64
typealias GLsync                                                      Ptr{Void}
typealias GLsizeiptrARB                                               ptrdiff_t
typealias GLint64EXT                                                  int64_t
typealias GLuint64EXT                                                 uint64_t
typealias GLvdpauSurfaceNV                                            GLintptr
export uint_least64_t, uint_least16_t, uint_fast32_t, int_fast8_t, int_fast16_t, int_least64_t, GLubyte, int16_t, GLfloat, int8_t, GLenum, GLboolean, uint16_t, GLclampd, GLvoid, uint_fast8_t, int64_t, int_fast64_t, uint32_t, GLint, uint_least32_t, GLsizeiptr, ptrdiff_t, int_least32_t, GLcharARB, GLuint, GLbyte, size_t, GLchar, uintptr_t, GLclampf, GLbitfield, uint_least8_t, GLhandleARB, GLeglImageOES, Pointer, GLdouble, int_least8_t, uint64_t, GLushort, uint8_t, int_fast32_t, intmax_t, wchar_t, GLhalfARB, GLhalf, intptr_t, GLsizei, GLuint64, uintmax_t, uint_fast16_t, int32_t, __gwchar_t, GLintptrARB, GLhalfNV, GLshort, GLint64, int_least16_t, GLintptr, uint_fast64_t, GLsync, GLsizeiptrARB, GLint64EXT, GLuint64EXT, GLvdpauSurfaceNV


text:
# Attempting to use finitediff instead of autodiff for linearizations
using FiniteDifferences
x0, u0 = getXU_0(j);

m = central_fdm(10,1)
A1 = FiniteDifferences.jacobian(m, z->fxdot(z,u0,j,g),x0)[1]
B1 = FiniteDifferences.jacobian(m, z->fxdot(x0,z,j,g),u0)[1]

A2, B2 = linearizeDiff(x0,u0,j)
##
println()
println("Aerr = ", norm(A1-A2))

println("Berr = ", norm(B1 - B2))

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "SPIRV_Tools"
version = v"2020.2"

# Collection of sources required to build IGC
sources = [
    GitSource("https://github.com/KhronosGroup/SPIRV-Tools.git", "fd8e130510a6b002b28eee5885a9505040a9bdc9"),
    # vendored dependencies
    GitSource("https://github.com/google/effcee.git", "6fa2a03cebb4fb18fbad086d53d1054928bef54e"), # 2019.0
    GitSource("https://github.com/google/re2.git", "209eda1b607909cf3c9ad084264039546155aeaa"), # 2020-04-01
]

# Bash recipe for building across all platforms
script = raw"""
# put vendored dependencies in places they will be picked up by the build system
mv effcee SPIRV-Tools/external/effcee
mv re2 SPIRV-Tools/external/re2

cd SPIRV-Tools
install_license LICENSE

CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Explicitly use our cmake toolchain file
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})

# Skip tests
CMAKE_FLAGS+=(-DSPIRV_SKIP_TESTS=ON)

# Don't use -Werror
CMAKE_FLAGS+=(-DSPIRV_WERROR=OFF)

# Point to the SPIRV Headers JLL
CMAKE_FLAGS+=(-DSPIRV-Headers_SOURCE_DIR=${prefix})

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("spirv-as", :spirv_as),
    ExecutableProduct("spirv-cfg", :spirv_cfg),
    ExecutableProduct("spirv-dis", :spirv_dis),
    ExecutableProduct("spirv-link", :spirv_link),
    ExecutableProduct("spirv-opt", :spirv_opt),
    ExecutableProduct("spirv-reduce", :spirv_reduce),
    ExecutableProduct("spirv-val", :spirv_val),
    LibraryProduct("libSPIRV-Tools-shared", :libSPIRV_Tools),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("SPIRV_Headers_jll")
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)