text:
type Class_agg_ls <: abstract_ls_info
    step_size_P::Float64
    step_size_D::Float64
    num_steps::Int64
    frac_progress::Float64
    predict_red::Float64
    actual_red::Float64
    do_ls::Bool

    function Class_agg_ls(iter::Class_iterate, dir::Class_point, pars::Class_parameters)
        this = new()

        this.step_size_D = 0.0
        this.step_size_P = 0.0
        this.num_steps = 0
        this.predict_red = -1.0
        this.frac_progress = NaN
        this.actual_red = NaN
        #mu_ = get_mu(iter) - dir.mu
        #this.do_ls = 0.1 < minimum(comp_predicted_scaled(iter,dir,1.0)) && maximum(comp_predicted(iter,dir,1.0)) .< 10 * get_mu(iter)
        #comp_max = (get_mu(iter) - dir.mu) / pars.ls.comp_feas_agg + get_mu(iter) / 4.0
        #this.do_ls = true
        eta = - dir.mu / get_mu(iter)
        gamma = 1.0 - eta
        y_tilde = (gamma * get_mu(iter) - eta * iter.point.y .* get_primal_res(iter)) ./ iter.point.s
        this.do_ls = dot(eval_grad_lag(iter, get_mu(iter) * gamma, y_tilde), dir.x) < 0.0

        #this.do_ls = true
        #norm(comp_predicted(iter,dir,1.0),Inf) < get_mu(iter) / pars.ls.comp_feas_agg
        return this
    end
end

function accept_func!(accept::Class_agg_ls, intial_it::Class_iterate, candidate::Class_iterate, dir::Class_point, step_size::Float64, filter::Array{Class_filter,1},  pars::Class_parameters, timer::class_advanced_timer)
    tau = get_mu(candidate) / (scaled_dual_feas(candidate, pars) * (1.0 - pars.agg_protection_factor))
    max_thres = pars.agg_protection_factor
    #@show get_mu(candidate), get_mu(intial_it)
    #@show get_mu(candidate), get_mu(intial_it)
    if get_mu(candidate) / get_mu(intial_it) >= 1.0 - max_thres || tau >= 1.0
      #@show get_mu(candidate) / get_mu(intial_it), tau
      return :success, step_size
    else
      #println("failure")
      suggested_step_size = max(max_thres^2, accept.step_size_P * tau^2)
      return :not_enough_progress, suggested_step_size
    end
end


text:
using Test, Mineos

@testset "All tests" begin
    include("eigenfrequencies.jl")
end

ext: jl


text:
using Random, Printf, Colors

import Luxor
const lx = Luxor

function triangle(points, degree)
    lx.sethue(cols[degree])
    lx.poly(points, :fill)
end

function sierpinski(points, degree)
    triangle(points, degree)
    if degree > 1
        p1, p2, p3 = points
        sierpinski([p1, lx.midpoint(p1, p2),
                        lx.midpoint(p1, p3)], degree-1)
        sierpinski([p2, lx.midpoint(p1, p2),
                        lx.midpoint(p2, p3)], degree-1)
        sierpinski([p3, lx.midpoint(p3, p2),
                        lx.midpoint(p1, p3)], degree-1)
    end
end

function draw(n)
    lx.circle(lx.O, 100, :clip)
    points = lx.ngon(lx.O, 150, 3, -pi/2, vertices=true)
    sierpinski(points, n)
end

lx.Drawing(400, 250)
lx.background("white")
lx.origin()

depth = 8
cols = distinguishable_colors(depth) # from Colors.jl
draw(depth)

lx.finish()
lx.preview()

import Winston
const wn = Winston

t = collect(range(0,stop=4pi,length=1000));
f(x::Array) = 10x .* exp.(-0.3x) .* sin.(3x);
g(x::Array) = 0.1x.*(2pi .- x).*(4pi .- x);
h(x::Array) = 10.0 ./ (1 .+ x.*x);  # Note need spaces tp clarify 1./ and 1 ./

y1 = f(t)
y2 = g(t)
y3 = h(t)
wn.plot(t,y1,"b",t,y2,"r--",t,y3,"k;")

x = randn(10000)
wn.plothist(x, nbins=50)

p = wn.FramedPlot( title="Frequency vs Time", ylabel="\\Theta",xlabel="\\Tau");
wn.add(p,wn.FillBetween(t,y1,t,y2));
wn.add(p,wn.Curve(t,y1,color="red"));
wn.add(p,wn.Curve(t,y2,color="blue"));

wn.display(p)

wn.semilogx(t,y1)
wn.title("log(t) vs 10x * exp(-0.3x) * sin(3x)")

wn.semilogy(y2,y3)
wn.title("0.1*(2\\pi - x)*(4\\pi - x)  vs  log(1 /( 1 + x*x))")

p = wn.FramedPlot(aspect_ratio=1,xrange=(-10,110),yrange=(-10,110));
n = 21;
x = collect(range(0.0, length=n, stop=100.0));

# Create a set of random variates

yA =  10.0*randn(n) .+ 40.0;
yB =  x .+ 5.0*randn(n);

# Set labels and symbol styles

a = wn.Points(x, yA, kind="circle");
wn.setattr(a,label="'a' points");

b = wn.Points(x, yB);
wn.setattr(b,label="'b' points");
wn.style(b, kind="filled circle");

# Plot a line which 'fits' through the yB points
# and add a legend in the top LHS part of the graph

s = wn.Slope(1, (0,0), kind="dotted");
wn.setattr(s, label="slope");
lg = wn.Legend(.1, .9, Any[a,b,s] );
wn.add(p, s, a, b, lg);

wn.display(p)

wn.savefig(p, "MyWPlot.png")

import Cairo
const ca = Cairo

c  = ca.CairoRGBSurface(512, 128);
cr = ca.CairoContext(c);
ca.save(cr);

ca.set_source_rgb(cr, 0.8, 0.8, 0.8);
ca.rectangle(cr, 0.0, 0.0, 512.0, 128.0); # background
ca.fill(cr);
ca.restore(cr);
ca.save(cr);

x0=51.2;  y0=64.0;
x1=204.8; y1=115.4;
x2=307.2; y2=12.8;
x3=460.8; y3=64.0;

ca.move_to(cr, x0, y0);
ca.curve_to(cr, x1, y1, x2, y2, x3, y3);
ca.set_line_width(cr, 10.0);
ca.stroke(cr);
ca.restore(cr);

ca.move_to(cr, 12.0, 12.0);
ca.set_source_rgb(cr, 0, 0, 0);
ca.show_text(cr,"Cairo Curve")
ca.write_to_png(c,"cstroke.png");

run(`display cstroke.png`)

import PyPlot
const py = PyPlot

x = collect(range(0.0,stop=2pi,length=1000))
y = sin.(3*x + 4*cos.(2*x));
py.title("A sinusoidally modulated sinusoid");
py.plot(x, y, color="red", linewidth=2.0, linestyle="--");

