text:
A = zeros(10,10)
A[3,4] = 1
A[5,6] = 1
A[10,1] = 10
p, err = track(A, (5,5), 3)
@test norm(p .- [4,5]) < sqrt(eps())


text:
using ITensorNetworkAD
using Test

@testset "ITensorNetworks.jl" begin
  for filename in [
    "lattice.jl",
    "peps.jl",
    "models.jl",
    "mpstensor.jl",
    "treetensor.jl",
    "itensor_network.jl",
    "projectors.jl",
    "interface.jl",
  ]
    println("Running $filename in ITensorNetworks.jl")
    include(filename)
  end
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PentasphericalVariogram(range=r, sill=s, nugget=n)
    PentasphericalVariogram(ball; sill=s, nugget=n)

A pentaspherical variogram with range `r`, sill `s` and nugget `n`.
Optionally, use a custom metric `ball`.
"""
struct PentasphericalVariogram{T,B} <: Variogram
  sill::T
  nugget::T
  ball::B
end

PentasphericalVariogram(ball; sill=1.0, nugget=0.0) =
  PentasphericalVariogram(sill, nugget, ball)

PentasphericalVariogram(; range=1.0, sill=1.0, nugget=0.0) =
  PentasphericalVariogram(sill, nugget, MetricBall(range))

function (γ::PentasphericalVariogram)(h::T) where {T}
  r = radius(γ.ball)
  s = γ.sill
  n = γ.nugget

  # constants
  c1 = T(15) / T(8)
  c2 = T(5)  / T(4)
  c3 = T(3)  / T(8)
  s1 = c1*(h/r) - c2*(h/r)^3 + c3*(h/r)^5
  s2 = T(1)

  (h < r) * (s - n) * s1 +
  (h ≥ r) * (s - n) * s2 +
  (h > 0) * n
end

isstationary(::Type{<:PentasphericalVariogram}) = true


text:
"Create infeasible state trajectory initialization problem from problem"
function infeasible_problem(prob::Problem{T},R_inf::T=1.0) where T
    N = prob.N; n = prob.model.n; m = prob.model.m
    @assert all([prob.obj[k] isa QuadraticCost for k = 1:N]) #TODO generic cost

    # modify problem with slack control
    obj_inf = CostFunction[]
    for k = 1:N-1
        cost_inf = copy(prob.obj[k])
        cost_inf.R = cat(cost_inf.R,R_inf*Diagonal(I,prob.model.n)/prob.dt,dims=(1,2))
        cost_inf.r = [cost_inf.r; zeros(prob.model.n)]
        cost_inf.H = [cost_inf.H; zeros(prob.model.n,prob.model.n)]
        push!(obj_inf,cost_inf)
    end
    push!(obj_inf,copy(prob.obj[N]))

    model_inf = add_slack_controls(prob.model)
    u_slack = slack_controls(prob)
    con_inf = infeasible_constraints(prob.model.n,prob.model.m)

    con_prob = ConstraintSet[]
    constrained = is_constrained(prob)
    for k = 1:N-1
        _con = GeneralConstraint[]
        constrained ? append!(_con,update_constraint_set_jacobians(prob.constraints.C[k],n,n,m)) : nothing
        push!(_con,con_inf)
        push!(con_prob,_con)
    end

    constrained ? push!(con_prob,prob.constraints.C[N]) : push!(con_prob,GeneralConstraint[])

    update_problem(prob,model=model_inf,obj=Objective(obj_inf),
        constraints=Constraints(con_prob),U=[[prob.U[k];u_slack[k]] for k = 1:prob.N-1])
end

"Return a feasible problem from an infeasible problem"
function infeasible_to_feasible_problem(prob::Problem{T},prob_altro::Problem{T},
        state::NamedTuple,opts::ALTROSolverOptions{T}) where T
    prob_altro_feasible = copy(prob)

    if state.minimum_time
        prob_altro_feasible = minimum_time_problem(prob_altro_feasible,opts.R_minimum_time,
            opts.dt_max,opts.dt_min)

        # initialize sqrt(dt) from previous solve
        for k = 1:prob.N-1
            prob_altro_feasible.U[k][end] = prob_altro.U[k][end]
            k != 1 ? prob_altro_feasible.X[k][end] = prob_altro.X[k][end] : prob_altro_feasible.X[k][end] = 0.0
        end
        prob_altro_feasible.X[end][end] = prob_altro.X[end][end]
    end

    if opts.dynamically_feasible_projection
        projection!(prob_altro_feasible, opts.opts_al.opts_uncon)
    end

    return prob_altro_feasible
end

"Calculate slack controls that produce infeasible state trajectory"
function slack_controls(prob::Problem{T}) where T
    N = prob.N
    n = prob.model.n
    m = prob.model.m

    dt = prob.dt

    u_slack = [zeros(n) for k = 1:N-1]#zeros(n,N-1)
    x = [zeros(n) for k = 1:N]#zeros(n,N)
    x[1] = prob.x0

    for k = 1:N-1
        evaluate!(x[k+1],prob.model,x[k],prob.U[k],dt)

        u_slack[k] = prob.X[k+1] - x[k+1]
        x[k+1] += u_slack[k]
    end
    return u_slack
end

function line_trajectory(x0::Vector,xf::Vector,N::Int)
    t = range(0,stop=N,length=N)
    slope = (xf .- x0)./N
    x_traj = [slope*t[k] for k = 1:N]
    x_traj[1] = x0
    x_traj[end] = xf
    x_traj
end

function line_trajectory!(X::VectorTrajectory,x0::Vector,xf::Vector)
    N = length(X)
    t = range(0,stop=N,length=N)
    slope = (xf .- x0)./N
    copyto!(X,[slope*t[k] for k = 1:N])
    X[1] .= x0
    X[end] .= xf
    return nothing
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
import CompactBasisFunctions: Basis, nodes, nnodes

@testset "$(rpad("Basis Tests",80))" begin

    struct BasisTest{T} <: Basis{T} end

    test_basis = BasisTest{Float64}()

    @test_throws ErrorException basis(test_basis)
    @test_throws ErrorException nodes(test_basis)

    @test_throws ErrorException nbasis(test_basis)
    @test_throws ErrorException nnodes(test_basis)

    @test_throws ErrorException order(test_basis)
    @test_throws ErrorException degree(test_basis)

    @test_throws ErrorException eachindex(test_basis)

end


text:
module GitHub

import JSON
using Compat, Compat.Pkg

const AUTH_NOTE = "Julia Package Manager"
const AUTH_DATA = Dict{Any,Any}(
    "scopes" => ["repo"],
    "note" => "$AUTH_NOTE: $(gethostname())",
    "note_url" => "http://docs.julialang.org/en/latest/manual/packages/",
)

function user()
    usr = LibGit2.getconfig("github.user", "")
    if isempty(usr) #TODO: add `config` command to Git REPL and change below info
        throw(Pkg.PkgError("""
        no GitHub user name configured; please configure with:

            PkgDev.config()

        """))
    end
    return usr
end

function curl(url::AbstractString, opts::Cmd=``)
    success(`curl --version`) || throw(Pkg.PkgError("using the GitHub API requires having `curl` installed"))
    out, proc = open(`curl -i -s -S $opts $url`,"r")
    head = readline(out)
    status = parse(Int,split(head,r"\s+";limit=3)[2])
    header = Dict{AbstractString,AbstractString}()
    for line in eachline(out)
        if !occursin(r"^\s*$", line)
            (k,v) = split(line, r":\s*"; limit=2)
            header[k] = v
            continue
        end
        wait(proc); return status, header, String(read(out))
    end
    throw(Pkg.PkgError("strangely formatted HTTP response"))
end
curl(url::AbstractString, data::Nothing, opts::Cmd=``) = curl(url,opts)
curl(url::AbstractString, data, opts::Cmd=``) =
    curl(url,`--data $(sprint(io->JSON.print(io,data))) $opts`)

function delete_token()
    tokfile = Pkg.Dir.path(".github","token")
    Base.rm(tokfile)
    Compat.@info("Could not authenticate with existing token. Deleting token and trying again.")
end

readtoken(tokfile=Pkg.Dir.path(".github","token")) = isfile(tokfile) ? strip(readchomp(tokfile)) : ""

function token(user::AbstractString=user())
    tokfile = Pkg.Dir.path(".github","token")
    tok = readtoken(tokfile)
    !isempty(tok) && return tok

    Compat.@info("""