y = collect(range(0,stop=3π,length=250))
py.surf(y, y, y .* sin.(y) .* cos.(y)' .* exp.(-0.4y))

py.xkcd()
x = collect(range(1, length=101, stop=10));
y = sin.(3x + cos.(5x))
p = py.plot(x,y)

import Gadfly
gd = Gadfly

dd  =  gd.plot(x =  rand(10),  y = rand(10));
gd.draw(SVG("random-pts.svg",  15cm, 12cm) , dd);

x = collect(1:100);
y1 = ones(100) - 2*rand(100);
y2 = randn(100);

gd.plot(
 gd.layer(x=x,y=y1,gd.Geom.line,gd.Theme(default_color=gd.colorant"red")),
 gd.layer(x=x,y=y2,gd.Geom.line,gd.Theme(default_color=gd.colorant"blue"))
)

using RDatasets, DataFrames;

mlmf = dataset("mlmRev","Gcsemv")
df = mlmf[completecases(mlmf), :]

names(df)

# gd.set_default_plot_size(20cm, 12cm);
gd.plot(df, x="Course", y="Written", color="Gender")

gd.plot((x,y) -> exp.(-x) .* cos.(y), -2, 0.4, 0, 4)

gd.plot((x,y) -> x .* exp.(-(x - floor.(x))).^2 .- y.^2, -8.0, 8, -2.0, 2.0)

import Compose
const cm =  Compose

function sierpinski(n)
  if n == 0
    cm.compose(cm.context(), cm.polygon([(1,1), (0,1), (1/2, 0)]));
  else
    t = sierpinski(n - 1);
    cm.compose(cm.context(), (cm.context( 1/4, 0, 1/2, 1/2), t),
                             (cm.context( 0, 1/2, 1/2, 1/2), t),
                             (cm.context( 1/2, 1/2, 1/2, 1/2), t));
  end
end

cx1 = compose(sierpinski(1), linewidth(0.2mm),fill(nothing), stroke("black"));
img = SVG("sierp1.svg", 10cm, 8.66cm); draw(img,cx1);

cx3 = compose(sierpinski(3), linewidth(0.2mm),fill(nothing), stroke("black"));
img = SVG("sierp3.svg", 10cm, 8.66cm); draw(img,cx3);

cx5 = compose(sierpinski(5), linewidth(0.2mm),fill(nothing), stroke("black"));
img = SVG("sierp5.svg", 10cm, 8.66cm); draw(img,cx5);

using PGFPlots

p = Axis( [ Plots.Linear(x -> sin.(3x) .* exp.(-0.3x), (0,8),
                         legendentry = L"$\sin(3x)*exp(-0.3x)$"),
            Plots.Linear(x -> sqrt.(x) ./ (1+x.^2), (0,8),
                         legendentry = L"$\sqrt{2x}/(1+x^2)$")
        ])

# This requires installation of the pdf2svg utility
save("linear-plots.svg", p);

fq  =  randn(10000);
p = Axis(Plots.Histogram(fq, bins=100), ymin=0)

save("histogram-plot.svg", p);


text:
#
# These types define the different connection states for the protocol.
#

# TODO: REname to ConnectionState, because the socket is the TCP Socket
"Enum value for the different states a WebSocket can be in."
struct SocketState
	v::Symbol
end

const STATE_CONNECTING     = SocketState(:connecting)
const STATE_OPEN           = SocketState(:open)
const STATE_CLOSING        = SocketState(:closing)
const STATE_CLOSED         = SocketState(:closed)

#
# These types define the input interface for the client logic.
#

"Abstract type for all commands sent to `ClientProtocol`.

These commands are sent as arguments to the different `handle` functions on `ClientProtocol`. Each
command represents an action on a WebSocket, such as sending a text frame, ping request, or closing
the connection."
abstract type ClientProtocolInput end

"Send a text frame, sent to `ClientProtocol`."
struct SendTextFrame <: ClientProtocolInput
	data::Vector{UInt8}
	# True if this is the final frame in the text message.
	isfinal::Bool
	# What WebSocket opcode should be used.
	opcode::Opcode

	SendTextFrame(data::Vector{UInt8}, isfinal::Bool, opcode::Opcode) = new(data, isfinal, opcode)
	SendTextFrame(data::String, isfinal::Bool, opcode::Opcode) = new(unsafe_wrap(Vector{UInt8}, data), isfinal, opcode)
end

"Send a binary frame, sent to `ClientProtocol`."
struct SendBinaryFrame <: ClientProtocolInput
	data::Array{UInt8, 1}
	# True if this is the final frame in the text message.
	isfinal::Bool
	# What WebSocket opcode should be used.
	opcode::Opcode
end

"Send a ping request to the server."
struct ClientPingRequest  <: ClientProtocolInput end

"A frame was received from the server."
struct FrameFromServer <: ClientProtocolInput
	frame::Frame
end

"A request to close the WebSocket."
struct CloseRequest <: ClientProtocolInput end

"Used when the underlying network socket was closed."
struct SocketClosed <: ClientProtocolInput end

"A pong reply was expected, but never received."
struct PongMissed <: ClientProtocolInput end

"No Close frame response was received by the client after a reasonable time"
struct AbnormalNoCloseResponseReceived <: ClientProtocolInput end

"The socket did not close the socket in a reasonable timeframe."
struct AbnormalSocketNotClosedByServer <: ClientProtocolInput end


text:
module SpringCollab2020LeaveTheoBehindTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/LeaveTheoBehindTrigger" LeaveTheoBehindTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight)

const placements = Ahorn.PlacementDict(
    "Leave Theo Behind (Spring Collab 2020)" => Ahorn.EntityPlacement(
        LeaveTheoBehindTrigger,
        "rectangle",
    ),
)

end


text:
@testset "check gradient for Linear block " begin
    clear()
    # [1] prepare input data and its label
    T = Array{Float64}
    x = randn(256, 62)
    l = rand(64, 62)
    l = l ./ sum(l,dims=1)
    m = Linear(256, 64; type=T)

    x = Variable(x; type=T);
    l = Variable(l; type=T);

    # [2] forward and backward propagation
    outs = forward(m, x)
    LOSS1 = mse(outs, l)
    COST1 = cost(LOSS1)
    backward()

    # [3] with a samll change of a weight
    GRAD = m.w.delta[1]
    DELTA = 1e-6
    m.w.value[1] += DELTA

    # [4] forward and backward propagation
    outs = forward(m, x)
    LOSS2 = mse(outs, l)
    COST2 = cost(LOSS2)
    backward()

    # [5] check if the auto-grad is true or not
    dLdW = (COST2 - COST1)/DELTA;   # numerical gradient
    err  = abs((dLdW-GRAD)/(GRAD+eps(Float64)))*100;  # relative error in %
    err  = err < 1e-1 ? 0.0 : err;
    @test err<1.0
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using FEASTSolver
using LinearAlgebra
using Printf
using Random
using MatrixMarket
using NonlinearEigenproblems: nep_gallery, compute_Mder
using DelimitedFiles

function bf()
    N = diagm(-1 => [1,1,1,1,1,1,1])
    Mh0 = 1/6*(4*I + N + N')
    Mh1 = N - N'
    Mh2 = -1*(2*I - N - N')
    Mh3 = Mh1
    Mh4 = -Mh2
    c = [0.6 1.3; 1.3 0.1; 0.1 1.2; 1.0 1.0; 1.2 1.0]
    I8 = Matrix(I, 8, 8)
    M0 = c[1,1] * kron(I8, Mh0) + c[1,2] * kron(Mh0, I8)
    M1 = c[2,1] * kron(I8, Mh1) + c[2,2] * kron(Mh1, I8)
    M2 = c[3,1] * kron(I8, Mh2) + c[3,2] * kron(Mh2, I8)
    M3 = c[4,1] * kron(I8, Mh3) + c[4,2] * kron(Mh3, I8)
    M4 = c[5,1] * kron(I8, Mh4) + c[5,2] * kron(Mh4, I8)
    function butterfly(A0, A1, A2, A3, A4)
        f(z) =  z^4*A4 + z^3*A3 + z^2*A2 + z*A1 + A0
    end

    return butterfly(M0, M1, M2, M3, M4)
end

function deficient()
	# Random.seed!(1234)
	# A0 = complex.(rand(15,15))
	# A1 = complex.(rand(15,15))
	# A0[:,1] .= 0
	A0 = mmread("data/quadraticM0.mtx")
	A1 = mmread("data/quadraticM1.mtx")

	function T(z::ComplexF64)
		return (z + 0.2)*(z - 0.1).*A1 .+ A0
	end

	return T
end

function quadratic()
	A0 = mmread("data/system5A0.mtx")
	A1 = mmread("data/system5A1.mtx")
	A2 = mmread("data/system5A2.mtx")

	function T(z::ComplexF64)
		return (A2 .* (z*z)) .+ (A1 .* z) .+ A0
	end

	return T
end

function gun()
	nep_gun = nep_gallery("nlevp_native_gun")
	T(z) = compute_Mder(nep_gun, z)
	return T
end

function hadeler()
	nep_hadeler = nep_gallery("nlevp_native_hadeler", 100, 200)
	T(z) = compute_Mder(nep_hadeler, z)
	return T
end

function info(Λ, X, residuals, c, r)
    inside(x) = in_contour(x, c, r)
    in_eig = Λ[inside.(Λ)]
    in_res = residuals[inside.(Λ)]
    print("number inside : ")
    println(size(Λ[inside.(Λ)])[1])
    if sum(inside.(Λ)) > 0
        in_res_conv = in_res[in_res .<= 1e-3]
        in_eig_conv = in_eig[in_res .<= 1e-3]
        print("number inside converged : ")
        println(size(in_eig_conv)[1])
        print("max res inside: ")
        println(maximum(residuals[inside.(Λ)]))
        if size(in_res_conv, 1) > 0
            print("max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
    # display(in_eig_conv)
end

function res_table(feast_call, dim, max_iter, max_nodes, spurious, c, r, err_tol)
	all_res = zeros((max_iter, max_nodes))
	println(repeat("-", max_nodes*8*2))
	for j=1:max_nodes
	    print(2^(j+1))
	    print("\t\t")
	end
	println("\n"*repeat("-", max_nodes*8*2))
	for i=1:max_iter
	    for j=1:max_nodes
			if i>1 && minimum(all_res[1:i-1,j]) < err_tol
				# all_res[i,j] = all_res[i-1,j]
				all_res[i,j] = -Inf
				print("\t\t")
			else
		        inside(z) = in_contour(z, c, r)
		        e, v, res = feast_call(rand(ComplexF64, dim), 2^(j+1), i-1, err_tol)
		        in_eig = e[inside.(e)]
		        in_res = res[inside.(e)]
		        in_res_conv = in_res[in_res .<= spurious]
		        if size(in_res_conv, 1) == 0
		            all_res[i,j] = maximum(res[inside.(e)])
		        else
		            all_res[i,j] = maximum(in_res_conv)
		        end
		        @printf("%2.2e\t", all_res[i,j])
			end
	    end
	    println()
	end
	println(repeat("-", max_nodes*8*2))
	return all_res
end

err_tol_all = 10e-16
max_iter_all = 10

# for i=2:8
# 	println("\nBeyn N=$(2^i)")
# 	e, v, res = nlfeast!(bf(), rand(ComplexF64,64,30), 2^i, 0, c=complex(1.0,1.0), r=0.5)
# 	info(e, v, res, complex(1.0,1.0), 0.5)
# end

Random.seed!(1729)

# println("\nButterfly Problem")
# max_iter_butterfly = 25
# bf_feast_call(X, N, iter, err_tol) = nlfeast!(bf(), X, N, iter, c=complex(1.0,1.0), r=0.5, ϵ=err_tol)
# bf_res = res_table(bf_feast_call, (64,20), max_iter_butterfly, 7, 1e-0, complex(1.0,1.0), 0.5, err_tol_all)
# writedlm("test/data/bf.dat", zip(1:max_iter_butterfly, eachcol(bf_res)...))
#
#
# println("\nDeficient Quadratic Problem")
#
# deficient_feast_call(X, N, iter, err_tol) = nlfeast_moments!(deficient(), X, N, iter, c=complex(0.0,0.0), r=0.25, ϵ=err_tol, moments=2, store=true)
# def_res = res_table(deficient_feast_call, (15,4), max_iter_all, 7, 1e-5, complex(0.0,0.0), 0.25, err_tol_all)
# writedlm("test/data/deficient.dat", zip(1:max_iter_all, eachcol(def_res)...))
#
#
# println("\nLarger Quadratic Problem")
# quad_feast_call(X, N, iter, err_tol) = nlfeast_moments!(quadratic(), X, N, iter, c=complex(-1.55,0.0), r=0.05, ϵ=err_tol, moments=2, store=true)
# quad_res = res_table(quad_feast_call, (1000,70), max_iter_all, 6, 1e-6, complex(-1.55,0.0), 0.05, err_tol_all)
# writedlm("test/data/quad.dat", zip(1:max_iter_all, eachcol(quad_res)...))
#
#
# println("\nGun Cavity Problem")
# gun_feast_call(X, N, iter, err_tol) = nlfeast!(gun(), X, N, iter, c=complex(140000.0, 0.0), r=30000, ϵ=err_tol, store=false)
# gun_res = res_table(gun_feast_call, (9956,32), max_iter_all, 6, 1e-6, complex(140000.0, 0.0), 30000, err_tol_all)
# writedlm("test/data/gun.dat", zip(1:max_iter_all, eachcol(gun_res)...))

println("\nHadeler Problem")
hadeler_feast_call(X, N, iter, err_tol) = nlfeast!(hadeler(), X, N, iter, c=complex(-30.0, 0.0), r=10.0, ϵ=err_tol, store=true)
hadeler_res = res_table(hadeler_feast_call, (200,15), max_iter_all, 6, 1e-4, complex(-30.0, 0.0), 10.0, err_tol_all)
writedlm("test/data/hadeler.dat", zip(1:max_iter_all, eachcol(hadeler_res)...))


text:
export   LinearTransform,
         MinMaxTransform,
         updateLinearTransform,
         transform!,
         transform,
         itransform!,
         itransform

 """
 $(TYPEDEF)

 This type holds an affine transformation of the form x-> μ + Lx, where μ is an offset and L is a linear operator.
 We typically use this type to normalize samples to have zero mean and unit covariance of element-wise standard deviation.

 ## Fields

 $(TYPEDFIELDS)

 """
struct LinearTransform
    Nx::Int64

    # Flag for using diagonal or dense transformation
    μ::Array{Float64,1}

    L::Union{Diagonal, LowerTriangular}

    diag::Bool
end

# Identity transformation
LinearTransform(Nx::Int64) = LinearTransform(Nx, zeros(Nx), Diagonal(ones(Nx)), true)

function LinearTransform(X::AbstractMatrix{Float64}; diag::Bool=true, factor::Float64=1.0)
    # The entry `factor` allows to set the entry-wise standard deviation to a value different from 1
    Nx, Ne = size(X)
    μ = deepcopy(mean(X; dims = 2)[:,1])

    if diag == true || Nx == 1
        σ = std(X; dims = 2)[:,1]
        L = deepcopy((1.0/factor)*Diagonal(σ))
        diag = true
    else #create a dense transformation from the Cholesky factor
        @assert Nx>1 "Only works for Nx>1, otherwise use first method"
        L = deepcopy((1.0/factor)*cholesky(cov(X')).L)
    end

    return LinearTransform(Nx, μ, L, diag)
end

function updateLinearTransform!(Lin::LinearTransform, X::AbstractMatrix{Float64}; diag::Bool=true)
    Nx, Ne = size(X)
    Lin.μ .= copy(mean(X; dims = 2)[:,1])

    if diag == true || Nx == 1
        σ = copy(std(X; dims = 2)[:,1])
        Lin.L .= Diagonal(σ)

    else #create a dense transformation from the Cholesky factor
        @assert Nx>1 "Only works for Nx>1, otherwise use diagonal scaling"
        Lin.L = copy(cholesky(cov(X')).L)
    end
end

# For scalar vectors
function transform!(L::LinearTransform, xout::AbstractVector{Float64}, xin::AbstractVector{Float64})
    @assert size(xout,1) == size(xin,1) "Input and output dimensions don't match"
    @assert size(xin,1) == L.Nx "Input dimension is incorrect"

    copy!(xout, xin)
    xout -= L.μ

    ldiv!(L.L, xout)
end

function transform!(L::LinearTransform, x::AbstractVector{Float64})
    @assert size(x,1) == L.Nx "Input dimension is incorrect"
    x -= L.μ
    ldiv!(L.L, x)
    return x
end

transform(L::LinearTransform, x::AbstractVector{Float64}) = transform!(L, zero(x), x)

function itransform!(L::LinearTransform, xout::AbstractVector{Float64}, xin::AbstractVector{Float64})
    @assert size(xout) == size(xin) "Input and output dimensions don't match"
    @assert size(xin,1) == L.Nx "Input dimension is incorrect"

    copy!(xout, xin)
    mul!(xout, L.L, xout)
    xout += L.μ

    return xout
end

function itransform!(L::LinearTransform, x::AbstractVector{Float64}, idx::Union{UnitRange{Int64}, Array{Int64,1}})
    @assert size(x,1) == L.Nx "Input dimension is incorrect"

    if typeof(L.L)<:Diagonal
        mul!(x, Diagonal(view(L.L.diag,idx)), x)
    else
        error("Not yet implemented")
    end

    # mul!(X, L.L, X)
    x += L.μ

    return x
end

function itransform!(L::LinearTransform, x::AbstractVector{Float64})
    @assert size(x,1) == L.Nx "Input dimension is incorrect"

    mul!(x, L.L, x)
    x += L.μ

    return x
end

itransform(L::LinearTransform, x::AbstractVector{Float64}) = itransform!(L, zero(x), x)

# For ensemble matrix
function transform!(L::LinearTransform, Xout::AbstractArray{Float64}, Xin::AbstractMatrix{Float64})
    @assert size(Xout,1) == size(Xin,1) "Input and output dimensions don't match"
    @assert size(Xin,1) == L.Nx "Input dimension is incorrect"

    copy!(Xout, Xin)
    Xout .-= L.μ

    ldiv!(L.L, Xout)
    # return Xout
end

function transform!(L::LinearTransform, X::AbstractMatrix{Float64}, idx::Union{UnitRange{Int64}, Array{Int64,1}})
    @assert size(X,1) == length(idx) "Input dimension is incorrect"
    X .-= view(L.μ,idx)
    if typeof(L.L)<:Diagonal
        ldiv!(Diagonal(view(L.L.diag,idx)), X)
    else
        error("Not yet implemented")
    end
    return X
end

function transform!(L::LinearTransform, X::AbstractMatrix{Float64})
    @assert size(X,1) == L.Nx "Input dimension is incorrect"
    X .-= L.μ
    ldiv!(L.L, X)
    return X
end

transform(L::LinearTransform, X::AbstractMatrix{Float64}) = transform!(L, zero(X), X)

transform(X::AbstractMatrix{Float64}; diag::Bool = true) = transform!(LinearTransform(X; diag = diag), zero(X), X)

function itransform!(L::LinearTransform, Xout::AbstractMatrix{Float64}, Xin::AbstractMatrix{Float64})
    @assert size(Xout) == size(Xin) "Input and output dimensions don't match"
    @assert size(Xin,1) == L.Nx "Input dimension is incorrect"

    copy!(Xout, Xin)
    mul!(Xout, L.L, Xout)
    Xout .+= L.μ

    return Xout
end

function itransform!(L::LinearTransform, X::AbstractMatrix{Float64}, idx::Union{UnitRange{Int64}, Array{Int64,1}})
    @assert size(X,1) == L.Nx "Input dimension is incorrect"

    if typeof(L.L)<:Diagonal
        mul!(X, Diagonal(view(L.L.diag,idx)), X)
    else
        error("Not yet implemented")
    end

    # mul!(X, L.L, X)
    X .+= L.μ

    return X
end

function itransform!(L::LinearTransform, X::AbstractMatrix{Float64})
    @assert size(X,1) == L.Nx "Input dimension is incorrect"

    mul!(X, L.L, X)
    X .+= L.μ

    return X
end

itransform(L::LinearTransform, X::AbstractMatrix{Float64}) = itransform!(L, zero(X), X)

# Define a min-max transformation to use with the package c3
struct MinMaxTransform
    Nx::Int64

    L::Diagonal #(maximum(X[i,:] - minimum(X[i,:]))

    b::Array{Float64,1} #minimum(X[i,:])/(maximum(X[i,:] - minimum(X[i,:]))
end

function MinMaxTransform(X::AbstractMatrix{Float64})
    Nx = size(X,1)
     b = zeros(Nx)
     L = Diagonal(zeros(Nx))

     @inbounds for i=1:Nx
         xi = view(X,i,:)
         Mi = maximum(xi)
         mi = minimum(xi)
         @assert Mi>mi "Minimum and maximum are equal"
         L.diag[i] = (Mi - mi)/2
         b[i] = (Mi + mi)/(Mi - mi)
     end
    return MinMaxTransform(Nx, L, b)
end

function transform!(L::MinMaxTransform, Xout::AbstractMatrix{Float64}, Xin::AbstractMatrix{Float64})
    @assert size(Xout,1) == size(Xin,1) "Input and output dimensions don't match"
    @assert size(Xin,1) == L.Nx "Input dimension is incorrect"

    copy!(Xout, Xin)
    ldiv!(L.L, Xout)
    Xout .-= L.b

    # return Xout
end

function transform!(L::MinMaxTransform, X::AbstractMatrix{Float64})
    @assert size(X,1) == L.Nx "Input dimension is incorrect"
    ldiv!(L.L, X)
    X .-= L.b
    return X
end

transform(L::MinMaxTransform, X::AbstractMatrix{Float64}) = transform!(L, zero(X), X)

function itransform!(L::MinMaxTransform, Xout::AbstractMatrix{Float64}, Xin::AbstractMatrix{Float64})
    @assert size(Xout,1) == size(Xin,1) "Input and output dimensions don't match"
    @assert size(Xin,1) == L.Nx "Input dimension is incorrect"

    copy!(Xout, Xin)
    Xout .+= L.b
    mul!(Xout, L.L, Xout)
    return Xout
end

function itransform!(L::MinMaxTransform, X::AbstractMatrix{Float64})
    X .+= L.b
    mul!(X, L.L, X)
    return X
end

itransform(L::MinMaxTransform, X::AbstractMatrix{Float64}) = itransform!(L, zero(X), X)


text:
using Sobol
using NLopt
using PyPlot
using SpecialFunctions
using SparseArrays
using LinearAlgebra
using BenchmarkTools
using PyPlot
using Statistics
using MultistartOptimization

### ALL FUNCTIONS NEEDED FOR SOBOL PT GENERATION AND MLSL

# Generate an N-dimentional Sobol sequence of points
function gen_sobol(Npts, Dims,lb,ub)
    sobol_seq = SobolSeq(lb, ub)
    int_arr = [reshape(next!(sobol_seq),1,Dims) for i = 1:Npts]
    sobol_array = permutedims(reshape(hcat(int_arr...),(length(int_arr[1]), length(int_arr))))

    return sobol_array
end

# Generate 9-Dimentional schwefel test function:
# Global min is: f(x) = 0 at x = (420.9687,...,420.9687)
function schwefel_funcND(x)
    sum = 0
    Dims = length(x)
    for ii = 1:Dims
        xi = x[ii]
        sum = sum + xi*sin(sqrt(abs(xi)))
    end
    y = 418.982887*Dims - sum
    return y
end

# derivative of 9-Dimentional schwefel test function:
function schwefel_funcND_der(x)
    hold = similar(x)
    for ii = 1:size(hold)[1]
        hold[ii] = -sin(sqrt(abs(x[ii]))) - (x[ii]*cos(sqrt(abs(x[ii])))*sign(x[ii])/(2*sqrt(abs(x[ii]))))
    end
    return hold
end

#  test function for 9d schwefel
function testfunc_Nd_shwef(x::Vector, grad::Vector)
    if length(grad) > 0
        for i = 1:length(grad)
            grad[i] = schwefel_funcND_der(x)[i]
        end
    end
    return schwefel_funcND(x)
end

# camel function for testing
function six_hump_camel(x)
    y = 4*x[1]^2 - 2.1*x[1]^4 + (1/3)x[1]^6 + x[1]*x[2] - 4*x[2]^2+ 4*x[2]^4
    return y
end

# 2D shubert function
function shubert_func(x)

    x1 = x[1]
    x2 = x[2]
    sum1 = 0
    sum2 = 0
    for ii = 1:5
        new1 = ii * cos((ii+1)*x1+ii)
        new2 = ii * cos((ii+1)*x2+ii)
        sum1 = sum1 + new1
        sum2 = sum2 + new2
    end
    return sum1 * sum2

end

# derivative of 2D shubert function
function shubert_func_deriv(x)
    x1 = x[1]
    x2 = x[2]
    sum1 = 0
    sum2 = 0
    for ii = 1:5
        new1 = ii * cos((ii+1)*x1+ii)
        sum1 = sum1 + new1
    end
    for jj = 1:5
        new2 = -(jj * (jj+1)* sin((jj+1)*x2+jj))
        sum2 = sum2 + new2
    end
    der2 = sum1*sum2
    sum1 = 0
    sum2 = 0
    for ii = 1:5
        new1 = ii * cos((ii+1)*x2+ii)
        sum1 = sum1 + new1
    end
    for jj = 1:5
        new2 = -(jj * (jj+1)* sin((jj+1)*x1+jj))
        sum2 = sum2 + new2
    end

    der1 = sum1*sum2

    return [der1 der2]
end

#  test function for 2d shubert
function testfunc_shub(x::Vector, grad::Vector)
    if length(grad) > 0
        # gradients for six hump camel
        grad[1] = shubert_func_deriv(x)[1]
        grad[2] = shubert_func_deriv(x)[1]
    end
    return shubert_func(x)
end

# tes function for 2d six hump camel
function testfunc_6hc(x::Vector, grad::Vector)
    if length(grad) > 0
        # gradients for six hump camel
        grad[1] = 8*x[1] - 8.4*x[1]^3 + 2*x[1]^5 + x[2]
        grad[2] = x[1] - 8*x[2] + 16*x[2]^3
    end
    return six_hump_camel(x)
end

## optimizers for local phase of MLSL
function local_search(x_init,lb,ub,dims,loc_alg,tol,testfunc)

    opt = Opt(loc_alg, dims)
    opt.lower_bounds = lb
    opt.upper_bounds = ub
    opt.xtol_rel = tol
    opt.maxeval = 1000

    opt.min_objective = testfunc

    (minf,minx,ret) = optimize(opt, x_init)
    numevals = opt.numevals # the number of function evaluations

    return minf,minx,numevals
    #println("got $minf at $minx after $numevals iterations (returned $ret)")

end

## test global optimizers from Nlopt
function glob_test(lb,ub,dims,loc_alg,tol,obj_func,glob_alg,maxvl)

    optloc = Opt(loc_alg, dims)
    optloc.lower_bounds = lb
    optloc.upper_bounds = ub
    optloc.xtol_rel = tol
    optloc.maxeval = maxvl

    optloc.min_objective = obj_func

    opt = Opt(glob_alg, dims)
    opt.lower_bounds = lb
    opt.upper_bounds = ub
    opt.xtol_rel = tol
    opt.maxeval = maxvl

    opt.min_objective = obj_func
    opt.population = 64*dims^2 # stay consistent with my version of MLSL
    opt.local_optimizer = optloc

    (minf,minx,ret) = optimize(opt,ones(dims))
    numevals = opt.numevals # the number of function evaluations

    return minf,minx,numevals

end

# sort the N Sobol points by objective function value and slice alpha*N points
function sort_and_reduce(pts,Npts,func,alpha)
    func_val_pts = Array{Float64, 1}(undef,Npts)
    for iter in 1:Npts
        func_val_pts[iter] = func(pts[iter,:])
    end
    perm = sortperm(func_val_pts)

    permed_pts = copy(pts)
    for iter in 1:Npts
        permed_pts[iter,:] = pts[perm[iter],:]
    end

    return permed_pts[1:round(Int64,alpha*Npts),:], func_val_pts[perm][1:round(Int64,alpha*Npts),:], round(Int64,alpha*Npts)

end

# return the critical distance
function distance_metric(k,n_dims,Nr,sigma,omega_n,lb,ub)
    return (((abs(prod(ub - lb)))/omega_n)*(sigma*log(k*Nr))/(k*Nr))^(1/n_dims)
end

# return TRUE if the euclidean distance between two points is smaller
# return FALSE if greater
function is_close_euclidean(pt1,pt2,crit_dist)
    if norm(pt1 - pt2) >= crit_dist
        return false
    else
        return true
    end
end

# check if this local min has already been found and is it the set arr
function check_if_same_min(arr, f_cand)
    return f_cand in arr
end

# Basically just test a pure random search
function PRS(Npts,Dims,objfunc,lb,ub,loc_alg,tol,testfunc)
    func_evals = 0 # return total number of objective function evaluations
    holder = Array{Float64, 2}(undef,Npts,Dims+1)
    sobol_vect = gen_sobol(Npts, Dims,lb,ub)
    for iter in 1:Npts
        x_init = sobol_vect[iter,:]

        # perform local search on Sobol pt.
        holder[iter,1],holder[iter,2:end],loc_func_evals = local_search(x_init,lb,ub,Dims,loc_alg,tol,testfunc)

        func_evals += loc_func_evals
    end

    glob_min = argmin(holder[:,1])

    return holder[glob_min,1], holder[glob_min,2:end], func_evals
end

function MLSL(N_dims,num_objective_points,bb_lb,bb_ub,alpha_factor,ls_func,iter_MAX,sigma,loc_alg,tol,testfuncc)

    #### Set some important parameters ###
    omega_n = (pi^(N_dims/2))/(gamma(1+(N_dims/2))) # precompute omega factor for crit dist
    ###  ###  ###  ###  ###  ###  ###  ###

    W = 0
    N_evals = 0
    W_exp = 0
    K = 0

    ### generate Sobol objective points
    points_array = gen_sobol(num_objective_points, N_dims,bb_lb,bb_ub)
    reduced_pts,reduced_fvals,reduced_num_pts = sort_and_reduce(points_array,num_objective_points,ls_func,alpha_factor)

    ### generate a (1,Nr) array of cluster assignments and keep updating
    ### an index val of 0 means that no cluster has yet to be assigned
    clusters = zeros(Int8, reduced_num_pts)

    ### initialize a dynamic array of current local minima
    local_mins = Array{Float64}(undef,1)
    local_min_locs = Array{Float64,2}(undef,1,N_dims)

    ### initialize a dynamic array of errors and function evaluation numbers
    # return total number of objective function evaluations:
    # NOTE that 'N_sobol_pts' func vals have already been evaluated during set reduction
    func_evals = num_objective_points

    while (K < iter_MAX) && (W_exp < W + 0.5)

        K += 1
        crit_dist = distance_metric(K,N_dims,reduced_num_pts,sigma,omega_n,bb_lb,bb_ub) # calculate critical dist for iteration

        for i = 1:reduced_num_pts
            if clusters[i] == 0
                if W == 0 ## BASE CASE
                    lm,loc,fevals = local_search(reduced_pts[i,:],bb_lb,bb_ub,N_dims,loc_alg,tol,testfuncc)
                    func_evals += fevals
                    N_evals += 1
                    local_mins[1] = lm
                    local_min_locs[1,:] = loc
                    clusters[i] = 1   # assign first cluster
                    W += 1
                else
                    ### sort through all cluster-assigned points in Nr and see if
                    ### canditate falls in cluster via r_crit.
                    for clus_num = 1:W
                        for (ind, val) in enumerate(findall(x -> x==clus_num,clusters))
                            if (is_close_euclidean(reduced_pts[i,:],reduced_pts[val,:],crit_dist) == true) && (reduced_fvals[i] >= reduced_fvals[val])
                                clusters[i] = clus_num
                                break
                            end
                        end
                    end

                end

                # if no cluster assignment found
                if clusters[i] == 0
                    lm,loc,fevals = local_search(reduced_pts[i,:],bb_lb,bb_ub,N_dims,loc_alg,tol,testfuncc)
                    func_evals += fevals
                    N_evals += 1
                    local_mins = [local_mins;lm]
                    local_min_locs = [local_min_locs;loc']
                    clusters[i] = 1   # assign next cluster
                    if !check_if_same_min(local_mins, lm) # if minimum already in set, so not increment stopval
                        W += 1
                    end
                end

            else
                continue
            end

        end
        W_exp = (W*(N_evals-1))/(N_evals-W-2)
    end

    glob_min = argmin(local_mins)
    return local_mins[glob_min], local_min_locs[glob_min,:], func_evals
end


text:
using Test
using VortexLattice
using LinearAlgebra

ztol = sqrt(eps())

# reflects a vector across the x-z plane
flipy(x) = [x[1], -x[2], x[3]]

# constructs a normal vector the way AVL does
function avl_normal_vector(dr, theta)
    st, ct = sincos(theta)
    bhat = dr/norm(dr) # bound vortex vector
    shat = [0, -dr[3], dr[2]]/sqrt(dr[2]^2+dr[3]^2) # chordwise strip normal vector
    chat = [ct, -st*shat[2], -st*shat[3]] # camberline vector
    ncp = cross(chat, dr) # normal vector perpindicular to camberline and bound vortex
    return ncp / norm(ncp) # normal vector used by AVL
end

@testset "AVL - Run 1 - Wing with Uniform Spacing" begin

    # Simple Wing with Uniform Spacing

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # adjust chord length so x-chord length matches AVL
    chord = @. chord/cos(theta)

    # vortex rings with symmetry
    mirror = false
    symmetric = true

    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.24324, atol=1e-3)
    @test isapprox(CD, 0.00243, atol=1e-5)
    @test isapprox(CDiff, 0.00245, atol=1e-5)
    @test isapprox(Cm, -0.02252, atol=1e-4)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)

    # vortex rings with mirrored geometry
    mirror = true
    symmetric = false

    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.24324, atol=1e-3)
    @test isapprox(CD, 0.00243, atol=1e-5)
    @test isapprox(CDiff, 0.00245, atol=1e-5)
    @test isapprox(Cm, -0.02252, atol=1e-4)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 2 - Wing with Cosine Spacing" begin

    # Run 2: Simple Wing with Cosine Spacing

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Cosine()
    spacing_c = Uniform()
    mirror = true
    symmetric = false

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # adjust chord length so x-chord length matches AVL
    chord = @. chord/cos(theta)

    # vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.23744, atol=1e-3)
    @test isapprox(CD, 0.00254, atol=1e-5)
    @test isapprox(CDiff, 0.00243, atol=1e-5)
    @test isapprox(Cm, -0.02165, atol=1e-4)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 3 - Wing at High Angle of Attack" begin

    # Simple Wing at High Angle of Attack

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = true

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 8.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # adjust chord length so x-chord length matches AVL
    chord = @. chord/cos(theta)

    # vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.80348, atol=1e-3)
    @test isapprox(CD, 0.02651, atol=1e-4)
    @test isapprox(CDiff, 0.02696, atol=1e-5)
    @test isapprox(Cm, -0.07399, atol=1e-3)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 4 - Wing with Dihedral" begin

    # Simple Wing with Dihedral

    # NOTE: There is some interaction between twist, dihedral, and chordwise
    # position which causes the normal vectors found by AVL to differ from those
    # computed by this package.  We therefore manually overwrite the normal
    # vectors when this occurs in order to get a better comparison.

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 3.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = true

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # adjust chord length so x-chord length matches AVL
    chord = @. chord/cos(theta)

    # also get normal vector as AVL defines it
    ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

    # vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    for (ip, p) in enumerate(surface)
        # check that our normal vector is approximately the same as AVL's
        @test isapprox(p.ncp, ncp, rtol=0.01)
        # replace our normal vector with AVL's normal vector for this test
        surface[ip] = set_normal(p, ncp)
    end

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.24787, atol=1e-3)
    @test isapprox(CD, 0.00246, atol=1e-4)
    @test isapprox(CDiff, 0.00245, atol=1e-5)
    @test isapprox(Cm, -0.02395, atol=1e-4)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 5 - Wing with Dihedral at Very High Angle of Attack" begin

    # Simple Wing with Dihedral at Very High Angle of Attack

    # NOTE: this test case is nonphysical, so it just tests the numerics

    # NOTE: There is some interaction between twist, dihedral, and chordwise
    # position which causes the normal vectors found by AVL to differ from those
    # computed by this package.  We therefore manually overwrite the normal
    # vectors when this occurs in order to get a better comparison.

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 3.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = true

    # adjust chord length to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 20.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

    # vortex rings, untwisted geometry
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    for (ip, p) in enumerate(surface)
        # check that our normal vector is approximately the same as AVL's
        @test isapprox(p.ncp, ncp, rtol=0.01)
        # replace our normal vector with AVL's normal vector for this test
        surface[ip] = set_normal(p, ncp)
    end

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 1.70982, rtol=0.02)
    @test isapprox(CD, 0.12904, rtol=0.02)
    @test isapprox(CDiff, 0.11502, rtol=0.02)
    @test isapprox(Cm, -0.45606, rtol=0.02)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 6 - Wing and Tail without Finite Core Model" begin

    # Wing and Tail without Finite Core Model

    # NOTE: AVL's finite-core model is turned off for these tests

    # NOTE: There is some interaction between twist, dihedral, and chordwise
    # position which causes the normal vectors found by AVL to differ from those
    # computed by this package.  We therefore manually overwrite the normal
    # vectors when this occurs in order to get a better comparison.

    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 6
    nc_h = 1
    spacing_s_h = Uniform()
    spacing_c_h = Uniform()
    mirror_h = false

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 5
    nc_v = 1
    spacing_s_v = Uniform()
    spacing_c_v = Uniform()
    mirror_v = false

    # adjust chord lengths to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)
    chord_h = @. chord_h/cos(theta_h)
    chord_v = @. chord_v/cos(theta_v)

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = [true, true, false]

    ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

    # vortex rings - finite core deactivated
    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    for (ip, p) in enumerate(wing)
        # check that our normal vector is approximately the same as AVL's
        @test isapprox(p.ncp, ncp, rtol=0.01)
        # replace our normal vector with AVL's normal vector for this test
        wing[ip] = set_normal(p, ncp)
    end

    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    translate!(htail, [4.0, 0.0, 0.0])

    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    translate!(vtail, [4.0, 0.0, 0.0])

    surfaces = [wing, htail, vtail]
    surface_id = [1, 1, 1]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.60408, atol=1e-3)
    @test isapprox(CD, 0.01058, atol=1e-4)
    @test isapprox(CDiff, 0.010378, atol=1e-4)
    @test isapprox(Cm, -0.02778, atol=2e-3)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 7 - Wing and Tail with Finite Core Model" begin

    # Wing and Tail with Finite Core Model

    # NOTE: There is some interaction between twist, dihedral, and chordwise
    # position which causes the normal vectors found by AVL to differ from those
    # computed by this package.  We therefore manually overwrite the normal
    # vectors when this occurs in order to get a better comparison.

    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 6
    nc_h = 1
    spacing_s_h = Uniform()
    spacing_c_h = Uniform()
    mirror_h = false

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 5
    nc_v = 1
    spacing_s_v = Uniform()
    spacing_c_v = Uniform()
    mirror_v = false

    # adjust chord lengths to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)
    chord_h = @. chord_h/cos(theta_h)
    chord_v = @. chord_v/cos(theta_v)

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = [true, true, false]

    ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c, fcore = (c,Δs)->max(c/4, Δs/2))

    for (ip, p) in enumerate(wing)
        # check that our normal vector is approximately the same as AVL's
        @test isapprox(p.ncp, ncp, rtol=0.01)
        # replace our normal vector with AVL's normal vector for this test
        wing[ip] = set_normal(p, ncp)
    end

    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    translate!(htail, [4.0, 0.0, 0.0])

    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    translate!(vtail, [4.0, 0.0, 0.0])

    surfaces = [wing, htail, vtail]
    surface_id = [1, 2, 3]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, derivatives = false)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.60562, atol=1e-3)
    @test isapprox(CD, 0.01058, atol=1e-4)
    @test isapprox(CDiff, 0.0104855, atol=1e-4)
    @test isapprox(Cm, -0.03377, atol=2e-3)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "AVL - Run 8 - Wing with Chordwise Panels" begin

    # Simple Wing with Chordwise Panels

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 6
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false

    # adjust chord length to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = true

    # vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.24454, atol=1e-3)
    @test isapprox(CD, 0.00247, atol=1e-5)
    @test isapprox(CDiff, 0.00248, atol=1e-5)
    @test isapprox(Cm, -0.02091, atol=1e-4)
    @test isapprox(CY, 0.0, atol=1e-16)
    @test isapprox(Cl, 0.0, atol=1e-16)
    @test isapprox(Cn, 0.0, atol=1e-16)
end

@testset "AVL - Run 9 - Wing with Cosine-Spaced Spanwise and Chordwise Panels" begin

    # Simple Wing with Cosine-Spaced Spanwise and Chordwise Panels

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 6
    spacing_s = Cosine()
    spacing_c = Cosine()
    mirror = false

    # adjust chord length to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = true

    # vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.23879, atol=1e-3)
    @test isapprox(CD, 0.00249, atol=1e-5)
    @test isapprox(CDiff, 0.0024626, atol=1e-5)
    @test isapprox(Cm, -0.01995, atol=1e-4)
    @test isapprox(CY, 0.0, atol=1e-16)
    @test isapprox(Cl, 0.0, atol=1e-16)
    @test isapprox(Cn, 0.0, atol=1e-16)
end

@testset "AVL - Run 10 - Wing with Sideslip" begin

    # Simple Wing with Sideslip

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()

    # adjust chord length to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 15.0*pi/180
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

    # vortex rings
    mirror = true
    symmetric = false

    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.22695, atol=1e-3)
    @test isapprox(CD, 0.00227, atol=1e-5)
    @test isapprox(CDiff, 0.0022852, atol=1e-5)
    @test isapprox(Cm, -0.02101, atol=1e-4)
    @test isapprox(CY, 0.0, atol=1e-5)
    @test isapprox(Cl, 0.00644, atol=1e-4)
    @test isapprox(Cn, -0.00012, atol=2e-4)
end

@testset "AVL - Run 11 - Wing Stability Derivatives" begin

    # Run 11: Simple Wing Stability Derivatives

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = true
    symmetric = false

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    dCF, dCM = stability_derivatives(system)

    CDa, CYa, CLa = dCF.alpha
    Cla, Cma, Cna = dCM.alpha
    CDb, CYb, CLb = dCF.beta
    Clb, Cmb, Cnb = dCM.beta
    CDp, CYp, CLp = dCF.p
    Clp, Cmp, Cnp = dCM.p
    CDq, CYq, CLq = dCF.q
    Clq, Cmq, Cnq = dCM.q
    CDr, CYr, CLr = dCF.r
    Clr, Cmr, Cnr = dCM.r

    @test isapprox(CLa, 4.638088, rtol=0.01)
    @test isapprox(CLb, 0.0, atol=ztol)
    @test isapprox(CYa, 0.0, atol=ztol)
    @test isapprox(CYb, -0.000007, atol=1e-4)
    @test isapprox(Cla, 0.0, atol=ztol)
    @test isapprox(Clb, 0.025749, rtol=0.015)
    @test isapprox(Cma, -0.429247, rtol=0.01)
    @test isapprox(Cmb, 0.0, atol=ztol)
    @test isapprox(Cna, 0.0, atol=ztol)
    # @test isapprox(Cnb, -0.000466, atol=1e-4)
    @test isapprox(Clp, -0.518725, rtol=0.01)
    @test isapprox(Clq, 0.0, atol=ztol)
    @test isapprox(Clr, 0.064243, rtol=0.01)
    @test isapprox(Cmp, 0.0, atol=ztol)
    @test isapprox(Cmq, -0.517094, rtol=0.01)
    @test isapprox(Cmr, 0.0, atol=ztol)
    @test isapprox(Cnp, -0.019846, rtol=0.01)
    @test isapprox(Cnq, 0.0, atol=ztol)
    @test isapprox(Cnr, -0.000898, rtol=0.01)
end

@testset "Lifting Line Coefficients" begin
    # Simple Wing with Uniform Spacing

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # adjust chord length so x-chord length matches AVL
    chord = @. chord/cos(theta)

    # vortex rings with symmetry
    mirror = false
    symmetric = true

    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    grids = [grid]
    surfaces = [surface]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)

    r_ll, c_ll = lifting_line_geometry(grids)

    cf, cm = lifting_line_coefficients(system, r_ll, c_ll; frame=Stability())

    cl_avl = [0.2618, 0.2646, 0.2661, 0.2664, 0.2654, 0.2628, 0.2584, 0.2513,
        0.2404, 0.2233, 0.1952, 0.1434]
    cd_avl = [0.0029, 0.0024, 0.0023, 0.0023, 0.0023, 0.0023, 0.0024, 0.0024,
        0.0025, 0.0026, 0.0026, 0.0022]
    cm_avl = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

    @test isapprox(cf[1][3,:], cl_avl, atol=1e-3, norm=(x)->norm(x, Inf))
    @test isapprox(cf[1][1,:], cd_avl, atol=1e-4, norm=(x)->norm(x, Inf))
    @test isapprox(cm[1][2,:], cm_avl, atol=1e-4, norm=(x)->norm(x, Inf))
end

@testset "Geometry Generation" begin

    # Tests of the geometry generation functions

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 6
    spacing_s = Uniform()
    spacing_c = Uniform()

    # Symmetric Geometry

    halfgrid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        spacing_s=spacing_s, spacing_c=spacing_c)

    halfgrid2, surface2 = grid_to_surface_panels(halfgrid1)

    halfgrid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc;
        spacing_s=spacing_s, spacing_c=spacing_c)

    surface4 = similar(surface1)
    VortexLattice.update_surface_panels!(surface4, halfgrid1)

    @test isapprox(halfgrid1, halfgrid2)

    for I in CartesianIndices(surface1)
        for field in fieldnames(SurfacePanel)
            @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field))
        end
    end

    @test isapprox(halfgrid1, halfgrid3)

    for I in CartesianIndices(surface1)
        for field in fieldnames(SurfacePanel)
            @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field))
        end
    end

    for I in CartesianIndices(surface1)
        for field in fieldnames(SurfacePanel)
            @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field))
        end
    end

    # Mirrored Geometry

    mirror = true

    grid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    grid2, surface2 = grid_to_surface_panels(halfgrid1; mirror = mirror)

    grid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; mirror=mirror,
        spacing_s=spacing_s, spacing_c=spacing_c)

    surface4 = similar(surface1)
    VortexLattice.update_surface_panels!(surface4, grid1)

    @test isapprox(grid1, grid2)

    for I in CartesianIndices(surface1)
        for field in fieldnames(SurfacePanel)
            @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field))
        end
    end

    @test isapprox(grid1, grid3)

    for I in CartesianIndices(surface1)
        for field in fieldnames(SurfacePanel)
            @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field))
        end
    end

    for I in CartesianIndices(surface1)
        for field in fieldnames(SurfacePanel)
            @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field))
        end
    end