Creating a personal access token for Julia Package Manager on GitHub.
\tYou will be asked to provide credentials to your GitHub account.""")
    params = merge(AUTH_DATA, Dict("fingerprint" => randstring(40)))
    status, header, content = curl("https://api.github.com/authorizations",params,`-u $user`)
    tfa = false

    # Check for two-factor authentication
    if status == 401 && get(header, "X-GitHub-OTP", "") |> x->startswith(x, "required") && isinteractive()
        tfa = true
        Compat.@info("Two-factor authentication in use.  Enter auth code.  (You may have to re-enter your password.)")
        print(STDERR, "Authentication code: ")
        code = readline(STDIN) |> chomp
        status, header, content = curl("https://api.github.com/authorizations",params,`-H "X-GitHub-OTP: $code" -u $user`)
    end

    if status == 422
        error_code = JSON.parse(content)["errors"][1]["code"]
        throw(Pkg.PkgError("GitHub returned validation error (422): $error_code: $(JSON.parse(content)["message"])"))
    else
        (status != 401 && status != 403) || throw(Pkg.PkgError("$status: $(JSON.parse(content)["message"])"))
        tok = JSON.parse(content)["token"]
    end

    mkpath(dirname(tokfile))
    open(io->println(io,tok),tokfile,"w")
    chmod(tokfile, 0o600)
    return tok
end

function req(resource::AbstractString, data, opts::Cmd=``)
    url = "https://api.github.com/$resource"
    status, header, content = curl(url,data,`-u $(token()):x-oauth-basic $opts`)
    if (status == 302) || (status == 307) # Temporary redirect
        url = chomp(header["Location"])
        status, header, content = curl(url,data,`-u $(token()):x-oauth-basic $opts`)
    end
    response = JSON.parse(content)
    status, response
end

GET(resource::AbstractString, data, opts::Cmd=``) = req(resource,data,opts)
HEAD(resource::AbstractString, data, opts::Cmd=``) = req(resource,data,`-I $opts`)
PUT(resource::AbstractString, data, opts::Cmd=``) = req(resource,data,`-X PUT $opts`)
POST(resource::AbstractString, data, opts::Cmd=``) = req(resource,data,`-X POST $opts`)
PATCH(resource::AbstractString, data, opts::Cmd=``) = req(resource,data,`-X PATCH $opts`)
DELETE(resource::AbstractString, data, opts::Cmd=``) = req(resource,data,`-X DELETE $opts`)

for m in (:GET,:HEAD,:PUT,:POST,:PATCH,:DELETE)
    @eval $m(resource::AbstractString, opts::Cmd=``) = $m(resource,nothing,opts)
end

function pushable(owner::AbstractString, repo::AbstractString, user::AbstractString=user())
    status, response = HEAD("repos/$owner/$repo")
    status == 404 && throw(Pkg.PkgError("repo $owner/$repo does not exist"))
    status, response = GET("repos/$owner/$repo/collaborators/$user")
    status == 204 && return true
    status == 404 && return false
    throw(Pkg.PkgError("unexpected API status code: $status – $(response["message"])"))
end

function fork(owner::AbstractString, repo::AbstractString)
    status, response = POST("repos/$owner/$repo/forks")
    if status == 401
        delete_token()
        status, response = POST("repos/$owner/$repo/forks")
    end
    status == 202 || throw(Pkg.PkgError("forking $owner/$repo failed: $(response["message"])"))
    return response
end

function credentials()
    username = user()
    password = token(username)
    return Nullable(LibGit2.UserPasswordCredentials(username, password))
end

end # module


text:
@testset "maxintfloat $T" for T in (Double16, Double32, Double64)
    @test isinteger(maxintfloat(T))

    # Previous integer is representable, next integer is not
    @test maxintfloat(T) == (maxintfloat(T) - one(T)) + one(T)
    @test maxintfloat(T) != (maxintfloat(T) + one(T)) - one(T)

    @test maxintfloat(T) < floatmax(T)
end

@testset "Inf and NaN generation $T" for T in (Double16, Double32, Double64)
    @test isinf(T(Inf)) == isinf(inf(T))
    @test isinf(T(Inf)) == isposinf(posinf(T))
    @test isinf(T(-Inf)) == isneginf(neginf(T))
    @test isnan(T(NaN)) == isnan(nan(T))

    @test isinf(T(Inf32)) == isinf(inf(T))
    @test isinf(T(Inf32)) == isposinf(posinf(T))
    @test isinf(T(-Inf32)) == isneginf(neginf(T))
    @test isnan(T(NaN32)) == isnan(nan(T))

    @test isinf(T(Inf16)) == isinf(inf(T))
    @test isinf(T(Inf16)) == isposinf(posinf(T))
    @test isinf(T(-Inf16)) == isneginf(neginf(T))
    @test isnan(T(NaN16)) == isnan(nan(T))
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module SISModel
using DifferentialEquations

# +
function sis_ode(du, u, parms, t)
    beta = parms[1]
    gamma = parms[2]
    susceptible = init[1]
    infected_ind = init[2]

    du[1] = -beta * susceptible * infected_ind + gamma * infected_ind
    du[2] = beta * susceptible * infected_ind - parms[2] * init[2]

end
# -

parms = [0.1,0.05]
init = [0.99,0.01,0.0]
tspan = (0.0,200.0)

sis_prob = ODEProblem(sis_ode, init, tspan,parms)

sis_sol = solve(sis_prob);

#Visualization
using Plots
plot(sis_sol,xlabel="Time",ylabel="Number")

end


text:
include("plot_utils.jl")

function BmFig(ρ, u, f, ff_1, ff_2)
    DetectedPopulation = 10^9
    ft0 = 1 - ff_1
    ft1 = ff_1 - ff_2
    ft2 = ff_2
    if ft0 <= 0
        println("Wrong selected frequencies. Exiting....")
        return
    end

    N(s, ρ, u) = exp((ρ / u - 1 / u) * s) + 2 * sinh((1 / u - ρ / u) * s) / ρ

    td0 = findtd(u, ρ, DetectedPopulation * ft0) # time till the number of exclusively type-0 cells are at least DetectedPopulation*ft0
    ts1 = Int.(floor(td0 / 3)) # time of appearance of type-1 cells
    ts2 = Int.(ceil(2 * td0 / 3)) # time of appearance of type-2 cells
    s_1 = finds_i(u, ρ, td0 - ts1, DetectedPopulation * ft1) # calculates the selective growth advantage necessary for the new type to reach
    # at least DetectedPopulation*ft1 number of cells after td0 - ts1 time.
    # this is a helper function to allow us to select the values that make the pattern visually apparent for the example.
    ρ1 = ρ / (s_1 + 1)
    s_2 = finds_i(u, ρ1, td0 - ts2, DetectedPopulation * ft2)
    ρ2 = ρ1 / (s_2 + 1)
#     @show td0, ts1, s_1, ts2, s_2, ρ1, ρ2
    overallPopulation = (N(td0, ρ, u) + N(td0 - ts1, ρ1, u) + N(td0 - ts2, ρ2, u))
    type1Ratio = N(td0 - ts1, ρ1, u) / overallPopulation
    type2Ratio = N(td0 - ts2, ρ2, u) / overallPopulation

    println("========================================SUMMARY=========================================== ")
    println(
        "The growth parameters of the example are:\n\t ρ = ",
        ρ,
        "\n\t mutation probability: u = ",
        u,
        "\n\t detectability frequency: f = ",
        f,
        "\n\t rho of type 1: ρ_1 = ",
        ρ1,
        "\n\t rho of type 2: ρ_2 = ",
        ρ2,
        "",
    )

    println("The population at detection time is ", overallPopulation)
    println(
        "The type-1 mutation appears in ",
        100 * (type1Ratio + type2Ratio),
        "% of the population",
    )
    println("The type-2 mutation appears in ", 100 * type2Ratio, "% of the population")
    println("========================================================================================== \n\n")

    # Create the 0-->ts1 part of Bm
    xs, ys, error, options = BmMutIndexPlotPrepare(f / ft0, td0, u, ρ)

    # Create the ts1-->ts2 part of Bm
    xs1, ys1, error1, options1 = BmMutIndexPlotPrepare(f / ft1, td0 - ts1, u, ρ1)

    # Create the ts1-->detection part of the plot
    xs2, ys2, error2, options2 = BmMutIndexPlotPrepare(f / ft2, td0 - ts2, u, ρ2)

    # @show length(ys[1]), length(ys1[1]), length(ys2[1]), td0-ts2

    xs_pop = [popEstimate(s, u, ρ) for s = 1:td0]
    [xs_pop[s] += popEstimate(s - ts1, u, ρ1) for s = ts1:ts2]
    [xs_pop[s] += popEstimate(s - ts2, u, ρ2) for s = ts2:td0]

    # Probabilities are a function of the population

    # Construct final ys
    _y0 = deepcopy(ys[1])
    _y1 = deepcopy(ys1[1])
    _y2 = deepcopy(ys2[1])
    # after the first type is introduced
    #p0, p1, p2 = calculate_sampling_probabilities(u, ts1, ts2, td0, ρ, ρ1, ρ2 )
    p0, p1, p2 = calculate_sampling_probabilities_populations(u, ts1, ts2, td0, ρ, ρ1, ρ2)
    h_timeseries = plotQuick(1:length(p0), [p0, p1, p2], "Sampling Probabilities over Time")
    [_y0[s+ts1] = p1[s+ts1] * _y1[s] + p0[s+ts1] * _y0[s+ts1] for s = 1:(td0-ts1)]
    # after the second type is introduced
    [
        _y0[ts2+s] =
            p0[ts2+s] * _y0[ts2+s] + p1[ts2+s] * _y1[(ts2-ts1)+s] + p2[ts2+s] * _y2[s]
        for s = 1:td0-ts2
    ]
    y_out = [_y0]
    _m = deepcopy(_y0)

    ## Seeding randomly from the living population  (3rd strategy. See paper for details)
    title = "Seeding randomly from the living population"
    error = zeros(length(_y0))
    [
        error[s+ts1] = sqrt(
            (_m[s+ts1] - _y0[s+ts1])^2 * p0[s+ts1] + (_m[s+ts1] - _y1[s])^2 * p1[s+ts1],
        ) for s = 1:(td0-ts1)
    ]
    [
        error[s+ts2] = sqrt(
            (_m[s+ts2] - _y0[s+ts2])^2 * p0[s+ts2] +
            (_m[s+ts2] - _y1[s+ts2-ts1])^2 * p1[s+ts2] +
            (_m[s+ts2] - _y2[s])^2 * p2[s+ts2],
        ) for s = 1:td0-ts2
    ]
    axs1 = plotBmPop(xs_pop, y_out, error, title)

    ## Only New type can Seed
    title = "Seeding from the most advanced clone"
    # Construct final ys
    _y0 = deepcopy(ys[1])
    _y1 = deepcopy(ys1[1])
    _y2 = deepcopy(ys2[1])
    # after the first type is introduced
    p1 = 1
    p0 = 0
    [_y0[s+ts1] = p1 * _y1[s] + p0 * _y0[s+ts1] for s = 1:td0-ts1]
    # after the second type is introduced
    p0 = 0
    p1 = 0
    p2 = 1
    [_y0[ts2+s] = p0 * _y0[ts2+s] + p1 * _y1[(ts2-ts1)+s] + p2 * _y2[s] for s = 1:(td0-ts2)]
    y_out = [_y0]
    _m = deepcopy(_y0)
    error = zeros(length(_y0))
    [
        error[s+ts1] = sqrt((_m[s+ts1] - _y0[s+ts1])^2 * p0 + (_m[s+ts1] - _y1[s])^2 * p1)
        for s = 1:ts2-ts1
    ]
    [
        error[s+ts2] = sqrt(
            (_m[s+ts2] - _y0[s+ts2])^2 * p0 +
            (_m[s+ts2] - _y1[s+ts2-ts1])^2 * p1 +
            (_m[s+ts2] - _y2[s])^2 * p2,
        ) for s = 1:td0-ts2
    ]
    axs2 = plotBmPop(xs_pop, y_out, error, title)
    # @show _m[ts2+1:end] - _y2

    ## Probabilities Proportional to fitness
    title = "Seeding proportional to subclone fitness"
    # Construct final ys
    _y0 = deepcopy(ys[1])
    _y1 = deepcopy(ys1[1])
    _y2 = deepcopy(ys2[1])
    # after the first type is introduced
    p1 = (1 + s_1) / (2 + s_1)
    p0 = (1) / (2 + s_1)
    [_y0[s+ts1] = p1 * _y1[s] + p0 * _y0[s+ts1] for s = 1:(td0-ts1)]
    # after the second type is introduced
    p0 = 1 / (2 + s_1 + (1 + s_1) * (1 + s_2))
    p1 = (1 + s_1) / (2 + s_1 + (1 + s_1) * (1 + s_2))
    p2 = (1 + s_1) * (1 + s_2) / (2 + s_1 + (1 + s_1) * (1 + s_2))
    [
        _y0[ts2+s] = p0 * _y0[ts2+s] + p1 * _y1[(ts2-ts1)+s] + p2 * _y2[s]
        for s = 1:length(_y2)
    ]
    y_out = [_y0]
    _m = deepcopy(_y0)
    error = zeros(length(_y0))
    [
        error[s+ts1] = sqrt((_m[s+ts1] - _y0[s+ts1])^2 * p0 + (_m[s+ts1] - _y1[s])^2 * p1)
        for s = 1:(td0-ts1)
    ]
    [
        error[s+ts2] = sqrt(
            (_m[s+ts2] - _y0[s+ts2])^2 * p0 +
            (_m[s+ts2] - _y1[s+ts2-ts1])^2 * p1 +
            (_m[s+ts2] - _y2[s])^2 * p2,
        ) for s = 1:length(_y2)
    ]
    axs3 = plotBmPop(xs_pop, y_out, error, title)

    h = @pgf GroupPlot(
        {
            group_style =
                {group_size = "3 by 1", horizontal_sep = "2cm", vertical_sep = "3cm"},
        },
        axs2,
        axs3,
        axs1,
    )

    return h, h_timeseries, axs1, axs2, axs3
end

# Z is the total number of mutations comprising Bp at detection time. It is assumed to be a measurable quantity
function BpFig(Z, ρ, u, f, ff1, ff2)
    DetectedPopulation = 10^9
    ft0 = 1 - ff_1
    ft1 = ff_1 - ff_2
    ft2 = ff_2
    if ft0 <= 0
        println("Wrong selected frequencies. Exiting....")
        return
    end
    f_bp = f
    f_bm = f
    N(s, ρ, u) = exp((ρ / u - 1 / u) * s) + 2 * sinh((1 / u - ρ / u) * s) / ρ

    td0 = findtd(u, ρ, DetectedPopulation * ft0) # time till the number of exclusively type-0 cells
                                                 # are at least DetectedPopulation*ft0
    ts1 = Int.(floor(td0 / 3)) # time of appearance of type-1 cells
    ts2 = Int.(ceil(2 * td0 / 3)) # time of appearance of type-2 cells
    s_1 = finds_i(u, ρ, td0 - ts1, DetectedPopulation * ft1) # calculates the selective growth advantage necessary for the new type to reach
    # at least DetectedPopulation*ft1 number of cells after td0 - ts1 time.
    # this is a helper function to allow us to select the values that make the pattern visually apparent for the example.
    ρ1 = ρ / (s_1 + 1)
    s_2 = finds_i(u, ρ1, td0 - ts2, DetectedPopulation * ft2)
    ρ2 = ρ1 / (s_2 + 1)
    # @show ft0, ft1, ft2
    overallPopulation = (N(td0, ρ, u) + N(td0 - ts1, ρ1, u) + N(td0 - ts2, ρ2, u))
    type1Ratio = N(td0 - ts1, ρ1, u) / overallPopulation
    type2Ratio = N(td0 - ts2, ρ2, u) / overallPopulation

    # Create the 0-->ts1 part of Bm
#     @show f, ft0 ,ft1, ft2
    xs, ys, error, options = BmMutIndexPlotPrepare(f / ft0, td0, u, ρ)

    # Create the ts1-->ts2 part of Bm
    xs1, ys1, error1, options1 = BmMutIndexPlotPrepare(f / ft1, td0 - ts1, u, ρ1)

    # Create the ts1-->detection part of the plot
    xs2, ys2, error2, options2 = BmMutIndexPlotPrepare(f / ft2, td0 - ts2, u, ρ2)

    # @show length(ys[1]), length(ys1[1]), length(ys2[1]), td0-ts2

    xs_pop = [popEstimate(s, u, ρ) for s = 1:td0]
    [xs_pop[s] += popEstimate(s - ts1, u, ρ1) for s = ts1:ts2]
    [xs_pop[s] += popEstimate(s - ts2, u, ρ2) for s = ts2:td0]

    ## Seeding randomly from the population
    title = "Seeding randomly from the population"
    # Construct final ys
    _y0 = deepcopy(ys[1])
    _y1 = deepcopy(ys1[1])
    _y2 = deepcopy(ys2[1])
    # after the first type is introduced
    p0, p1, p2 = calculate_sampling_probabilities_populations(u, ts1, ts2, td0, ρ, ρ1, ρ2)

    h_timeseries = plotQuick(1:length(p0), [p0, p1, p2], "Sampling Probabilities over Time")
    [_y0[s+ts1] = p1[s+ts1] * _y1[s] + p0[s+ts1] * _y0[s+ts1] for s = 1:(td0-ts1)]
    # after the second type is introduced
    [
        _y0[ts2+s] =
            p0[ts2+s] * _y0[ts2+s] + p1[ts2+s] * _y1[(ts2-ts1)+s] + p2[ts2+s] * _y2[s]
        for s = 1:td0-ts2
    ]
    y_out = [_y0]
    _m = deepcopy(_y0)
    error = zeros(length(_y0))
    [
        error[s+ts1] = sqrt(
            (_m[s+ts1] - _y0[s+ts1])^2 * p0[s+ts1] + (_m[s+ts1] - _y1[s])^2 * p1[s+ts1],
        ) for s = 1:(td0-ts1)
    ]
    [
        error[s+ts2] = sqrt(
            (_m[s+ts2] - _y0[s+ts2])^2 * p0[s+ts2] +
            (_m[s+ts2] - _y1[s+ts2-ts1])^2 * p1[s+ts2] +
            (_m[s+ts2] - _y2[s])^2 * p2[s+ts2],
        ) for s = 1:td0-ts2
    ]
    axs1 = plotBmPop(xs_pop, y_out, error, title)

    # Calculate Bp
    bpys = [
        [_m[k_seed] - k_seed + Z - 1 for k_seed = 1:td0]
    ]
    Bpaxs1 = plotBp(xs_pop, bpys, error, title)

    ## Only New type can Seed
    title = "Seeding from the most advanced clone"
    # Construct final ys
    _y0 = deepcopy(ys[1])
    _y1 = deepcopy(ys1[1])
    _y2 = deepcopy(ys2[1])
    # after the first type is introduced
    p1 = 1
    p0 = 0
    [_y0[s+ts1] = p1 * _y1[s] + p0 * _y0[s+ts1] for s = 1:td0-ts1]
    # after the second type is introduced
    p0 = 0
    p1 = 0
    p2 = 1
    [_y0[ts2+s] = p0 * _y0[ts2+s] + p1 * _y1[(ts2-ts1)+s] + p2 * _y2[s] for s = 1:(td0-ts2)]
    y_out = [_y0]
    _m = deepcopy(_y0)
    error = zeros(length(_y0))
    [
        error[s+ts1] = sqrt((_m[s+ts1] - _y0[s+ts1])^2 * p0 + (_m[s+ts1] - _y1[s])^2 * p1)
        for s = 1:ts2-ts1
    ]
    [
        error[s+ts2] = sqrt(
            (_m[s+ts2] - _y0[s+ts2])^2 * p0 +
            (_m[s+ts2] - _y1[s+ts2-ts1])^2 * p1 +
            (_m[s+ts2] - _y2[s])^2 * p2,
        ) for s = 1:td0-ts2
    ]
    axs2 = plotBmPop(xs_pop, y_out, error, title)

    # Calculate Bp
    bpys = [
        [_m[k_seed] - k_seed + Z - 1 for k_seed = 1:td0]
    ]
    Bpaxs2 = plotBp(xs_pop, bpys, error, title)

    # Probabilities Proportional to fitness
    title = "Seeding proportional to subclone fitness"
    # Construct final ys
    _y0 = deepcopy(ys[1])
    _y1 = deepcopy(ys1[1])
    _y2 = deepcopy(ys2[1])
    # after the first type is introduced
    p1 = (1 + s_1) / (2 + s_1)
    p0 = (1) / (2 + s_1)
    [_y0[s+ts1] = p1 * _y1[s] + p0 * _y0[s+ts1] for s = 1:(td0-ts1)]
    # after the second type is introduced
    p0 = 1 / (2 + s_1 + (1 + s_1) * (1 + s_2))
    p1 = (1 + s_1) / (2 + s_1 + (1 + s_1) * (1 + s_2))
    p2 = (1 + s_1) * (1 + s_2) / (2 + s_1 + (1 + s_1) * (1 + s_2))
    [
        _y0[ts2+s] = p0 * _y0[ts2+s] + p1 * _y1[(ts2-ts1)+s] + p2 * _y2[s]
        for s = 1:length(_y2)
    ]
    y_out = [_y0]
    _m = deepcopy(_y0)

    error = zeros(length(_y0))
    [
        error[s+ts1] = sqrt((_m[s+ts1] - _y0[s+ts1])^2 * p0 + (_m[s+ts1] - _y1[s])^2 * p1)
        for s = 1:(td0-ts1)
    ]
    [
        error[s+ts2] = sqrt(
            (_m[s+ts2] - _y0[s+ts2])^2 * p0 +
            (_m[s+ts2] - _y1[s+ts2-ts1])^2 * p1 +
            (_m[s+ts2] - _y2[s])^2 * p2,
        ) for s = 1:length(_y2)
    ]
    axs3 = plotBmPop(xs_pop, y_out, error, title)

    # Calculate Bp
    bpys = [
        [_m[k_seed] - k_seed + Z - 1 for k_seed = 1:td0]
    ]
    Bpaxs3 = plotBp(xs_pop, bpys, error, title)

    h = @pgf GroupPlot(
        {
            group_style =
                {group_size = "3 by 2", horizontal_sep = "2cm", vertical_sep = "3cm"},
        },
        axs1,
        axs3,
        axs2,
        Bpaxs1,
        Bpaxs3,
        Bpaxs2,
    )

    return h, h_timeseries, Bpaxs1, Bpaxs2, Bpaxs3
end

function plotBmNeutral(xs, ys, title)
    ggplot2_axis_theme = @pgf {
    }

    ggplot2_plot_theme =
        @pgf {mark = "square", mark_size = 0, mark_options = "solid", line_width = "1.64pt"}

    n = length(ys)

    # Evenly spread out colors
    colors = range(HSL(colorant"red"), stop=HSL(colorant"green"), length=n)
    @pgf SemiLogXAxis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = minimum(minimum(ys)),
            ymax = maximum(maximum(ys)),
            title = title,
            xlabel = "Seeding time",
            ylabel = L"$\mathbb{E}[B^{\mathit{neutral}}_{md}]$",
            legend_pos = "north west",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
        Legend([L"$f = 0.005$", L"$f = 0.01$", L"$f = 0.05$", L"$f = 0.1$"]),
    )
end

max_stars_repo_name: SunPathLab/Comet

max_issues_repo_name: SunPathLab/Comet

max_forks_repo_name: SunPathLab/Comet


text:
#=
Computes a sequence of marginal densities for a continuous state space
Markov chain :math:`X_t` where the transition probabilities can be represented
as densities. The estimate of the marginal density of X_t is

    1/n sum_{i=0}^n p(X_{t-1}^i, y)

This is a density in y.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-01

References
----------

http://quant-econ.net/jl/stationary_densities.html
=#

"""
A look ahead estimator associated with a given stochastic kernel p and a vector
of observations X.