end

@testset "Grid Input" begin

    # Tests for inputting a grid rather than a surface

    xle = [0.0, 0.4]
    yle = [0.0, 7.5]
    zle = [0.0, 0.0]
    chord = [2.2, 1.8]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 6
    spacing_s = Uniform()
    spacing_c = Uniform()

    Sref = 30.0
    cref = 2.0
    bref = 15.0
    rref = [0.50, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 1.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # adjust chord length so x-chord length matches AVL
    chord = @. chord/cos(theta)

    mirror = false
    symmetric = true

    halfgrid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    # steady analysis, single grid

    grids = [halfgrid]

    system = steady_analysis(grids, ref, fs; symmetric=symmetric)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    @test isapprox(CL, 0.24454, atol=1e-3)
    @test isapprox(CD, 0.00247, atol=1e-5)
    @test isapprox(CDiff, 0.00248, atol=1e-5)
    @test isapprox(Cm, -0.02091, atol=1e-4)
    @test isapprox(CY, 0.0, atol=1e-16)
    @test isapprox(Cl, 0.0, atol=1e-16)
    @test isapprox(Cn, 0.0, atol=1e-16)

    # steady analysis multiple grids

    # NOTE: AVL's finite-core model is turned off for these tests

    # NOTE: There is some interaction between twist, dihedral, and chordwise
    # position which causes the normal vectors found by AVL to differ from those
    # computed by this package.  We therefore manually overwrite the normal
    # vectors when this occurs in order to get a better comparison.

    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 6
    nc_h = 1
    spacing_s_h = Uniform()
    spacing_c_h = Uniform()
    mirror_h = false

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 5
    nc_v = 1
    spacing_s_v = Uniform()
    spacing_c_v = Uniform()
    mirror_v = false

    # adjust chord lengths to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)
    chord_h = @. chord_h/cos(theta_h)
    chord_v = @. chord_v/cos(theta_v)

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = [true, true, false]

    ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180)

    # vortex rings - finite core deactivated
    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    translate!(hgrid, [4.0, 0.0, 0.0])
    translate!(htail, [4.0, 0.0, 0.0])

    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    translate!(vgrid, [4.0, 0.0, 0.0])
    translate!(vtail, [4.0, 0.0, 0.0])

    grids = [wgrid, hgrid, vgrid]
    surface_id = [1, 1, 1]

    system = steady_analysis(grids, ref, fs; symmetric=symmetric, surface_id=surface_id)

    CF, CM = body_forces(system; frame=Stability())

    CDiff = far_field_drag(system)

    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    # some differences are expected since we don't set the normal vector in
    # VortexLattice equal to that used by AVL

    @test isapprox(CL, 0.60408, rtol=0.01)
    @test isapprox(CD, 0.01058, rtol=0.01)
    @test isapprox(CDiff, 0.010378, atol=0.02)
    @test isapprox(Cm, -0.02778, rtol=0.01)
    @test isapprox(CY, 0.0, atol=ztol)
    @test isapprox(Cl, 0.0, atol=ztol)
    @test isapprox(Cn, 0.0, atol=ztol)
end

@testset "Update Trailing Edge Coefficients" begin

    # This test checks whether the function which updates the trailing edge
    # coefficients `update_trailing_edge_coefficients!` results in the same
    # trailing edge coefficients as `influnece_coefficients!`

    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 6
    nc_h = 1
    spacing_s_h = Uniform()
    spacing_c_h = Uniform()
    mirror_h = false

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 5
    nc_v = 1
    spacing_s_v = Uniform()
    spacing_c_v = Uniform()
    mirror_v = false

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = [true, true, false]

    # horseshoe vortices
    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    translate!(htail, [4.0, 0.0, 0.0])

    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    translate!(vtail, [4.0, 0.0, 0.0])

    surfaces = [wing, htail, vtail]
    surface_id = [1, 2, 2]

    # number of panels
    N = nc*ns + nc_h*ns_h + nc_v*ns_v
    AIC1 = zeros(N, N)

    VortexLattice.influence_coefficients!(AIC1, surfaces;
        symmetric = symmetric,
        trailing_vortices = fill(false, length(surfaces)),
        surface_id = surface_id)

    AIC2 = copy(AIC1)

    VortexLattice.update_trailing_edge_coefficients!(AIC2, surfaces;
        symmetric = symmetric,
        trailing_vortices = fill(false, length(surfaces)),
        surface_id = surface_id)

    @test isapprox(AIC1, AIC2)