##### Fields

- `p::Function`: The stochastic kernel. Signature is `p(x, y)` and it should be
vectorized in both inputs
- `X::Matrix`: A vector containing observations. Note that this can be passed as
any kind of `AbstractArray` and will be coerced into an `n x 1` vector.

"""
type LAE
    p::Function
    X::Matrix

    function LAE(p::Function, X::AbstractArray)
        n = length(X)
        new(p, reshape(X, n, 1))
    end
end

"""
A vectorized function that returns the value of the look ahead estimate at the
values in the array y.

##### Arguments

- `l::LAE`: Instance of `LAE` type
- `y::Array`: Array that becomes the `y` in `l.p(l.x, y)`

##### Returns

- `psi_vals::Vector`: Density at `(x, y)`

"""
function lae_est{T}(l::LAE, y::AbstractArray{T})
    k = length(y)
    v = l.p(l.X, reshape(y, 1, k))
    psi_vals = mean(v, 1)
    return squeeze(psi_vals, 1)
end


text:
@testset "ReifiedConstraint" begin
    @testset "Basic ==" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([3, 4]))
        @variable(m, b, Bin)
        @constraint(m, b := {x + y == 7})
        @objective(m, Max, b)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 1.0
        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value(x) ≈ 4
        @test JuMP.value(y) ≈ 3
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "Basic >=" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 3]))
        @variable(m, b, Bin)
        @constraint(m, b := {x + y >= 6.1})
        @objective(m, Max, b)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 1.0
        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value(x) ≈ 4
        @test JuMP.value(y) ≈ 3
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "Basic >=" begin
        m = Model(CSJuMPTestOptimizer(; branch_strategy = :ABS))
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 3]))
        @variable(m, b, Bin)
        # missing { }
        @test_macro_throws ErrorException begin
            @constraint(m, b := x + y >= 6.1)
        end

        # no vectorization support
        @test_macro_throws ErrorException begin
            @constraint(m, b := {[x, y] .>= 6.1})
        end
    end

    @testset "Basic == where not active" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4, 5]))
        @variable(m, y, CS.Integers([2, 4]))
        @variable(m, b, Bin)
        @constraint(m, b := {x + y == 7})
        @objective(m, Max, b + x + y)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 9.0
        @test JuMP.value(b) ≈ 0.0
        @test JuMP.value(x) ≈ 5
        @test JuMP.value(y) ≈ 4
        com = CS.get_inner_model(m)
        @test !com.constraints[1].activator_in_inner
        @test is_solved(com)
    end

    @testset "Alldifferent" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 3, 4]))
        @variable(m, b, Bin)
        @constraint(m, b := {[x, y] in CS.AllDifferent()})
        @objective(m, Max, 0.9 * b + x + y)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 8.0
        @test JuMP.value(b) ≈ 0.0
        @test JuMP.value(x) ≈ 4
        @test JuMP.value(y) ≈ 4
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "Alldifferent minimize" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([1, 2, 3, 4]))
        @variable(m, b, Bin)
        @constraint(m, !b := {[x, y] in CS.AllDifferent()})
        @objective(m, Min, 0.9 * b + x + y)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 2.9
        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value(x) ≈ 1
        @test JuMP.value(y) ≈ 1
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "Alldifferent where active" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 4]))
        @variable(m, b, Bin)
        @constraint(m, b := {[x, y] in CS.AllDifferent()})
        @objective(m, Max, 5b + x + y)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 11.0
        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value(x) ≈ 2 || JuMP.value(y) ≈ 2
        @test JuMP.value(x) ≈ 4 || JuMP.value(y) ≈ 4
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "TableConstraint where active" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "keep_logs" => true,
            "logging" => [],
        ))
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 4]))
        @variable(m, 1 <= z <= 10, Int)
        @variable(m, b, Bin)
        @constraint(
            m,
            b := {
                [x, y, z] in CS.TableSet([
                    1 2 3
                    1 4 3
                    2 4 3
                    3 4 3
                    1 2 4
                    1 4 4
                    2 4 4
                    3 4 4
                ]),
            }
        )
        @objective(m, Max, 5b + x + y)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 5 + 2 + 4
        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value(x) ≈ 2
        @test JuMP.value(y) ≈ 4
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "TableConstraint where inactive" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "keep_logs" => true,
            "logging" => [],
        ))
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 4]))
        @variable(m, b, Bin)
        @constraint(m, !b := {[x, y] in CS.TableSet([
            1 2
            2 4
            3 4
            4 4
        ])})
        @objective(m, Max, 5b + x + y)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 5 + 4 + 2
        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value(x) ≈ 4
        @test JuMP.value(y) ≈ 2
        com = CS.get_inner_model(m)
        @test is_solved(com)
    end

    @testset "TableConstraint with optimization" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => [],
            "lp_optimizer" => cbc_optimizer,
            "keep_logs" => true,
        ))
        @variable(m, 1 <= a <= 100, Int)
        @variable(m, 1 <= b <= 100, Int)
        @variable(m, 1 <= c <= 100, Int)
        @variable(m, reified, Bin)
        table = zeros(Int, (52, 3))
        r = 1
        for i in 1:100, j in i:100, k in j:100
            if i^2 + j^2 == k^2
                table[r, :] .= [i, j, k]
                r += 1
            end
        end
        @constraint(m, reified := {[a, b, c] in CS.TableSet(table)})
        @objective(m, Max, 5 * reified + b + c)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 5 + 96 + 100
        @test JuMP.value(a) ≈ 28
        @test JuMP.value(b) ≈ 96
        @test JuMP.value(c) ≈ 100
        @test JuMP.value(reified) ≈ 1
        com = CS.get_inner_model(m)
        @test is_solved(com)
        logs = CS.get_logs(com)
        @test CS.sanity_check_log(logs[:tree])
    end

    @testset "all different != 0 (Issue 202)" begin
        n = 2
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_solutions" => true,
            "logging" => [],
        ))
        @variable(model, 0 <= x[1:n] <= n, Int)
        b_len = length([1 for i in 2:n for j in 1:(i - 1) for k in 1:3])
        @variable(model, bs[1:b_len], Bin)
        c = 1
        for i in 2:n, j in 1:(i - 1)
            # b1: bs[c]
            # b2: bs[c+1]
            # b3: bs[c+2]
            @constraint(model, bs[c] := {x[i] != 0})
            @constraint(model, bs[c + 1] := {x[j] != 0})
            @constraint(model, bs[c + 2] := {bs[c] + bs[c + 1] == 2})
            @constraint(model, bs[c + 2] => {x[i] != x[j]})
            c += 3
        end

        optimize!(model)

        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        @test num_sols == 7

        xx_set = Set()
        for sol in 1:num_sols
            xx = convert.(Integer, JuMP.value.(x, result = sol))
            bss = convert.(Integer, JuMP.value.(bs, result = sol))
            push!(xx_set, (xx[1], xx[2]))
        end
        @test length(xx_set) == num_sols
        @test (0, 0) in xx_set
        @test (0, 1) in xx_set
        @test (0, 2) in xx_set
        @test (1, 0) in xx_set
        @test (1, 2) in xx_set
        @test (2, 0) in xx_set
        @test (2, 1) in xx_set
    end
end


text:
using Test
using VertexFinder: Queue, enqueue!, dequeue!, isempty

q = Queue(Int, 1)
@test isempty(q) == true

enqueue!(q, 1)
enqueue!(q, 2)
@test isempty(q) == false

@test dequeue!(q) == 1
@test dequeue!(q) == 2
@test isempty(q) == true

@test_throws Exception dequeue!(q)


text:
function measure_mean!(tuning_run::Run, x::Configuration)
    configurations = Array{Configuration}(undef, tuning_run.cost_evaluations)
    fill!(configurations, deepcopy(x))

    pmap_cost(x::Configuration) = tuning_run.cost(x, tuning_run.cost_arguments)
    results = pmap(pmap_cost, configurations)

    for i = 1:tuning_run.cost_evaluations
        tuning_run.cost_values[i] = results[i]
    end

    mean(results)
end

function sequential_measure_mean!(tuning_run::Run, x::Configuration)
    for i = 1:tuning_run.cost_evaluations
        tuning_run.cost_values[i] = tuning_run.cost(x, tuning_run.cost_arguments)
    end
    mean(tuning_run.cost_values)
end


text:
"""
    near_field_forces!(properties, surfaces, wakes, reference, freestream, Γ;
        dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, surface_id,
        wake_finite_core, wake_shedding_locations, trailing_vortices, xhat)

Calculate local panel forces in the body frame.
"""
@inline function near_field_forces!(props, surfaces, wakes, ref, fs, Γ;
    dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, surface_id,
    wake_finite_core, wake_shedding_locations, trailing_vortices, xhat)

    # number of surfaces
    nsurf = length(surfaces)

    # loop through receiving surfaces
    iΓ = 0 # index for accessing Γ
    for isurf = 1:nsurf

        receiving = surfaces[isurf]
        nr = length(receiving)
        nr1, nr2 = size(receiving)
        cr = CartesianIndices(receiving)

        # loop through receiving panels
        for i = 1:length(receiving)

            # get panel cartesian index
            I = cr[i]

            # --- Calculate forces on the panel bound vortex --- #

            # bound vortex location
            rc = top_center(receiving[i])

            # freestream velocity
            Vi = freestream_velocity(fs)

            # rotational velocity
            Vi += rotational_velocity(rc, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                Vi += additional_velocity(rc)
            end

            # velocity due to surface motion
            if !isnothing(Vh)
                Vi += Vh[isurf][i]
            end

            # induced velocity from surfaces and wakes
            jΓ = 0 # index for accessing Γ
            for jsurf = 1:nsurf

                # number of panels on sending surface
                sending = surfaces[jsurf]
                Ns = length(sending)

                # see if wake panels are being used
                wake_panels = nwake[jsurf] > 0

                # check if we need to shift shedding locations
                if isnothing(wake_shedding_locations)
                    shedding_locations = nothing
                else
                    shedding_locations = wake_shedding_locations[jsurf]
                end

                # extract circulation values corresonding to the sending surface
                vΓ = view(Γ, jΓ+1:jΓ+Ns)

                # induced velocity from this surface
                if isurf == jsurf
                    # induced velocity on self
                    Vi += induced_velocity(I, surfaces[jsurf], vΓ;
                        finite_core = surface_id[isurf] != surface_id[jsurf],
                        wake_shedding_locations = shedding_locations,
                        symmetric = symmetric[jsurf],
                        trailing_vortices = trailing_vortices[jsurf] && !wake_panels,
                        xhat = xhat)
                else
                    # induced velocity on another surface
                    Vi += induced_velocity(rc, surfaces[jsurf], vΓ;
                        finite_core = surface_id[isurf] != surface_id[jsurf],
                        wake_shedding_locations = shedding_locations,
                        symmetric = symmetric[jsurf],
                        trailing_vortices = trailing_vortices[jsurf] && !wake_panels,
                        xhat = xhat)
                end

                # induced velocity from corresponding wake
                if wake_panels
                    Vi += induced_velocity(rc, wakes[jsurf];
                        finite_core = wake_finite_core[jsurf] || (surface_id[isurf] != surface_id[jsurf]),
                        symmetric = symmetric[jsurf],
                        nc = nwake[jsurf],
                        trailing_vortices = trailing_vortices[jsurf],
                        xhat = xhat)
                end

                jΓ += Ns # increment Γ index for sending panels
            end

            # steady part of Kutta-Joukowski theorem
            Γi = I[1] == 1 ? Γ[iΓ+i] : Γ[iΓ+i] - Γ[iΓ+i-1] # net circulation
            Δs = top_vector(receiving[I]) # bound vortex vector
            tmp = cross(Vi, Δs)
            Fbi = RHO*Γi*tmp

            if !isnothing(dΓdt)
                # unsteady part of Kutta-Joukowski theorem

                #TODO: decide whether to divide by perpindicular velocity like
                # Drela does in ASWING?

                dΓdti = I[1] == 1 ? dΓdt[iΓ+i] : (dΓdt[iΓ+i] + dΓdt[iΓ+i-1])/2
                c = receiving[I].chord
                Fbi += RHO*dΓdti*c*tmp

            end

            # --- Calculate forces on the left bound vortex --- #

            # bound vortex location
            rc = left_center(receiving[I])

            # freestream velocity
            Veff = freestream_velocity(fs)

            # rotational velocity
            Veff += rotational_velocity(rc, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                Veff += additional_velocity(rc)
            end

            # velocity due to surface motion
            if !isnothing(Vv)
                Veff += Vv[isurf][I[1], I[2]]
            end

            # NOTE: We don't include induced velocity in the effective velocity
            # for the vertical segments because its influence is likely negligible
            # once we take the cross product with the bound vortex vector. This
            # is also assumed in AVL. This could change in the future.

            # steady part of Kutta-Joukowski theorem
            Γli = Γ[iΓ+i]
            Δs = left_vector(receiving[I])
            Fbli = RHO*Γli*cross(Veff, Δs)

            # --- Calculate forces on the right bound vortex --- #

            rc = right_center(receiving[I])

            # freestream velocity
            Veff = freestream_velocity(fs)

            # rotational velocity
            Veff += rotational_velocity(rc, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                Veff += additional_velocity(rc)
            end

            # velocity due to surface motion
            if !isnothing(Vv)
                Veff += Vv[isurf][I[1], I[2]+1]
            end

            # NOTE: We don't include induced velocity in the effective velocity
            # for the vertical segments because its influence is likely negligible
            # once we take the cross product with the bound vortex vector. This
            # is also assumed in AVL. This could change in the future.

            # steady part of Kutta-Joukowski theorem
            Γri = Γ[iΓ+i]
            Δs = right_vector(receiving[I])
            Fbri = RHO*Γri*cross(Veff, Δs)

            # store panel circulation, velocity, and forces
            q = 1/2*RHO*ref.V^2

            props[isurf][i] = PanelProperties(Γ[iΓ+i]/ref.V, Vi/ref.V,
                Fbi/(q*ref.S), Fbli/(q*ref.S), Fbri/(q*ref.S))
        end

        # increment Γ index for receiving panels
        iΓ += nr
    end

    return props
end

"""
    near_field_forces_derivatives!(properties, dproperties, surfaces, reference,
        freestream, Γ, dΓ; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake,
        surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat)