end

@testset "Wake Induced Velocity" begin

    # This test constructs two surfaces, one using surface panels,
    # and one using wake panels.  It then tests that the wake panel
    # implementations of `induced_velocity` yields identical results to the
    # surface panel implementations

    # generate the surface/wake geometry
    nc = 5
    ns = 10

    x = range(0, 1, length = nc+1)
    y = range(0, 2, length = ns+1)
    z = range(0, 3, length = ns+1)

    surface = Matrix{SurfacePanel{Float64}}(undef, nc, ns)
    wake = Matrix{WakePanel{Float64}}(undef, nc, ns)
    Γ = rand(nc, ns)

    for i = 1:nc, j = 1:ns
        rtl = [x[i], y[j], z[j]]
        rtr = [x[i], y[j+1], z[j+1]]
        rbl = [x[i+1], y[j], z[j]]
        rbr = [x[i+1], y[j+1], z[j+1]]
        rcp = (rtl + rtr + rbl + rbr)/4
        ncp = cross(rcp - rtr, rcp - rtl)
        core_size = 0.1
        chord = 0.0 # only used for unsteady simuulations

        surface[i,j] = SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord)
        wake[i,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, Γ[i, j])
    end

    # Test induced velocity calculation at an arbitrary point in space:
    rcp = [10, 11, 12]

    # no finite core model, no symmetry, no trailing vortices

    Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:];
        finite_core = false,
        symmetric = false,
        trailing_vortices = false,
        xhat = [1, 0, 0])

    Vw = VortexLattice.induced_velocity(rcp, wake;
        finite_core = false,
        symmetric = false,
        trailing_vortices = false,
        xhat = [1, 0, 0])

    @test isapprox(Vs, Vw)

    # finite core model, symmetry, and trailing vortices

    Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:];
        finite_core = true,
        symmetric = true,
        trailing_vortices = true,
        xhat = [1, 0, 0])

    Vw = VortexLattice.induced_velocity(rcp, wake;
        finite_core = true,
        symmetric = true,
        trailing_vortices = true,
        xhat = [1, 0, 0])

    @test isapprox(Vs, Vw)

    # Test induced velocity calculation at a trailing edge point:
    is = 3 # trailing edge index
    I = CartesianIndex(nc+1, is) # index on surface

    # no finite core model, no symmetry, no trailing vortices
    Vs = VortexLattice.induced_velocity(is, surface, Γ[:];
        finite_core = false,
        symmetric = false,
        trailing_vortices = false,
        xhat = [1, 0, 0])

    Vw = VortexLattice.induced_velocity(I, wake;
        finite_core = false,
        symmetric = false,
        trailing_vortices = false,
        xhat = [1, 0, 0])

    @test isapprox(Vs, Vw)

    # finite core model, symmetry, and trailing vortices
    Vs = VortexLattice.induced_velocity(is, surface, Γ[:];
        finite_core = true,
        symmetric = true,
        trailing_vortices = true,
        xhat = [1, 0, 0])

    Vw = VortexLattice.induced_velocity(I, wake;
        finite_core = true,
        symmetric = true,
        trailing_vortices = true,
        xhat = [1, 0, 0])

    @test isapprox(Vs, Vw)

end

@testset "Unsteady Vortex Lattice Method - Rectangular Wing" begin

    Uinf = 1.0
    AR = 4

    # reference parameters
    cref = 1.0
    bref = AR
    Sref = bref*cref
    rref = [0.0, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    # freestream parameters
    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    # geometry
    xle = [0.0, 0.0]
    yle = [-bref/2, bref/2]
    zle = [0.0, 0.0]
    chord = [cref, cref]
    theta = [0.0, 0.0]
    phi = [0.0, 0.0]
    ns = 13
    nc = 4
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false
    symmetric = false

    # non-dimensional time
    t = range(0.0, 10.0, step=0.2)
    dt = [t[i+1]-t[i] for i = 1:length(t)-1]

    # create vortex rings
    grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    surfaces = [surface]

    # run analysis
    system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt;
        symmetric=symmetric)

    # extract forces at each time step
    CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind())
end

@testset "Unsteady Vortex Lattice Method - Wing + Tail" begin

    # Unsteady Wing and Tail

    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 12
    nc = 1
    spacing_s = Uniform()
    spacing_c = Uniform()
    mirror = false

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 6
    nc_h = 1
    spacing_s_h = Uniform()
    spacing_c_h = Uniform()
    mirror_h = false

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 5
    nc_v = 1
    spacing_s_v = Uniform()
    spacing_c_v = Uniform()
    mirror_v = false

    # adjust chord lengths to match AVL (which uses chord length in the x-direction)
    chord = @. chord/cos(theta)
    chord_h = @. chord_h/cos(theta_h)
    chord_v = @. chord_v/cos(theta_v)

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = Freestream(Vinf, alpha, beta, Omega)

    symmetric = [true, true, false]

    # horseshoe vortices
    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)

    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    translate!(htail, [4.0, 0.0, 0.0])

    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    translate!(vtail, [4.0, 0.0, 0.0])

    surfaces = [wing, htail, vtail]
    surface_id = [1, 2, 3]

    # t
    t = range(0.0, 10.0, step=0.2)
    dt = t[2:end] - t[1:end-1]

    system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; symmetric)

    # extract forces at each time step
    CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind())
end


text:
"""
    simpleflux(Ca::Float,Co::Float,pisvel::Float)

Compute flux entering the ocean per unit area

```
mld=10 #mixed layer depth (m)
tim=86400*30 #relaxation time scale (s)
pisvel=mld/tim #piston velocity (m/s)
Co=0.0 #ocean value (e.g. concentration of some compound)
Ca=1.0 #atmospeheric value (e.g. equivalent compound concentration)

simpleflux(Ca,Co,pisvel)
```
"""
function simpleflux(Ca,Co,pisvel)
    return flx=pisvel*(Ca-Co)
end


text:
using BinDeps
using Compat
using Conda

@BinDeps.setup

# Use x.y.z for downloading binaries, but only check for x.y in shared libraries
# so that patch version increases can be found in system
libsymengine_version = "0.3.0"
libsymengine_soversion = join(split(libsymengine_version, ".")[1:2], ".")

if is_windows()
   libsymengine_soname = "symengine-$(libsymengine_soversion).dll"
elseif is_apple()
   libsymengine_soname = "libsymengine.$(libsymengine_soversion).dylib"
else
   libsymengine_soname = "libsymengine.so.$(libsymengine_soversion)"
end

# Use the dummy name libsymengine-dummy to avoid finding libsymengine.so
# We only want to find the versioned library
libdep = library_dependency("libsymengine_dummy", aliases=[libsymengine_soname])

path = abspath(dirname(@__FILE__), "usr")
if is_windows()
    isdir(path) || mkdir(path)
    if Sys.WORD_SIZE == 64
        suffix = "x86_64"
    else
        suffix = "x86"
    end
    url = "https://github.com/symengine/symengine/releases/download/"
    url *= "v$(libsymengine_version)/symengine-$(libsymengine_version)-binaries-msvc-$(suffix).tar.bz2"
    provides(Binaries, URI(url), libdep, unpacked_dir="symengine-$(libsymengine_version)/bin")
else
    env = Symbol(path)
    EnvManagerType = Conda.EnvManager{env}
    # Conda's method will install miniconda to check that a package exists.
    # This will indicate to BinDeps that the packages for this env exists unconditionally.
    BinDeps.package_available(m::EnvManagerType) = true
    # Adding the channel to conda will install miniconda even if conda is not the chosen LibraryProvider
    # Override the command for install in this env so that channels are added.
    function BinDeps.generate_steps(dep::BinDeps.LibraryDependency, manager::EnvManagerType, opts)
        Conda.add_channel("conda-forge", env)
        Conda.add("$(manager.packages[1])", env)
    end
    provides(EnvManagerType, "symengine=$(libsymengine_version)", [libdep])
end

@BinDeps.install Dict([(:libsymengine_dummy, :libsymengine)])


text:
using LuxurySparse
using TupleTools
@static if hasmethod(TupleTools.diff, Tuple{Tuple{}})
    tuple_diff(args...) = TupleTools.diff(args...)
else
    tuple_diff(v::Tuple{}) = () # similar to diff([])
    tuple_diff(v::Tuple{Any}) = ()
    tuple_diff(v::Tuple) = (v[2] - v[1], tuple_diff(Base.tail(v))...)
end

"""
    sort_unitary(U, locations::NTuple{N, Int}) -> U

Return an sorted unitary operator according to the locations.
"""
function sort_unitary(U::AbstractMatrix, locs::NTuple{N,Int}) where {N}
    if all(each > 0 for each in tuple_diff(locs))
        return U
    else
        return reorder(U, TupleTools.sortperm(locs))
    end