Version of [`near_field_forces!`](@ref) that also calculates the derivatives of
the local panel forces with respect to the freestream variables.
"""
near_field_forces_derivatives!

@inline function near_field_forces_derivatives!(props, dprops, surfaces, wakes,
    ref, fs, Γ, dΓ; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake,
    surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat)

    # unpack derivatives
    props_a, props_b, props_p, props_q, props_r = dprops
    Γ_a, Γ_b, Γ_p, Γ_q, Γ_r = dΓ

    # number of surfaces
    nsurf = length(surfaces)

    # loop through receiving surfaces
    iΓ = 0 # index for accessing Γ
    for isurf = 1:nsurf

        receiving = surfaces[isurf]
        nr = length(receiving)
        nr1, nr2 = size(receiving)
        cr = CartesianIndices(receiving)

        # loop through receiving panels
        for i = 1:length(receiving)

            I = cr[i]

            # --- Calculate forces on the panel bound vortex -- #
            rc = top_center(receiving[i])

            # freestream velocity
            Vi, dVi = freestream_velocity_derivatives(fs)
            Vi_a, Vi_b = dVi

            # rotational velocity
            Vrot, dVrot = rotational_velocity_derivatives(rc, fs, ref)
            Vi += Vrot
            Vi_p, Vi_q, Vi_r = dVrot

            # additional velocity field
            if !isnothing(additional_velocity)
                Vi += additional_velocity(rc)
            end

            # velocity due to surface motion
            if !isnothing(Vh)
                Vi += Vh[isurf][i]
            end

            # induced velocity from surfaces and wakes
            jΓ = 0 # index for accessing Γ
            for jsurf = 1:nsurf

                # number of panels on sending surface
                sending = surfaces[jsurf]
                Ns = length(sending)

                # see if wake panels are being used
                wake_panels = nwake[jsurf] > 0

                # check if we need to shift shedding locations
                if isnothing(wake_shedding_locations)
                    shedding_locations = nothing
                else
                    shedding_locations = wake_shedding_locations[jsurf]
                end

                # extract circulation values corresonding to the sending surface
                vΓ = view(Γ, jΓ+1:jΓ+Ns)

                vΓ_a = view(Γ_a, jΓ+1:jΓ+Ns)
                vΓ_b = view(Γ_b, jΓ+1:jΓ+Ns)
                vΓ_p = view(Γ_p, jΓ+1:jΓ+Ns)
                vΓ_q = view(Γ_q, jΓ+1:jΓ+Ns)
                vΓ_r = view(Γ_r, jΓ+1:jΓ+Ns)

                vdΓ = (vΓ_a, vΓ_b, vΓ_p, vΓ_q, vΓ_r)

                # induced velocity from this surface
                if isurf == jsurf
                    # induced velocity on self
                    Vind, dVind = induced_velocity_derivatives(I, surfaces[jsurf], vΓ, vdΓ;
                        finite_core = surface_id[isurf] != surface_id[jsurf],
                        wake_shedding_locations = shedding_locations,
                        symmetric = symmetric[jsurf],
                        trailing_vortices = trailing_vortices[jsurf] && !wake_panels,
                        xhat = xhat)
                else
                    # induced velocity on another surface
                    Vind, dVind = induced_velocity_derivatives(rc, surfaces[jsurf], vΓ, vdΓ;
                        finite_core = surface_id[isurf] != surface_id[jsurf],
                        wake_shedding_locations = shedding_locations,
                        symmetric = symmetric[jsurf],
                        trailing_vortices = trailing_vortices[jsurf] && !wake_panels,
                        xhat = xhat)
                end

                Vind_a, Vind_b, Vind_p, Vind_q, Vind_r = dVind

                Vi += Vind

                Vi_a += Vind_a
                Vi_b += Vind_b
                Vi_p += Vind_p
                Vi_q += Vind_q
                Vi_r += Vind_r

                # induced velocity from corresponding wake
                if wake_panels
                    Vi += induced_velocity(rc, wakes[jsurf];
                        finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf],
                        symmetric = symmetric[jsurf],
                        nc = nwake[jsurf],
                        trailing_vortices = trailing_vortices[jsurf],
                        xhat = xhat)
                end

                jΓ += Ns # increment Γ index for sending panels
            end

            # steady part of Kutta-Joukowski theorem
            if I[1] == 1
                Γi = Γ[iΓ+i]

                Γi_a = Γ_a[iΓ+i]
                Γi_b = Γ_b[iΓ+i]
                Γi_p = Γ_p[iΓ+i]
                Γi_q = Γ_q[iΓ+i]
                Γi_r = Γ_r[iΓ+i]
            else
                Γi = Γ[iΓ+i] - Γ[iΓ+i-1]

                Γi_a = Γ_a[iΓ+i] - Γ_a[iΓ+i-1]
                Γi_b = Γ_b[iΓ+i] - Γ_b[iΓ+i-1]
                Γi_p = Γ_p[iΓ+i] - Γ_p[iΓ+i-1]
                Γi_q = Γ_q[iΓ+i] - Γ_q[iΓ+i-1]
                Γi_r = Γ_r[iΓ+i] - Γ_r[iΓ+i-1]
            end

            # bound vortex vector
            Δs = top_vector(receiving[I])

            tmp = cross(Vi, Δs)

            Fbi = RHO*Γi*tmp

            Fbi_a = RHO*(Γi_a*tmp + Γi*cross(Vi_a, Δs))
            Fbi_b = RHO*(Γi_b*tmp + Γi*cross(Vi_b, Δs))
            Fbi_p = RHO*(Γi_p*tmp + Γi*cross(Vi_p, Δs))
            Fbi_q = RHO*(Γi_q*tmp + Γi*cross(Vi_q, Δs))
            Fbi_r = RHO*(Γi_r*tmp + Γi*cross(Vi_r, Δs))

            if !isnothing(dΓdt)
                # unsteady part of Kutta-Joukowski theorem

                #TODO: decide whether to divide by perpindicular velocity like
                # Drela does in ASWING?

                dΓdti = I[1] == 1 ? dΓdt[iΓ+i] : (dΓdt[iΓ+i] + dΓdt[iΓ+i-1])/2
                c = receiving[I].chord
                Fbi += RHO*dΓdti*c*tmp

            end

            # --- Calculate forces for the left bound vortex --- #

            rc = left_center(receiving[I])

            # freestream velocity
            Vfs, dVfs = freestream_velocity_derivatives(fs)
            Veff = Vfs
            Veff_a, Veff_b = dVfs

            # rotational velocity
            Vrot, dVrot = rotational_velocity_derivatives(rc, fs, ref)
            Veff += Vrot
            Veff_p, Veff_q, Veff_r = dVrot

            # additional velocity field
            if !isnothing(additional_velocity)
                Veff += additional_velocity(rc)
            end

            # velocity due to surface motion
            if !isnothing(Vv)
                Veff += Vv[isurf][I[1], I[2]]
            end

            # NOTE: We don't include induced velocity in the effective velocity
            # for the vertical segments because its influence is likely negligible
            # once we take the cross product with the bound vortex vector. This
            # is also assumed in AVL. This could change in the future.

            # steady part of Kutta-Joukowski theorem
            Γli = Γ[iΓ+i]

            Γli_a = Γ_a[iΓ+i]
            Γli_b = Γ_b[iΓ+i]
            Γli_p = Γ_p[iΓ+i]
            Γli_q = Γ_q[iΓ+i]
            Γli_r = Γ_r[iΓ+i]

            Δs = left_vector(receiving[I])

            tmp = cross(Veff, Δs)

            Fbli = RHO*Γli*tmp

            Fbli_a = RHO*(Γli_a*tmp + Γli*cross(Veff_a, Δs))
            Fbli_b = RHO*(Γli_b*tmp + Γli*cross(Veff_b, Δs))
            Fbli_p = RHO*(Γli_p*tmp + Γli*cross(Veff_p, Δs))
            Fbli_q = RHO*(Γli_q*tmp + Γli*cross(Veff_q, Δs))
            Fbli_r = RHO*(Γli_r*tmp + Γli*cross(Veff_r, Δs))

            # --- Calculate forces on the right bound vortex --- #

            rc = right_center(receiving[I])

            # freestream velocity
            Vfs, dVfs = freestream_velocity_derivatives(fs)
            Veff = Vfs
            Veff_a, Veff_b = dVfs

            # rotational velocity
            Vrot, dVrot = rotational_velocity_derivatives(rc, fs, ref)
            Veff += Vrot
            Veff_p, Veff_q, Veff_r = dVrot

            # additional velocity field
            if !isnothing(additional_velocity)
                Veff += additional_velocity(rc)
            end

            # velocity due to surface motion
            if !isnothing(Vv)
                Veff += Vv[isurf][I[1], I[2]+1]
            end

            # NOTE: We don't include induced velocity in the effective velocity
            # for the vertical segments because its influence is likely negligible
            # once we take the cross product with the bound vortex vector. This
            # is also assumed in AVL. This could change in the future.

            # steady part of Kutta-Joukowski theorem
            Γri = Γ[iΓ+i]

            Γri_a = Γ_a[iΓ+i]
            Γri_b = Γ_b[iΓ+i]
            Γri_p = Γ_p[iΓ+i]
            Γri_q = Γ_q[iΓ+i]
            Γri_r = Γ_r[iΓ+i]

            Δs = right_vector(receiving[I])

            tmp = cross(Veff, Δs)

            Fbri = RHO*Γri*tmp

            Fbri_a = RHO*(Γri_a*tmp + Γri*cross(Veff_a, Δs))
            Fbri_b = RHO*(Γri_b*tmp + Γri*cross(Veff_b, Δs))
            Fbri_p = RHO*(Γri_p*tmp + Γri*cross(Veff_p, Δs))
            Fbri_q = RHO*(Γri_q*tmp + Γri*cross(Veff_q, Δs))
            Fbri_r = RHO*(Γri_r*tmp + Γri*cross(Veff_r, Δs))

            # store panel circulation, velocity, and forces
            q = 1/2*RHO*ref.V^2

            props[isurf][I] = PanelProperties(Γ[iΓ+i]/ref.V, Vi/ref.V, Fbi/(q*ref.S),
                Fbli/(q*ref.S), Fbri/(q*ref.S))

            props_a[isurf][I] = PanelProperties(Γ_a[iΓ+i]/ref.V, Vi_a/ref.V, Fbi_a/(q*ref.S),
                Fbli_a/(q*ref.S), Fbri_a/(q*ref.S))
            props_b[isurf][I] = PanelProperties(Γ_b[iΓ+i]/ref.V, Vi_b/ref.V, Fbi_b/(q*ref.S),
                Fbli_b/(q*ref.S), Fbri_b/(q*ref.S))
            props_p[isurf][I] = PanelProperties(Γ_p[iΓ+i]/ref.V, Vi_p/ref.V, Fbi_p/(q*ref.S),
                Fbli_p/(q*ref.S), Fbri_p/(q*ref.S))
            props_q[isurf][I] = PanelProperties(Γ_q[iΓ+i]/ref.V, Vi_q/ref.V, Fbi_q/(q*ref.S),
                Fbli_q/(q*ref.S), Fbri_q/(q*ref.S))
            props_r[isurf][I] = PanelProperties(Γ_r[iΓ+i]/ref.V, Vi_r/ref.V, Fbi_r/(q*ref.S),
                Fbli_r/(q*ref.S), Fbri_r/(q*ref.S))
        end

        # increment Γ index for receiving panels
        iΓ += nr
    end

    return props, dprops
end

"""
    body_forces(system; kwargs...)

Return the body force coefficients given the panel properties for `surfaces`

Note that this function assumes that a near-field analysis has already been
performed to obtain the panel forces.

# Arguments
 - `system`: Object of type [`System`](@ref) which holds system properties

# Keyword Arguments
 - `frame`: frame in which to return `CF` and `CM`, options are [`Body()`](@ref) (default),
   [`Stability()`](@ref), and [`Wind()`](@ref)`
"""
function body_forces(system::System{TF}; frame = Body()) where TF

    @assert system.near_field_analysis[] "Near field analysis required"

    # unpack parameters stored in `system`
    surfaces = system.surfaces # surface panels defining each surface
    properties = system.properties # panel properties
    ref = system.reference[] # reference parameters
    fs = system.freestream[] # freestream parameters
    symmetric = system.symmetric # symmetric flag for each surface

    return body_forces(surfaces, properties, ref, fs, symmetric, frame)
end

"""
    body_forces(surfaces, properties, reference, freestream, symmetric; kwargs...)

Return the body force coefficients given the panel properties for `surfaces`

Note that this function assumes that a near-field analysis has already been
performed to obtain the panel forces.

# Arguments:
 - `surfaces`: Collection of surfaces, where each surface is represented by a
    matrix of surface panels (see [`SurfacePanel`](@ref)) of shape (nc, ns)
    where `nc` is the number of chordwise panels and `ns` is the number of
    spanwise panels
 - `properties`: Surface properties for each surface, where surface
    properties for each surface are represented by a matrix of panel properties
    (see [`PanelProperties`](@ref)) of shape (nc, ns) where `nc` is the number
    of chordwise panels and `ns` is the number of spanwise panels
 - `reference`: Reference parameters (see [`Reference`](@ref))
 - `freestream`: Freestream parameters (see [`Freestream`]@ref)
 - `symmetric`: (required) Flag for each surface indicating whether a mirror image
   (across the X-Z plane) was used when calculating induced velocities
 - `frame`: frame in which to return `CF` and `CM`, options are [`Body()`](@ref) (default),
   [`Stability()`](@ref), and [`Wind()`](@ref)
"""
function body_forces(surfaces, properties, ref, fs, symmetric, frame)

    TF = eltype(eltype(eltype(properties)))

    # initialize body force coefficients
    CF = @SVector zeros(TF, 3)
    CM = @SVector zeros(TF, 3)

    # loop through all surfaces
    for isurf = 1:length(surfaces)

        # initialize surface contribution to body force coefficients
        CFi = @SVector zeros(TF, 3)
        CMi = @SVector zeros(TF, 3)

        # loop through all panels on this surface
        for i = 1:length(surfaces[isurf])

            # top bound vortex
            rc = top_center(surfaces[isurf][i])
            Δr = rc - ref.r
            cf = properties[isurf][i].cfb
            CFi += cf
            CMi += cross(Δr, cf)

            # left bound vortex
            rc = left_center(surfaces[isurf][i])
            Δr = rc - ref.r
            cf = properties[isurf][i].cfl
            CFi += cf
            CMi += cross(Δr, cf)

            # right bound vortex
            rc = right_center(surfaces[isurf][i])
            Δr = rc - ref.r
            cf = properties[isurf][i].cfr
            CFi += cf
            CMi += cross(Δr, cf)
        end

        # adjust forces from this surface to account for symmetry
        if symmetric[isurf]
            CFi = SVector(2*CFi[1], 0.0, 2*CFi[3])
            CMi = SVector(0.0, 2*CMi[2], 0.0)
        end

        # add to body forces
        CF += CFi
        CM += CMi

    end

    # add reference length in moment normalization
    reference_length = SVector(ref.b, ref.c, ref.b)
    CM = CM ./ reference_length

    # switch to specified frame
    CF, CM = body_to_frame(CF, CM, ref, fs, frame)

    return CF, CM
end

"""
    body_forces_derivatives(system)

Return the body force coefficients for the `system` and their derivatives with
respect to the freestream variables

Note that this function assumes that a near-field analysis has already been
performed to obtain the panel forces.

# Arguments:
 - `system`: Object of type `System` which holds system properties
"""
@inline function body_forces_derivatives(system::System)

    # float number type
    TF = eltype(system)

    @assert system.near_field_analysis[] "Near field analysis required"
    @assert system.derivatives[] "Derivative computations required"

    # unpack parameters stored in `system`
    surfaces = system.surfaces # surface panels defining each surface
    ref = system.reference[] # reference parameters
    fs = system.freestream[] # freestream parameters
    symmetric = system.symmetric # symmetric flag for each surface
    properties = system.properties
    props_a, props_b, props_p, props_q, props_r = system.dproperties

    # initialize body force coefficients
    CF = @SVector zeros(TF, 3)
    CM = @SVector zeros(TF, 3)

    CF_a = @SVector zeros(TF, 3)
    CF_b = @SVector zeros(TF, 3)
    CF_p = @SVector zeros(TF, 3)
    CF_q = @SVector zeros(TF, 3)
    CF_r = @SVector zeros(TF, 3)

    CM_a = @SVector zeros(TF, 3)
    CM_b = @SVector zeros(TF, 3)
    CM_p = @SVector zeros(TF, 3)
    CM_q = @SVector zeros(TF, 3)
    CM_r = @SVector zeros(TF, 3)

    # loop through all surfaces
    for isurf = 1:length(surfaces)

        # initialize surface contribution to body force coefficients
        CFi = @SVector zeros(TF, 3)
        CMi = @SVector zeros(TF, 3)

        CFi_a = @SVector zeros(TF, 3)
        CFi_b = @SVector zeros(TF, 3)
        CFi_p = @SVector zeros(TF, 3)
        CFi_q = @SVector zeros(TF, 3)
        CFi_r = @SVector zeros(TF, 3)

        CMi_a = @SVector zeros(TF, 3)
        CMi_b = @SVector zeros(TF, 3)
        CMi_p = @SVector zeros(TF, 3)
        CMi_q = @SVector zeros(TF, 3)
        CMi_r = @SVector zeros(TF, 3)

        for i = 1:length(surfaces[isurf])

            # top bound vortex
            rc = top_center(surfaces[isurf][i])
            Δr = rc - ref.r
            cf = properties[isurf][i].cfb
            CFi += cf
            CMi += cross(Δr, cf)

            cf_a = props_a[isurf][i].cfb
            CFi_a += cf_a
            CMi_a += cross(Δr, cf_a)

            cf_b = props_b[isurf][i].cfb
            CFi_b += cf_b
            CMi_b += cross(Δr, cf_b)

            cf_p = props_p[isurf][i].cfb
            CFi_p += cf_p
            CMi_p += cross(Δr, cf_p)

            cf_q = props_q[isurf][i].cfb
            CFi_q += cf_q
            CMi_q += cross(Δr, cf_q)

            cf_r = props_r[isurf][i].cfb
            CFi_r += cf_r
            CMi_r += cross(Δr, cf_r)

            # left bound vortex
            rc = left_center(surfaces[isurf][i])
            Δr = rc - ref.r
            cf = properties[isurf][i].cfl
            CFi += cf
            CMi += cross(Δr, cf)

            cf_a = props_a[isurf][i].cfl
            CFi_a += cf_a
            CMi_a += cross(Δr, cf_a)

            cf_b = props_b[isurf][i].cfl
            CFi_b += cf_b
            CMi_b += cross(Δr, cf_b)

            cf_p = props_p[isurf][i].cfl
            CFi_p += cf_p
            CMi_p += cross(Δr, cf_p)

            cf_q = props_q[isurf][i].cfl
            CFi_q += cf_q
            CMi_q += cross(Δr, cf_q)

            cf_r = props_r[isurf][i].cfl
            CFi_r += cf_r
            CMi_r += cross(Δr, cf_r)

            # right bound vortex
            rc = right_center(surfaces[isurf][i])
            Δr = rc - ref.r
            cf = properties[isurf][i].cfr
            CFi += cf
            CMi += cross(Δr, cf)

            cf_a = props_a[isurf][i].cfr
            CFi_a += cf_a
            CMi_a += cross(Δr, cf_a)

            cf_b = props_b[isurf][i].cfr
            CFi_b += cf_b
            CMi_b += cross(Δr, cf_b)

            cf_p = props_p[isurf][i].cfr
            CFi_p += cf_p
            CMi_p += cross(Δr, cf_p)

            cf_q = props_q[isurf][i].cfr
            CFi_q += cf_q
            CMi_q += cross(Δr, cf_q)

            cf_r = props_r[isurf][i].cfr
            CFi_r += cf_r
            CMi_r += cross(Δr, cf_r)
        end

        # adjust forces from this surface to account for symmetry
        if symmetric[isurf]
            CFi = SVector(2*CFi[1], 0.0, 2*CFi[3])
            CMi = SVector(0.0, 2*CMi[2], 0.0)

            CFi_a = SVector(2*CFi_a[1], 0.0, 2*CFi_a[3])
            CFi_b = SVector(2*CFi_b[1], 0.0, 2*CFi_b[3])
            CFi_p = SVector(2*CFi_p[1], 0.0, 2*CFi_p[3])
            CFi_q = SVector(2*CFi_q[1], 0.0, 2*CFi_q[3])
            CFi_r = SVector(2*CFi_r[1], 0.0, 2*CFi_r[3])

            CMi_a = SVector(0.0, 2*CMi_a[2], 0.0)
            CMi_b = SVector(0.0, 2*CMi_b[2], 0.0)
            CMi_p = SVector(0.0, 2*CMi_p[2], 0.0)
            CMi_q = SVector(0.0, 2*CMi_q[2], 0.0)
            CMi_r = SVector(0.0, 2*CMi_r[2], 0.0)
        end

        # add to body forces
        CF += CFi
        CM += CMi

        CF_a += CFi_a
        CF_b += CFi_b
        CF_p += CFi_p
        CF_q += CFi_q
        CF_r += CFi_r

        CM_a += CMi_a
        CM_b += CMi_b
        CM_p += CMi_p
        CM_q += CMi_q
        CM_r += CMi_r

    end

    # add reference length in moment normalization
    reference_length = SVector(ref.b, ref.c, ref.b)
    CM = CM ./ reference_length

    CM_a = CM_a ./ reference_length
    CM_b = CM_b ./ reference_length
    CM_p = CM_p ./ reference_length
    CM_q = CM_q ./ reference_length
    CM_r = CM_r ./ reference_length

    # pack up derivatives
    dCF = (CF_a, CF_b, CF_p, CF_q, CF_r)
    dCM = (CM_a, CM_b, CM_p, CM_q, CM_r)

    return CF, CM, dCF, dCM
end

"""
    body_forces_history(system, surface_history, property_history; frame=Body())

Return the body force coefficients `CF`, `CM` at each time step in `property_history`.

# Arguments:
 - `system`: Object of type [`System`](@ref) which holds system properties
 - `surface_history`: Vector of surfaces at each time step, where each surface is
    represented by a matrix of surface panels (see [`SurfacePanel`](@ref)) of shape
    (nc, ns) where `nc` is the number of chordwise panels and `ns` is the number
    of spanwise panels
 - `property_history`: Vector of surface properties for each surface at each
    time step, where surface properties are represented by a matrix of panel
    properties (see [`PanelProperties`](@ref)) of shape (nc, ns) where `nc` is
    the number of chordwise panels and `ns` is the number of spanwise panels

# Keyword Arguments
 - `frame`: frame in which to return `CF` and `CM`, options are [`Body()`](@ref) (default),
   [`Stability()`](@ref), and [`Wind()`](@ref)`
"""
function body_forces_history(system, surface_history::AbstractVector{<:AbstractVector{<:AbstractMatrix}},
    property_history; frame=Body())

    # unpack system parameters
    symmetric = system.symmetric
    ref = system.reference[]
    fs = system.freestream[]

    # float type
    TF = eltype(system)

    # number of time steps
    nt = length(property_history)

    # convert single freestream input to vector
    if isa(fs, Freestream)
        fs = fill(fs, nt)
    end

    # initialize time history coefficients
    CF = Vector{SVector{3, TF}}(undef, nt)
    CM = Vector{SVector{3, TF}}(undef, nt)

    # populate time history coefficients
    for it = 1:nt
        CF[it], CM[it] = body_forces(surface_history[it], property_history[it],
            ref, fs[it], symmetric, frame)
    end

    return CF, CM
end

"""
    lifting_line_coefficients(system, r, c; frame=Body())

Return the force and moment coefficients (per unit span) for each spanwise segment
of a lifting line representation of the geometry.

This function requires that a near-field analysis has been performed on `system`
to obtain panel forces.

# Arguments
 - `system`: Object of type [`System`](@ref) that holds precalculated
    system properties.
 - `r`: Vector with length equal to the number of surfaces, with each element
    being a matrix with size (3, ns+1) which contains the x, y, and z coordinates
    of the resulting lifting line coordinates
 - `c`: Vector with length equal to the number of surfaces, with each element
    being a vector of length `ns+1` which contains the chord lengths at each
    lifting line coordinate.