end

"""
    swaprows!(v::VecOrMat, i::Int, j::Int[, f1, f2]) -> VecOrMat

swap row i and row j of v inplace, with f1, f2 factors applied on i and j (before swap).
"""
function swaprows! end

Base.@propagate_inbounds function swaprows!(
    v::AbstractMatrix{T},
    i::Int,
    j::Int,
    f1,
    f2,
) where {T}
    for c = 1:size(v, 2)
        temp = v[i, c]
        v[i, c] = v[j, c] * f2
        v[j, c] = temp * f1
    end
    v
end

Base.@propagate_inbounds function swaprows!(v::AbstractMatrix{T}, i::Int, j::Int) where {T}
    for c = 1:size(v, 2)
        temp = v[i, c]
        v[i, c] = v[j, c]
        v[j, c] = temp
    end
    v
end

"""
    swapcols!(v::VecOrMat, i::Int, j::Int[, f1, f2]) -> VecOrMat

swap col i and col j of v inplace, with f1, f2 factors applied on i and j (before swap).
"""
function swapcols! end

Base.@propagate_inbounds function swapcols!(
    v::AbstractMatrix{T},
    i::Int,
    j::Int,
    f1,
    f2,
) where {T}
    for c = 1:size(v, 1)
        temp = v[c, i]
        v[c, i] = v[c, j] * f2
        v[c, j] = temp * f1
    end
    v
end

Base.@propagate_inbounds function swapcols!(v::AbstractMatrix{T}, i::Int, j::Int) where {T}
    for c = 1:size(v, 1)
        temp = v[c, i]
        v[c, i] = v[c, j]
        v[c, j] = temp
    end
    v
end

Base.@propagate_inbounds swapcols!(v::AbstractVector, args...) = swaprows!(v, args...)

Base.@propagate_inbounds function swaprows!(v::AbstractVector, i::Int, j::Int, f1, f2)
    temp = v[i]
    v[i] = v[j] * f2
    v[j] = temp * f1
    v
end

Base.@propagate_inbounds function swaprows!(v::AbstractVector, i::Int, j::Int)
    temp = v[i]
    v[i] = v[j]
    v[j] = temp
    v
end

"""
    u1rows!(state::VecOrMat, i::Int, j::Int, a, b, c, d) -> VecOrMat

apply u1 on row i and row j of state inplace.
"""
function u1rows! end

Base.@propagate_inbounds function u1rows!(state::AbstractVector, i::Int, j::Int, a, b, c, d)
    w = state[i]
    v = state[j]
    state[i] = a * w + b * v
    state[j] = c * w + d * v
    state
end

Base.@propagate_inbounds function u1rows!(state::AbstractMatrix, i::Int, j::Int, a, b, c, d)
    for col = 1:size(state, 2)
        w = state[i, col]
        v = state[j, col]
        state[i, col] = a * w + b * v
        state[j, col] = c * w + d * v
    end
    state
end

"""
    mulrow!(v::AbstractVector, i::Int, f) -> VecOrMat

multiply row i of v by f inplace.
"""
function mulrow! end

Base.@propagate_inbounds function mulrow!(v::AbstractVector, i::Int, f)
    v[i] *= f
    return v
end

Base.@propagate_inbounds function mulrow!(v::AbstractMatrix, i::Int, f)
    for j = 1:size(v, 2)
        v[i, j] *= f
    end
    return v
end

"""
    mulcol!(v::AbstractVector, i::Int, f) -> VecOrMat

multiply col i of v by f inplace.
"""
function mulcol! end

Base.@propagate_inbounds function mulcol!(v::AbstractVector, i::Int, f)
    v[i] *= f
    return v
end

Base.@propagate_inbounds function mulcol!(v::AbstractMatrix, j::Int, f)
    for i = 1:size(v, 1)
        v[i, j] *= f
    end
    return v
end

"""
    matvec(x::VecOrMat) -> MatOrVec

Return vector if a matrix is a column vector, else untouched.
"""
function matvec end

matvec(x::AbstractMatrix) = size(x, 2) == 1 ? vec(x) : x
matvec(x::AbstractVector) = x

Base.@propagate_inbounds function unrows!(
    state::AbstractVector,
    inds::AbstractVector,
    U::AbstractMatrix,
)
    state[inds] = U * state[inds]
    return state
end

Base.@propagate_inbounds function unrows!(
    state::AbstractMatrix,
    inds::AbstractVector,
    U::AbstractMatrix,
)
    for k = 1:size(state, 2)
        state[inds, k] = U * state[inds, k]
    end
    return state
end

############# boost unrows! for sparse matrices ################
@inline unrows!(state::AbstractVector, inds::AbstractVector, U::IMatrix) = state

Base.@propagate_inbounds function unrows!(
    state::AbstractVector,
    inds::AbstractVector,
    U::SDDiagonal,
)
    for i = 1:length(U.diag)
        state[inds[i]] *= U.diag[i]
    end
    return state
end

Base.@propagate_inbounds function unrows!(
    state::AbstractMatrix,
    inds::AbstractVector,
    U::SDDiagonal,
)
    for j = 1:size(state, 2)
        for i = 1:length(U.diag)
            state[inds[i], j] *= U.diag[i]
        end
    end
    state
end

Base.@propagate_inbounds function unrows!(
    state::AbstractVector,
    inds::AbstractVector,
    U::SDPermMatrix,
)
    state[inds] = state[inds[U.perm]] .* U.vals
    state
end

Base.@propagate_inbounds function unrows!(
    state::AbstractMatrix,
    inds::AbstractVector,
    U::SDPermMatrix,
)
    for k = 1:size(state, 2)
        state[inds, k] = state[inds[U.perm], k] .* U.vals
    end
    state
end

Base.@propagate_inbounds function unrows!(
    state::AbstractVector,
    inds::AbstractVector,
    A::SDSparseMatrixCSC,
    work::AbstractVector,
)
    work .= 0
    for col = 1:length(inds)
        xj = state[inds[col]]
        for j = A.colptr[col]:(A.colptr[col+1]-1)
            work[A.rowval[j]] += A.nzval[j] * xj
        end
    end
    state[inds] = work
    state
end

Base.@propagate_inbounds function unrows!(
    state::AbstractMatrix,
    inds::AbstractVector,
    A::SDSparseMatrixCSC,
    work::Matrix,
)
    work .= 0
    for k = 1:size(state, 2)
        for col = 1:length(inds)
            xj = state[inds[col], k]
            for j = A.colptr[col]:(A.colptr[col+1]-1)
                work[A.rowval[j], k] += A.nzval[j] * xj
            end
        end
        state[inds, k] = view(work, :, k)
    end
    state
end

using LinearAlgebra: Transpose
Base.convert(::Type{Transpose{T,Matrix{T}}}, arr::AbstractMatrix{T}) where {T} =
    transpose(Matrix(transpose(arr)))
Base.convert(t::Type{Transpose{T,Matrix{T}}}, arr::Transpose{T}) where {T} =
    invoke(convert, Tuple{Type{Transpose{T,Matrix{T}}},Transpose}, t, arr)

function ispow(x::Integer, d::Integer)
    r = log(x) / log(d)
    return round(Int, r) ≈ r
end


text:
function gen_drmesh_exp(r_min, r_max, a, Nr)

    N = Nr - 1

    drmesh = zeros(Float64,N+1)
    SMALL = eps()

    if a < 0
        error("gen_drmesh_exp: a > 0 required")
    #
    elseif abs(a - 1) < SMALL
        error("gen_drmesh_exp: a == 1 not implemented")
    #
    else
        if N > 1
            β = log(a)/(N-1)
            α = (r_max - r_min) / (exp(β*N) - 1)
            for i in 1:N+1
                drmesh[i] = α * β * exp( β*(i-1) )
            end
        else
            error("mesh_exp_deriv: N > 1 required")
        end
    end
    return drmesh
end


text:
using AugmentedGPLikelihoods
using AugmentedGPLikelihoods.SpecialDistributions
using AugmentedGPLikelihoods.TestUtils: test_auglik
const AGPL = AugmentedGPLikelihoods
using Distributions
using GPLikelihoods
using GPLikelihoods.TestInterface: test_interface
using LogExpFunctions
using MeasureBase
using Random
using SpecialFunctions
using Test
using TupleVectors

@testset "AugmentedGPLikelihoods.jl" begin
    @info "Testing likelihoods"
    @testset "Likelihoods" begin
        include("likelihoods/bernoulli.jl")
        include("likelihoods/categorical.jl")
        include("likelihoods/laplace.jl")
        include("likelihoods/negativebinomial.jl")
        include("likelihoods/poisson.jl")
        include("likelihoods/studentt.jl")
    end

    @info "Testing SpecialDistributions"
    @testset "SpecialDistributions" begin
        include("SpecialDistributions/negativemultinomial.jl")
        include("SpecialDistributions/ntdist.jl")
        include("SpecialDistributions/polyagamma.jl")
        include("SpecialDistributions/polyagammanegativemultinomial.jl")
        include("SpecialDistributions/polyagammapoisson.jl")
    end

    include("utils.jl")
end


text:
# This test suite exploits the Weak Law of Large Numbers
# to verify that all of the functions defined on a
# distribution produce values within a close range of
# their theoretically predicted values.
#
# This includes tests of means, variances, skewness and kurtosis;
# as well as tests of more complex quantities like quantiles,
# entropy, etc...
#
# These tests are quite slow, but are essential to verifying the
# accuracy of our distributions

using Distributions
using Base.Test

probpts(n::Int) = ((1:n) - 0.5)/n
const pp  = float(probpts(100))
const lpp = log(pp)

# Use a large, odd number of samples for testing all quantities
n_samples = 5_000_001