# Keyword Arguments
 - `frame`: frame in which to return `cf` and `cm`, possible options are
    [`Body()`](@ref) (default), [`Stability()`](@ref), and [`Wind()`](@ref)`

# Return Arguments:
 - `cf`: Vector with length equal to the number of surfaces, with each element
    being a matrix with size (3, ns) which contains the x, y, and z direction
    force coefficients (per unit span) for each spanwise segment.
 - `cm`: Vector with length equal to the number of surfaces, with each element
    being a matrix with size (3, ns) which contains the x, y, and z direction
    moment coefficients (per unit span) for each spanwise segment.
"""
function lifting_line_coefficients(system, r, c; frame=Body())
    TF = promote_type(eltype(system), eltype(eltype(r)), eltype(eltype(c)))
    nsurf = length(system.surfaces)
    cf = Vector{Matrix{TF}}(undef, nsurf)
    cm = Vector{Matrix{TF}}(undef, nsurf)
    for isurf = 1:nsurf
        ns = size(system.surfaces[isurf], 2)
        cf[isurf] = Matrix{TF}(undef, 3, ns)
        cm[isurf] = Matrix{TF}(undef, 3, ns)
    end
    return lifting_line_coefficients!(cf, cm, system, r, c; frame)
end

"""
    lifting_line_coefficients!(cf, cm, system, r, c; frame=Body())

In-place version of [`lifting_line_coefficients`](@ref)
"""
function lifting_line_coefficients!(cf, cm, system, r, c; frame=Body())

    # number of surfaces
    nsurf = length(system.surfaces)

    # check that a near field analysis has been performed
    @assert system.near_field_analysis[] "Near field analysis required"

    # extract reference properties
    ref = system.reference[]
    fs = system.freestream[]

    # iterate through each lifting surface
    for isurf = 1:nsurf
        nc, ns = size(system.surfaces[isurf])
        # extract current surface panels and panel properties
        panels = system.surfaces[isurf]
        properties = system.properties[isurf]
        # loop through each chordwise set of panels
        for j = 1:ns
            # calculate segment length
            rls = SVector(r[isurf][1,j], r[isurf][2,j], r[isurf][3,j])
            rrs = SVector(r[isurf][1,j+1], r[isurf][2,j+1], r[isurf][3,j+1])
            ds = norm(rrs - rls)
            # calculate reference location
            rs = (rls + rrs)/2
            # calculate reference chord
            cs = (c[isurf][j] + c[isurf][j+1])/2
            # calculate section force and moment coefficients
            cfj = @SVector zeros(eltype(cf[isurf]), 3)
            cmj = @SVector zeros(eltype(cm[isurf]), 3)
            for i = 1:nc
                # add influence of bound vortex
                rb = top_center(panels[i,j])
                cfb = properties[i,j].cfb
                cfj += cfb
                cmj += cross(rb - rs, cfb)
                # add influence of left vortex leg
                rl = left_center(panels[i,j])
                cfl = properties[i,j].cfl
                cfj += cfl
                cmj += cross(rl - rs, cfl)
                # add influence of right vortex leg
                rr = right_center(panels[i,j])
                cfr = properties[i,j].cfr
                cfj += cfr
                cmj += cross(rr - rs, cfr)
            end
            # update normalization
            cfj *= ref.S/(ds*cs)
            cmj *= ref.S/(ds*cs^2)
            # change coordinate frame
            cfj, cmj = body_to_frame(cfj, cmj, ref, fs, frame)
            # save coefficients
            cf[isurf][:,j] = cfj
            cm[isurf][:,j] = cmj
        end
    end

    return cf, cm
end

"""
    body_to_frame(CF, CM, reference, freestream, frame)

Transform the coefficients `CF` and `CM` from the body frame to the frame
specified in `frame`
"""
body_to_frame

@inline body_to_frame(CF, CM, ref, fs, ::Body) = CF, CM

@inline function body_to_frame(CF, CM, ref, fs, ::Stability)
    R = body_to_stability(fs)
    return R*CF, R*CM
end

@inline function body_to_frame(CF, CM, ref, fs, ::Wind)
    # remove reference lengths
    reflen = SVector(ref.b, ref.c, ref.b)
    CM = CM .* reflen

    # rotate
    R = body_to_wind(fs)
    CF = R*CF
    CM = R*CM

    # reapply reference lengths
    CM = CM ./ reflen

    return CF, CM
end


text:
# In this script we run the model for a month
# for a few different flight interventions.
# In each case we assume the prevalence in the
# external population is 1e-6. Before the intervention,
# we assume 500 daily visitors to Mombassa and
# 1000 daily visitors to Nairobi.

push!(LOAD_PATH, "C:/Users/Joseph Hilton/Documents/GitHub/KenyaCoV/src")
using Plots,Parameters,Distributions
using KenyaCoV

# First do case with no controls:
u0,P,transport_matrix = model_ingredients_from_data("./src/2009_National_Estimates_of_Rural_and_Urban_Populations_by_County.csv","./src/flight_numbers_no_controls.csv","./src/projected_global_prevelance.csv",0.01)
#This method modifies the parameter set for changing the movement structure
KenyaCoV.transportstructure_params!(P,0.001,transport_matrix)
#Then you can modify other parameters
P.τ = 1/1. #e.g. increased treatment rate
#Define initial conditions by modifying the completely susceptible population
u0[30,3,1] += 1#One asymptomatic in Nairobi

#Create a JumpProblem which you can solve --- needs DifferentialEquations module for the solvers
jump_prob_tl = create_KenyaCoV_prob(u0,(0.,360.),P)
#Go straight to solution using solver compiled in the KenyaCoV module
@time sol_tl_no_controls = solve_KenyaCoV_prob(u0,(0.,360.),P,0.25)

# Now do no flights whatsoever:
u0,P,transport_matrix = model_ingredients_from_data("./src/2009_National_Estimates_of_Rural_and_Urban_Populations_by_County.csv","./src/flight_numbers_complete_shutdown.csv","./src/projected_global_prevelance.csv",0.01)
#This method modifies the parameter set for changing the movement structure
KenyaCoV.transportstructure_params!(P,0.001,transport_matrix)
#Then you can modify other parameters
P.τ = 1/1. #e.g. increased treatment rate
#Define initial conditions by modifying the completely susceptible population
u0[30,3,1] += 1#One asymptomatic in Nairobi

#Create a JumpProblem which you can solve --- needs DifferentialEquations module for the solvers
jump_prob_tl = create_KenyaCoV_prob(u0,(0.,360.),P)
#Go straight to solution using solver compiled in the KenyaCoV module
@time sol_tl_complete_shutdown = solve_KenyaCoV_prob(u0,(0.,360.),P,0.25)

# Flights halved after one week:
u0,P,transport_matrix = model_ingredients_from_data("./src/2009_National_Estimates_of_Rural_and_Urban_Populations_by_County.csv","./src/flight_numbers_halved_after_week.csv","./src/projected_global_prevelance.csv",0.01)
#This method modifies the parameter set for changing the movement structure
KenyaCoV.transportstructure_params!(P,0.001,transport_matrix)
#Then you can modify other parameters
P.τ = 1/1. #e.g. increased treatment rate
#Define initial conditions by modifying the completely susceptible population
u0[30,3,1] += 1#One asymptomatic in Nairobi

#Create a JumpProblem which you can solve --- needs DifferentialEquations module for the solvers
jump_prob_tl = create_KenyaCoV_prob(u0,(0.,360.),P)
#Go straight to solution using solver compiled in the KenyaCoV module
@time sol_tl_halved_after_week = solve_KenyaCoV_prob(u0,(0.,360.),P,0.25)

# Flights quartered after one week:
u0,P,transport_matrix = model_ingredients_from_data("./src/2009_National_Estimates_of_Rural_and_Urban_Populations_by_County.csv","./src/flight_numbers_quartered_after_week.csv","./src/projected_global_prevelance.csv",0.01)
#This method modifies the parameter set for changing the movement structure
KenyaCoV.transportstructure_params!(P,0.001,transport_matrix)
#Then you can modify other parameters
P.τ = 1/1. #e.g. increased treatment rate
#Define initial conditions by modifying the completely susceptible population
u0[30,3,1] += 1#One asymptomatic in Nairobi

#Create a JumpProblem which you can solve --- needs DifferentialEquations module for the solvers
jump_prob_tl = create_KenyaCoV_prob(u0,(0.,360.),P)
#Go straight to solution using solver compiled in the KenyaCoV module
@time sol_tl_quartered_after_week = solve_KenyaCoV_prob(u0,(0.,360.),P,0.25)

# Flights stopped completely after one week:
u0,P,transport_matrix = model_ingredients_from_data("./src/2009_National_Estimates_of_Rural_and_Urban_Populations_by_County.csv","./src/flight_numbers_shutdown_after_week.csv","./src/projected_global_prevelance.csv",0.01)
#This method modifies the parameter set for changing the movement structure
KenyaCoV.transportstructure_params!(P,0.001,transport_matrix)
#Then you can modify other parameters
P.τ = 1/1. #e.g. increased treatment rate
#Define initial conditions by modifying the completely susceptible population
u0[30,3,1] += 1#One asymptomatic in Nairobi

#Create a JumpProblem which you can solve --- needs DifferentialEquations module for the solvers
jump_prob_tl = create_KenyaCoV_prob(u0,(0.,360.),P)
#Go straight to solution using solver compiled in the KenyaCoV module
@time sol_tl_shutdown_after_week = solve_KenyaCoV_prob(u0,(0.,360.),P,0.25)


text:
struct Data{Tx<:Real} <: AbstractData
    ## entry::Vector{Float64}

    ## event times and indicators
    T::Vector{Tx}
    E::BitArray{1} ##Vector{Bool}

    ## fixed effects covariates
    Xp::Matrix{Tx} # class-membership model
    Xl::Array{Tx,3} ##Matrix{Tx} # linear mixed model
    Xe::Matrix{Tx} # hazard model

    ## random effects covariates
    Z::Array{Tx,3} ##Matrix{Tx}

    ## repeated measurements
    Y::Matrix{Tx}
    times::Matrix{Tx} # measurement times
    n::Vector{Int} ## number of repeated measurments

    cind::Vector{Int}
    b::Matrix{Tx}

    ## spline basis functions and their integrals evaluated at each event time
    spline::Spline.Spline1D
    haz_basis::Matrix{Tx}
    cumhaz_basis::Matrix{Tx}
end

length(data::Data) = length(data.T)

## required for xval
function slice_data(data::Data, sel)
    Data(data.T[sel], data.E[sel], data.Xp[sel,:], data.Xl[sel,:,:], data.Xe[sel,:], data.Z[sel,:,:], data.Y[sel,:], data.times[sel,:],
         data.n[sel], data.cind[sel], data.b[sel,:], data.spline, data.haz_basis[sel,:], data.cumhaz_basis[sel,:])
end


text:
abstract type Number end
abstract type Real          <: Number end
abstract type FloatingPoint <: Real end
abstract type Integer       <: Real end
abstract type Signed        <: Integer end
abstract type Unsigned      <: Integer end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#####################
# Aggregate flow circuits
# (like a flow circuit but aggregates flows over several examples and batches)
#####################

abstract type AggregateFlowΔNode{O,A} <: DecoratorΔNode{O} end
abstract type AggregateFlowLeafNode{O,A} <: AggregateFlowΔNode{O,A} end
abstract type AggregateFlowInnerNode{O,A} <: AggregateFlowΔNode{O,A} end

struct AggregateFlowLiteral{O,A} <: AggregateFlowLeafNode{O,A}
    origin::O
end

struct AggregateFlowConstant{O,A} <: AggregateFlowLeafNode{O,A}
    origin::O
end

struct AggregateFlow⋀{O,A} <: AggregateFlowInnerNode{O,A}
    origin::O
    children::Vector{<:AggregateFlowΔNode{<:O,A}}
end

mutable struct AggregateFlow⋁{O,A} <: AggregateFlowInnerNode{O,A}
    origin::O
    children::Vector{<:AggregateFlowΔNode{<:O,A}}
    aggr_flow::A
    aggr_flow_children::Vector{A}
end

const AggregateFlowΔ{O,A} = AbstractVector{<:AggregateFlowΔNode{O,A}}

#####################
# traits
#####################

@inline GateType(::Type{<:AggregateFlowLiteral}) = LiteralGate()
@inline GateType(::Type{<:AggregateFlowConstant}) = ConstantGate()

@inline GateType(::Type{<:AggregateFlow⋀}) = ⋀Gate()
@inline GateType(::Type{<:AggregateFlow⋁}) = ⋁Gate()

#####################
# constructors and conversions
#####################

function AggregateFlowΔ(circuit::Δ, ::Type{A}) where {A}
    cache = Dict{ΔNode, AggregateFlowΔNode}()
    sizehint!(cache, length(circuit)*4÷3)

    O = grapheltype(circuit) # type of node in the origin
    af_node(::LiteralGate, n::ΔNode) = AggregateFlowLiteral{O,A}(n)
    af_node(::ConstantGate, n::ΔNode) = AggregateFlowConstant{O,A}(n)

    af_node(::⋀Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        AggregateFlow⋀{O,A}(n, children)
    end

    af_node(::⋁Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        AggregateFlow⋁{O,A}(n, children, zero(A), some_vector(A, num_children(n)))
    end

    map(circuit) do node
        afn = af_node(GateType(node), node)
        cache[node] = afn
        afn
    end
end

#####################
# methods
#####################

@inline literal(n::AggregateFlowLiteral)::Lit  = literal(n.origin)
@inline constant(n::AggregateFlowConstant)::Bool = constant(n.origin)
@inline children(n::AggregateFlowInnerNode) = n.children

function collect_aggr_flows(afc::AggregateFlowΔ, batches::XBatches{Bool})
    reset_aggregate_flows(afc)
    accumulate_aggr_flows(afc, batches)
end

function collect_aggr_flows(fc::FlowΔ, batches::XBatches{Bool})
    @assert origin(fc) isa AggregateFlowΔ
    reset_aggregate_flows(origin(fc))
    accumulate_aggr_flows(fc, batches)
end

# set flow counters to zero
reset_aggregate_flows(afc::AggregateFlowΔ) = foreach(n->reset_aggregate_flow(n), afc)
reset_aggregate_flow(::AggregateFlowΔNode) = () # do nothing
reset_aggregate_flow(n::AggregateFlow⋁{O,A}) where {O,A} = (n.aggr_flow = zero(A) ; n.aggr_flow_children .= zero(A))

const opts_accumulate_flows = (flow_opts★..., compact⋁=false) #keep default options but insist on Bool flows

function accumulate_aggr_flows(afc::AggregateFlowΔ, batches::XBatches{Bool})
    fc = FlowΔ(afc, max_batch_size(batches), Bool, opts_accumulate_flows)
    accumulate_aggr_flows(fc, batches)
end

function accumulate_aggr_flows(fc::FlowΔ, batches::XBatches{Bool})
    for batch in batches
        accumulate_aggr_flows_batch(fc, batch)
    end
end

function accumulate_aggr_flows_batch(fc::FlowΔ, batch::XData{Bool})
    # pass a mini-batch through the flow circuit
    pass_up_down(fc, unweighted_data(batch))
    accumulate_aggr_flows_cached(fc, batch)
end

function accumulate_aggr_flows_cached(fc::FlowΔ, batch::XData{Bool})
    aggregate_data(xd::PlainXData, f::AbstractArray) = sum(f)
    aggregate_data(xd::PlainXData, f::AbstractArray{Bool}) = count(f)
    aggregate_data(xd::WXData, f::AbstractArray) = sum(f .* weights(xd))
    aggregate_data(xd::WXData, f::AbstractArray{Bool}) = sum(weights(xd)[f]) # see AggregateKernelBenchmark

    aggregate_data_factorized(xd::PlainXData, x1::BitVector, xs::BitVector...) = count_conjunction(x1, xs...)
    aggregate_data_factorized(xd::WXData, x1::BitVector, xs::BitVector...) = sum_weighted_product(weights(xd), x1, xs...)

    for n in fc
         # collect flows from mini-batch into aggregate statistics
        if n isa DownFlow⋁
            orig = grand_origin(n)::AggregateFlow⋁
            additional_flow = aggregate_data(batch,downflow(n))
            orig.aggr_flow += additional_flow
            # compute aggregate flows for all but the final child
            for i in 1:num_children(n)-1
                child = n.children[i]
                child_aggr_flow = aggregate_data_factorized(batch, downflow(n), pr_factors(origin(child))...)
                orig.aggr_flow_children[i] += child_aggr_flow
                additional_flow -= child_aggr_flow
            end
            # assing the remainder flow to the final child
            orig.aggr_flow_children[end] += additional_flow
            @assert isapprox(sum(orig.aggr_flow_children), orig.aggr_flow, rtol=0.0001) "Flow is leaking between parent and children: $(origin.aggr_flow) should equal $(origin.aggr_flow_children)"
            # no need to normalize flows here: we can leave this to parameter learning code
        end
    end
end


text:
######################################################################
# random_initializer.jl
######################################################################

"""
`ri = RandomInitializer()`

Random selection of points.
"""
mutable struct RandomInitializer <: Initializer end

function initialize(ri::RandomInitializer, em::ErgodicManager, tm::TrajectoryManager)
	xd = [[tm.x0[1], tm.x0[2]] for i = 1:tm.N+1]
	points = Array(Vector{Float64}, tm.N)
	for i = 1:tm.N
		points[i] = [rand(), rand()]
	end

	tsp_nn!(xd, points)

	ud = compute_controls(xd, tm.h)

	return xd, ud
end


text:
function onedShape(x,r,w)
#= """
  oned_shape.m - computes test functions and derivatives for a Lagrange
                 C0 element given element coordinates and Gauss points.
                 (assumes all nodes are uniformly distributed in the
                 element)

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    x_g, w_g, ϕ, ϕ_x, ϕ_xx = onedShape(x,r,w)
  ```

  Arguments:
  - `x`: Coordinates of the element nodes
  - `r`: Coordinates of Gauss points in (-1,1)
  - `w`: Gauss weights associated with r

  Output arguments:
  - `x_g`: Coordinates of Gauss points in the element
  - `w_g`: Gauss weights scaled by the element Jacobian
  - `ϕ`: Value of element shape functions at x_g
  - `ϕ_x`: First spatial derivatives of ϕ
  - `ϕ_xx`: Second spatial derivatives of ϕ