# Try out many parameterizations of any given distribution
for d in [Arcsine(),
          Beta(2.0, 2.0),
          Beta(3.0, 4.0),
          Beta(17.0, 13.0),
          # BetaPrime(3.0, 3.0),
          # BetaPrime(3.0, 5.0),
          # BetaPrime(5.0, 3.0),
          Cauchy(0.0, 1.0),
          Cauchy(10.0, 1.0),
          Cauchy(0.0, 10.0),
          Chi(12),
          Chisq(8),
          Chisq(12.0),
          Chisq(20.0),
          # Cosine(),
          # Empirical(),
          Erlang(1),
          Erlang(17.0),
          Exponential(1.0),
          Exponential(5.1),
          FDist(9, 9),
          FDist(9, 21),
          FDist(21, 9),
          Gamma(3.0, 2.0),
          Gamma(2.0, 3.0),
          Gamma(3.0, 3.0),
          Geometric(0.1),
          Geometric(0.5),
          Geometric(0.9),
          Gumbel(3.0, 5.0),
          Gumbel(5, 3),
          # HyperGeometric(1.0, 1.0, 1.0),
          # HyperGeometric(2.0, 2.0, 2.0),
          # HyperGeometric(3.0, 2.0, 2.0),
          # HyperGeometric(2.0, 3.0, 2.0),
          # HyperGeometric(2.0, 2.0, 3.0),
          InverseGaussian(1.0,1.0),
          InverseGaussian(2.0,7.0),
          InvertedGamma(1.0,1.0),
          InvertedGamma(2.0,3.0),
          Laplace(0.0, 1.0),
          Laplace(10.0, 1.0),
          Laplace(0.0, 10.0),
          Levy(0.0, 1.0),
          Levy(2.0, 8.0),
          Levy(3.0, 3.0),
          Logistic(0.0, 1.0),
          Logistic(10.0, 1.0),
          Logistic(0.0, 10.0),
          LogNormal(0.0, 1.0),
          LogNormal(10.0, 1.0),
          LogNormal(0.0, 10.0),
          NegativeBinomial(),
          NegativeBinomial(5, 0.6),
          NoncentralBeta(2,2,0),
          NoncentralBeta(2,6,5),
          NoncentralChisq(2,2),
          NoncentralChisq(2,5),
          NoncentralF(2,2,2),
          NoncentralF(8,10,5),
          NoncentralT(2,2),
          NoncentralT(10,2),
          Normal(0.0, 1.0),
          Normal(-1.0, 10.0),
          Normal(1.0, 10.0),
          # Pareto(),
          Poisson(2.0),
          Poisson(10.0),
          Poisson(51.0),
          Rayleigh(1.0),
          Rayleigh(5.0),
          Rayleigh(10.0),
          # Skellam(10.0, 2.0), # Entropy wrong
          # TDist(1), # Entropy wrong
          # TDist(28), # Entropy wrong
          Triangular(3.0, 1.0),
          Triangular(3.0, 2.0),
          Triangular(10.0, 10.0),
          Truncated(Normal(0, 1), -3, 3),
          # Truncated(Normal(-100, 1), 0, 1),
          Truncated(Normal(27, 3), 0, Inf),
          Uniform(0.0, 1.0),
          Uniform(3.0, 17.0),
          Uniform(3.0, 3.1),
          Weibull(2.3),
          Weibull(23.0),
          Weibull(230.0)]

    # NB: Uncomment if test fails
    # Mention distribution being run
    # println(d)

    n = length(pp)
    is_continuous = isa(d, Truncated) ? isa(d.untruncated, ContinuousDistribution) : isa(d, ContinuousDistribution)
    is_discrete = isa(d, Truncated) ? isa(d.untruncated, DiscreteDistribution) : isa(d, DiscreteDistribution)

    @assert is_continuous == !is_discrete
    sample_ty = is_continuous ? Float64 : Int

    # avoid checking high order moments for LogNormal and Logistic
    avoid_highord = isa(d, LogNormal) || isa(d, Logistic) || isa(d, Truncated)

    #####
    #
    #  Part 1: Capability of random number generation
    #
    #####

    # check that we can generate a single random draw
    draw = rand(d)

    # check that draw satifies insupport()
    @test insupport(d, draw)

    # check that we can generate many random draws at once
    x = rand(d, n)
    @test nsamples(typeof(d), x) == n

    # check that sequence of draws satifies insupport()
    @test insupport(d, x)

    # check that we can generate many random draws in-place
    rand!(d, x)

    #####
    #
    #  Part 2: Evaluation
    #  ----------------------
    #
    #  This part tests the integrity/consistency of following functions:
    #
    #  - pdf
    #  - cdf
    #  - ccdf
    #
    #  - logpdf
    #  - logcdf
    #  - logccdf
    #
    #  - quantile
    #  - cquantile
    #  - invlogcdf
    #  - invlogccdf
    #
    #####

    # evaluate by scalar

    x = zeros(sample_ty, n)
    r_cquan = zeros(sample_ty, n)

    r_pdf = zeros(n)
    r_cdf = zeros(n)
    r_ccdf = zeros(n)

    r_logpdf = zeros(n)
    r_logcdf = zeros(n)
    r_logccdf = zeros(n)

    r_invlogcdf = zeros(n)
    r_invlogccdf = zeros(n)

    for i in 1:n
        x[i] = quantile(d, pp[i])
        r_cquan[i] = cquantile(d, pp[i])

        xi = x[i]
        r_pdf[i] = pdf(d, xi)
        r_cdf[i] = cdf(d, xi)
        r_ccdf[i] = ccdf(d, xi)

        r_logpdf[i] = logpdf(d, xi)
        r_logcdf[i] = logcdf(d, xi)
        r_logccdf[i] = logccdf(d, xi)

        r_invlogcdf[i] = invlogcdf(d, lpp[i])
        r_invlogccdf[i] = invlogccdf(d, lpp[i])
    end

    # testing consistency between scalar evaluation and vectorized evaluation

    for i in 1:length(x)
        @test_approx_eq quantile(d, pp)  x
        @test_approx_eq cquantile(d, pp) r_cquan
        @test_approx_eq pdf(d, x)        r_pdf
        @test_approx_eq cdf(d, x)        r_cdf
        @test_approx_eq ccdf(d, x)       r_ccdf
        @test_approx_eq logpdf(d, x)     r_logpdf
        @test_approx_eq logcdf(d, x)     r_logcdf
        @test_approx_eq logccdf(d, x)    r_logccdf

        @test_approx_eq invlogcdf(d, lpp)  r_invlogcdf
        @test_approx_eq invlogccdf(d, lpp) r_invlogccdf
    end

    # # testing consistency between different functions

    @test_approx_eq logpdf(d, x) log(pdf(d, x))
    @test_approx_eq cquantile(d, 1 - pp) x

    if is_continuous
        @test_approx_eq cdf(d, x) pp
        @test_approx_eq ccdf(d, x) 1 - pp
        @test_approx_eq logcdf(d, x) lpp
        @test_approx_eq logccdf(d, x) lpp[end:-1:1]
        @test_approx_eq invlogcdf(d, lpp) x
        @test_approx_eq invlogccdf(d, lpp) x[end:-1:1]
    end

    # TODO: Test mgf, cf

    #####
    #
    #  Part 3: Other tests
    #
    #####

    # Test modes by looking at pdf(x +/- eps()) near a mode x

    if !isa(d, Uniform) && method_exists(modes,(typeof(d),))
        ms = modes(d)
        if isa(d, ContinuousUnivariateDistribution)
            if insupport(d, ms[1] + 0.1)
                @test pdf(d, ms[1]) > pdf(d, ms[1] + 0.1)
            end
            if insupport(d, ms[1] - 0.1)
                @test pdf(d, ms[1]) > pdf(d, ms[1] - 0.1)
            end
        end
    end
end


text:
import SpecialFunctions: gamma

abstract type AbstractWindResourceModel end

"""
    DiscritizedWindResource(wind_directions, wind_speeds, wind_probabilities, measurement_heights, air_density, ti_model, wind_shear_model)

    Struct defining a wind resource

# Arguments
- `wind_directions::Array{Float,1}(Nstates)`: an array of wind directions corresponding to each wind farm state in radians
- `wind_speeds::Array{Float,1}(Nstates)`: an array of wind speeds corresponding to each wind farm state in meters/second
- `wind_probabilities::Array{Float,1}(Nstates)`: an array of probabilities corresponding to each wind farm state with values in [0,1]
- `measurement_heights::Array{Float,1}(Nstates)`: an array of measurement heights corresponding to each wind farm state
- `air_density::Float`: the air density
- `ambient_ti::Array{Float,1}`: an array of the ambient turbulence intensity for each wind direction
- `wind_shear_model::Array{AbstractWindShearModel}(1)`: contains a struct defining the desired turbulence intensity model
"""
struct DiscretizedWindResource{AF, TF, ASM} <: AbstractWindResourceModel

    wind_directions::AF
    wind_speeds::AF
    wind_probabilities::AF
    measurement_heights::AF
    air_density::TF
    ambient_tis::AF
    wind_shear_model::ASM

end

"""
    rediscretize_windrose(windrosein::DiscretizedWindResource, ndirectionbins, nspeedbins)

    Function for re-interpreting a wind rose into a desired number of bins. Returns the new
    wind rose. Currently only works for windroses with a single speed in each direction.

# Arguments
- `windrosein::DiscretizedWindResource`: original wind rose
- `ndirectionbins::Integer`: number of direction bins for the new wind rose
- `start::Float`: direction for first bin in radians
- `averagespeed::Bool`: set whether or not to return the average wind speed as the speed for all bins
"""
function rediscretize_windrose(windrosein::DiscretizedWindResource, ndirectionbins; start=0.0, averagespeed=false)

    # create Akima spline interpolation of windrosein
    splinedirs = [-2*pi.+windrosein.wind_directions; windrosein.wind_directions; 2*pi.+windrosein.wind_directions]
    speedspline = Akima(splinedirs, [windrosein.wind_speeds; windrosein.wind_speeds; windrosein.wind_speeds])
    probabilityspline = Akima(splinedirs, [windrosein.wind_probabilities; windrosein.wind_probabilities; windrosein.wind_probabilities])
    heightspline = Akima(splinedirs, [windrosein.measurement_heights; windrosein.measurement_heights; windrosein.measurement_heights])
    ambienttispline = Akima(splinedirs, [windrosein.ambient_tis; windrosein.ambient_tis; windrosein.ambient_tis])

    # get new interpolated wind rose attributes
    directionsnew = collect(range(start,2*pi+start-2*pi/ndirectionbins,length=ndirectionbins))
    if averagespeed
        speedsnew = speedspline.(ones(ndirectionbins)*sum(directionsnew)/ndirectionbins)
    else
        speedsnew = speedspline.(directionsnew)
    end
    probabilitiesnew = probabilityspline.(directionsnew)
    heightsnew = heightspline.(directionsnew)
    ambient_tis_new = ambienttispline.(directionsnew)

    # re-normalize probabilites
    probabilitiesnew = probabilitiesnew./sum(probabilitiesnew)

    # create new wind rose
    windroseout = DiscretizedWindResource(directionsnew, speedsnew, probabilitiesnew, heightsnew, windrosein.air_density, ambient_tis_new, windrosein.wind_shear_model)

    return windroseout
end

function _dist_weibull(x,L;k=2.0)
    if L < 0.0001
        L = 0.0001
    end
    return (k/L)*(x/L)^(k-1)*exp(-(x/L)^k)
end

function setup_weibull_distribution(wind_directions,wind_probabilities,wind_speeds,nspeeds)
        if nspeeds == 1
            return wind_directions, wind_probabilities, wind_speeds
        end

        ndirections = length(wind_directions)
        dirs = zeros(ndirections*nspeeds)
        freqs = zeros(ndirections*nspeeds)
        speeds = zeros(ndirections*nspeeds)

        #direction loops
        k = 1
        for i = 1:ndirections
            for j = 1:nspeeds
                dirs[k] = wind_directions[i]
                k += 1
            end
        end

        #speed and frequency loops
        m = 1
        for i = 1:ndirections
            avg_speed = wind_speeds[i]
            speed_dist = range(25.0/(2.0*nspeeds)+0.001,stop=25.0-25.0/(2.0*nspeeds),length=nspeeds)
            dspeed = speed_dist[2]-speed_dist[1]
            num_int = 1000

            for j = 1:nspeeds
                speed_int = range(speed_dist[j]-dspeed/2.0,stop=speed_dist[j]+dspeed/2.0,length=num_int)
                k = 2.0
                scale = avg_speed/(gamma(1.0+1.0/k))
                freq_int = _dist_weibull.(speed_int,scale)

                speed_freq = trapz(speed_int,freq_int)
                speeds[m] = speed_dist[j]
                freqs[m] = speed_freq*wind_probabilities[i]
                m += 1
            end
        end

        return dirs, freqs, speeds
end