""" =#

  n  = size(x,1)
  ng = size(r,1)

  xg  = zeros(Float64,ng,n)
  ϕ   = zeros(Float64,ng,n)
  ϕ_x = zeros(Float64,ng,n)
  ϕ_xx = zeros(Float64,ng,n)

  if n == 2
    # Transform coordinates for linear elements
    c0 = ( x[n]-x[1] )/2.0
    c1 = ( x[n]+x[1] )/2.0

    x_g = c0.*r .+ c1

    ϕ[:,1] = ( 1.0 .- r )./2.0
    ϕ[:,2] = ( 1.0 .+ r )./2.0

    ϕ_x[:,2] = (0.5/c0)*ones(Float64,size(r))
    ϕ_x[:,1] = -ϕ_x[:,2]

    djac = c0

    w_g = djac.*w

  elseif (n==3)
    # Transform coordinates for quadratic elements
    c0 = ( x[n]-x[1] )/2.0
    c1 = ( x[n]+x[1] )/2.0

    x_g = c0.*r .+ c1

    ϕ[:,1] = 0.5.*r.*( r.-1.0 )
    ϕ[:,2] = -( r.+1.0 ).*( r.-1.0 )
    ϕ[:,3] = 0.5*r.*( r.+1.0 )

    ϕ_x[:,1] = ( r.-0.5 )/c0
    ϕ_x[:,2] = -2.0*r./c0
    ϕ_x[:,3] = ( r.+0.5 )./c0

    djac = c0

    w_g = djac.*w

    ϕ_xx[:,1] = ones(size(r))/c0^2
    ϕ_xx[:,2] =-2.0*ϕ_xx[:,1]
    ϕ_xx[:,3] = ϕ_xx[:,1]

  elseif (n==4)
    # Transform coordinates for (nonconforming) cubic elements
    c0 = ( x[n]-x[1] )/2.0
    c1 = ( x[n]+x[1] )/2.0

    x_g = c0.*r .+ c1

    r2  = r.*r
    r3  = r.*r2

    ϕ[:,1] = -( 9.0/16.0).*( r3 .- r2      .- r./9.0 .+ 1.0/9.0 )
    ϕ[:,2] =  (27.0/16.0).*( r3 .- r2./3.0 .- r      .+ 1.0/3.0 )
    ϕ[:,3] = -(27.0/16.0).*( r3 .+ r2./3.0 .- r      .- 1.0/3.0 )
    ϕ[:,4] =  ( 9.0/16.0).*( r3 .+ r2      .- r./9.0 .- 1.0/9.0 )

    ϕ_x[:,1] = -( 9.0/16.0).*( 3.0.*r2 .-       2.0.*r .- 1.0/9.0 )./c0
    ϕ_x[:,2] =  (27.0/16.0).*( 3.0.*r2 .- (2.0/3.0).*r .- 1.0     )./c0
    ϕ_x[:,3] = -(27.0/16.0).*( 3.0.*r2 .+ (2.0/3.0).*r .- 1.0     )./c0
    ϕ_x[:,4] =  ( 9.0/16.0).*( 3.0.*r2 .+       2.0.*r .- 1.0/9.0 )./c0

    ϕ_xx[:,1] = -( 9.0/16.0).*( 6.0.*r .- 2.0     )./c0^2
    ϕ_xx[:,2] =  (27.0/16.0).*( 6.0.*r .- 2.0/3.0 )./c0^2
    ϕ_xx[:,3] = -(27.0/16.0).*( 6.0.*r .+ 2.0/3.0 )./c0^2
    ϕ_xx[:,4] =  ( 9.0/16.0).*( 6.0.*r .+ 2.0     )./c0^2

    djac = c0

    w_g = djac.*w
#=
  elseif (n==4)
    # Transform coordinates for (nonconforming) cubic elements
    c0 = ( x[n]-x[1] )/2.0
    c1 = ( x[n]+x[1] )/2.0

    x_g = c0.*r .+ c1

    r2  = r.*r
    r3  = r.*r2

    ϕ[:,1] = -( 9.0/16.0).*( r3 .- r2      .- r./9.0 .+ 1.0/9.0 )
    ϕ[:,2] =  (27.0/16.0).*( r3 .- r2./3.0 .- r      .+ 1.0/3.0 )
    ϕ[:,3] = -(27.0/16.0).*( r3 .+ r2./3.0 .- r      .- 1.0/3.0 )
    ϕ[:,4] =  ( 9.0/16.0).*( r3 .+ r2      .- r./9.0 .- 1.0/9.0 )

    ϕ_r[:,1] = -( 9.0/16.0).*( 3.0.*r2 .-       2.0.*r .- 1.0/9.0 )
    ϕ_r[:,2] =  (27.0/16.0).*( 3.0.*r2 .- (2.0/3.0).*r .- 1.0     )
    ϕ_r[:,3] = -(27.0/16.0).*( 3.0.*r2 .+ (2.0/3.0).*r .- 1.0     )
    ϕ_r[:,4] =  ( 9.0/16.0).*( 3.0.*r2 .+       2.0.*r .- 1.0/9.0 )

    ϕ_rr[:,1] = -( 9.0/16.0).*( 6.0.*r .- 2.0     )
    ϕ_rr[:,2] =  (27.0/16.0).*( 6.0.*r .- 2.0/3.0 )
    ϕ_rr[:,3] = -(27.0/16.0).*( 6.0.*r .+ 2.0/3.0 )
    ϕ_rr[:,4] =  ( 9.0/16.0).*( 6.0.*r .+ 2.0     )

    dxdr = ϕ_r*x(:,1)
    djac = dxdr
    drdx = 1.0./djac

    ϕ_x[:,1] = ϕ_r[:,1].*drdx
    ϕ_x[:,2] = ϕ_r[:,2].*drdx
    ϕ_x[:,3] = ϕ_r[:,3].*drdx
    ϕ_x[:,4] = ϕ_r[:,4].*drdx
    w_g = djac.*w
=#

  else
    error("Elements higher than cubic not currently supported")

  end

  return x_g, w_g, ϕ, ϕ_x, ϕ_xx

end