text:
export mat_back_jacobian, apply_back_jacobian, ng
export test_apply_back, test_mat_back

function ng(f, θ, δ = 1e-5)
    res = []
    for i = 1:length(θ)
        push!(
            res,
            (f(ireplace(θ, i => θ[i] + δ / 2)) - f(ireplace(θ, i => θ[i] - δ / 2))) / δ,
        )
    end
    cat(res..., dims = 3)
end

ireplace(vec::Vector, pair::Pair) = (v = copy(vec); v[pair.first] = pair.second; v)
ireplace(vec::Number, pair::Pair) = pair.second

function mat_back_jacobian(T, block, θ; use_outeradj = false)
    dispatch!(block, θ)
    m = mat(T, block)
    N = size(m, 1)
    jac = zeros(T, size(m)..., length(θ))
    zm = use_outeradj ? OuterProduct(zeros(T, N), zeros(T, N)) : zero(m)
    for j = 1:size(m, 2)
        @inbounds for i = 1:size(m, 1)
            if m[i, j] != 0
                _setval(zm, i, j, 1)
                jac[i, j, :] = mat_back(ComplexF64, block, zm)
                _setval(zm, i, j, 1im)
                jac[i, j, :] += 1im * mat_back(ComplexF64, block, zm)
                _setval(zm, i, j, 0)
            end
        end
    end
    return jac
end
_setval(m::AbstractMatrix, i, j, v) = (m[i, j] = v; m)
_setval(m::OuterProduct, i, j, v) = (m.left[i] = v == 0 ? 0 : 1; m.right[j] = v; m)
Base.setindex!(m::PermMatrix, v, i, j) = m.perm[i] == j ? m.vals[i] = v : error()

function apply_back_jacobian(reg0::AbstractArrayReg{D}, block, θ; kwargs...) where D
    dispatch!(block, θ)
    out = apply!(copy(reg0), block)
    m = out.state
    zm = zero(m)
    jac = zeros(eltype(m), size(m)..., length(θ))
    for j = 1:size(m, 2)
        @inbounds for i = 1:size(m, 1)
            if m[i, j] != 0
                zm[i, j] = 1
                (in, inδ), col =
                    apply_back((copy(out), similar(reg0, copy(zm))), block; kwargs...)
                @assert in ≈ reg0
                jac[i, j, :] .= col
                zm[i, j] *= 1im
                (in, inδ), col = apply_back((copy(out), similar(reg0, copy(zm))), block)
                jac[i, j, :] .+= 1im .* col
                zm[i, j] = 0
            end
        end
    end
    return jac
end

function test_mat_back(
    T,
    block::AbstractBlock{D},
    param;
    δ = 1e-5,
    use_outeradj::Bool = false,
) where {D}
    function mfunc(param)
        dispatch!(block, param)
        mat(T, block)
    end
    # test loss is `real(sum(rand_matrix .* m))`
    got = mat_back_jacobian(T, block, param; use_outeradj = use_outeradj)
    num = ng(mfunc, param, δ)
    res = isapprox(got, num, atol = 10 * δ)
    if !res
        @show size(got)
        @show got
        @show num
    end
    return res
end

function test_apply_back(
    reg0,
    block::AbstractBlock{D},
    param;
    δ = 1e-5,
    kwargs...,
) where {D}
    function mfunc(param)
        block = dispatch(block, param)
        apply!(copy(reg0), block).state
    end
    # test loss is `real(sum(rand_matrix .* m))`
    got = apply_back_jacobian(reg0, block, param; kwargs...)
    num = ng(mfunc, param, δ)
    res = isapprox(got, num, atol = 10 * δ)
    if !res
        @show got
        @show num
    end
    return res
end


text:
"""
$(TYPEDSIGNATURES)

Generate `ODESystem`, dependent variables, and parameters from an `ODEProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
    prob.f isa DiffEqBase.AbstractParameterizedFunction &&
                            return (prob.f.sys, prob.f.sys.states, prob.f.sys.ps)
    @parameters t
    vars = reshape([Variable(:x, i)(t) for i in eachindex(prob.u0)],size(prob.u0))
    params = prob.p isa DiffEqBase.NullParameters ? [] :
             reshape([Variable(:α,i)() for i in eachindex(prob.p)],size(prob.p))
    @derivatives D'~t

    rhs = [D(var) for var in vars]

    if DiffEqBase.isinplace(prob)
        lhs = similar(vars, Any)
        prob.f(lhs, vars, params, t)
    else
        lhs = prob.f(vars, params, t)
    end

    eqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)
    de = ODESystem(eqs,t,vec(vars),vec(params))

    de
end

ext: jl


text:
using Graphs
using Test

testdir = dirname(@__FILE__)

modules = [
           "CDF",
           "Edgelist",
           "GML",
           "NET",
           "DOT",
           "GEXF",
           "Graph6",
           "GraphML",
           "LGCompressed"
          ]

include("graphio.jl")

# write your own tests here
@testset "GraphIO" begin
    for name in modules
        path = joinpath(testdir, name, "runtests.jl")
        include(path)
    end
end


text:
abstract type AbstractDemoCard end

"""
    democard(path::String)::T

Constructs a concrete AbstractDemoCard instance.

The return type `T` is determined by the extension of the path
to your demofile. Currently supported types are:

* [`MarkdownDemoCard`](@ref)
"""
function democard(path::String)::AbstractDemoCard
    validate_file(path)
    _, ext = splitext(path)
    if ext in markdown_exts
        return MarkdownDemoCard(path)
    else
        throw(ArgumentError("unrecognized democard format $(path)"))
    end
end

basename(x::AbstractDemoCard) = basename(x.path)

function get_default_id(card::AbstractDemoCard)
    name_without_ext = splitext(basename(card))[1]
    # default documenter id has -1 suffix
    replace(name_without_ext, ' ' => '-') * "-1"
end

function validate_id(id::AbstractString, card::AbstractDemoCard)
    if occursin(' ', id)
        throw(ArgumentError("invalid id in $(card.path), it should not contain spaces."))
    end
end

function load_config(card::T, key) where T <: AbstractDemoCard
    config = parse(card)

    if key == "cover"
        root = dirname(card.path)
        haskey(config, key) || return nothing

        cover_path = joinpath(root, config[key])
        isfile(cover_path) || throw(ArgumentError("$(cover_path) isn't a valid image file for cover."))
        return cover_path
    elseif key == "id"
        haskey(config, key) || return get_default_id(card)

        id = config[key]
        validate_id(id, card)
        return id
    elseif key == "title"
        return get(config, key) do
            name_without_ext = splitext(basename(card))[1]
            strip(replace(uppercasefirst(name_without_ext), "_" => " "))
        end
    elseif key == "description"
        return get(config, key, card.title)
    else
        throw(ArgumentError("Unrecognized key $(key) for $(T)"))
    end
end

### load concrete implementations

include("markdown.jl")


text:
function nuke_remote(gl::GitLink;
        verbose = true
    )

    ignorestatus = true

    # repo info
    gl_repo = repo_dir(gl)
    url = remote_url(gl)
    cbranch = _curr_branch(gl_repo)

    # set local repo
    _format_repo!(gl_repo, url; verbose)

    # check remote
    rhash0 = _remote_HEAD_hash(url)
    isempty(rhash0) && return false

    # reinit
    _rm(gl_repo)
    mkpath(gl_repo)
    _run("git -C $(gl_repo) init 2>&1"; verbose, ignorestatus)

    # first commit
    _up_push_dummy(gl) # Update push dummy (always push)
    _run("git -C $(gl_repo) add -A 2>&1"; verbose, ignorestatus)
    _run("git -C $(gl_repo) status 2>&1"; verbose, ignorestatus)
    user_name = get_global_config("user.email", "GitLink")
    user_email = get_global_config("user.email", "fake@email.com")
    _run("git -C $(gl_repo) -c user.name='$(user_name)' -c user.email='$(user_email)' commit -am 'Boooom at $(time())' 2>&1"; verbose, ignorestatus)

    # rename branch
    _run("git -C $(gl_repo) branch -m master $(cbranch)"; verbose, ignorestatus)

    # check before
    _run("git --no-pager -C $(gl_repo) log 2>&1"; verbose, ignorestatus)

    # Set remotes
    _run("git -C $(gl_repo) remote add $(_REMOTE_NAME) $(url) 2>&1"; verbose, ignorestatus)
    _run("git -C $(gl_repo) branch --set-upstream $(cbranch) 2>&1"; verbose, ignorestatus)

    # delete remote branch
    _run("git -C $(gl_repo) push -d $(_REMOTE_NAME) $(cbranch) 2>&1"; verbose, ignorestatus)

    # Force push
    _run("git -C $(gl_repo) push $(_REMOTE_NAME) $(cbranch) --force 2>&1"; verbose, ignorestatus)
    _run("git -C $(gl_repo) push $(_REMOTE_NAME) :$(cbranch) 2>&1"; verbose, ignorestatus) # delete obsolete remote branches
    # _run("git -C $(gl_repo) push $(_REMOTE_NAME) --mirror 2>&1"; verbose, ignorestatus)

    # check success
    rhash = _remote_HEAD_hash(url)
    chash = _HEAD_hash(gl_repo)
    if rhash != chash || rhash == rhash0
        _rm(gl_repo) # something fail
        return false
    end
    return true

end


text:
# using Dates

@testset "Spc.jl" begin
    filename = "data\\test.spc"
    dc = DDR2import.Spc.read(filename)
    @test dc["AFI"].name == "AFRICA"
    @test dc["AFI"].type == "AREA"
    @test dc["AFI"].sectors[1].name == "DA"
    @test dc["AFI"].sectors[2].type == "NAS"
    #LAST
    @test dc["BENELUX"].name == "EB/EH"
    @test dc["BENELUX"].sectors[end].name == "EH"
end


text:
"""
Submodule of `SirenSeq`.
Used for playing `SirenSeq` expressions using `csound` or writing them to audio files.
Must have `csound` installed.
"""
module Csound

export

CsHead, SoundFont, FluidInst, ChanMap,
makeOrc, makeScore

using SirenSeq

include("csoundgen.jl")

ChanMap = Dict{Int,CsInstr}

"""
	makeScore(chs::ChanMap, z; path::AbstractString=\"temp\", bpm::Int=defaultBpm)

Writes audio expression `z` to a new `.sco` csound score file;

`chs` mapping channel numbers to instruments \\
`path` is the file name (do not append \".sco\", that is done automatically), \\
`bpm` is the beats per minute, defaults to 120 (default can be changed with `setDefaultBpm`).
"""
function makeScore(chs::ChanMap, z; path::AbstractString="temp", bpm::Int=SirenSeq.defaultBpm)
	print("parsing expression ... \t")
	x = convert(Exp,atomOrExp(z))
	if isa(x,Atom) ; x = C(x) ; end
	SirenSeq.sortExp!(x)
	file = open("$(path).sco","w")
	for a in x.as
		csScoreAtom(file,chs,bpm,a)
	end
	close(file)
	println("score file written to:\n  $(path).sco")
	return
end

## write a single note atom as a csound note
function csScoreAtom(out::IOStream, chs::ChanMap, bpm::Int, a::SirenSeq.Note)
	if !haskey(chs,a.cha) ; throw(ErrorException("no instrument specified for channel $a.cha")) ; end
	n = chs[a.cha]
	ofs = convert(Float64,a.ofs*(240//bpm))
	dur = convert(Float64,a.dur*(240//bpm))
	pch = Scales.getPitch(a.sca,a.ocv,a.itv)
	vel = clamp(convert(Int,round(a.vel*127)),0,127)
	frq = 440*2^((pch-69)/12)
	write(out,"i \"$(instrName(n))\" $(@sprintf("%.6f",ofs)) $(@sprintf("%.6f",dur)) $(convert(Int,floor(frq))) $(vel)\n")
end

## most atoms are ignored as not yet supported
function csScoreAtom(chs::ChanMap, bpm::Int, a::Atom)
end

end


text:
using NeuralOperators
using Documenter
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"), sorting=:nyt)

DocMeta.setdocmeta!(NeuralOperators, :DocTestSetup, :(using NeuralOperators); recursive=true)

include("pages.jl")

makedocs(
    bib,
    modules=[NeuralOperators],
    authors="JingYu Ning <foldfelis@gmail.com> and contributors",
    repo="https://github.com/SciML/NeuralOperators.jl/blob/{commit}{path}#{line}",
    sitename="NeuralOperators.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="http://neuraloperators.sciml.ai",
        assets=String[],
    ),
    pages=pages,
)

deploydocs(;
    repo="github.com/SciML/NeuralOperators.jl",
)


text:
"""
    get_thumbnail(video_id, path)

Save thumbnail image to path. Path must not include extension name (e.g. `"/home/user/filename"`).

| id            | resolution  |
|:--------------|:------------|
| default       | 120x90      |
| mqdefault     | 320x180     |
| hqdefault     | 480x360     |
| sddefault     | 640x480     |
| maxresdefault | 1280x720    |

!!! note
    Some resolutions are different from original one.
    They have black border on Top-Bottom or Left-Right.
"""
function get_thumbnail(v_id, path)
    res = (
        "maxresdefault",
        "sddefault",
        "hqdefault",
        "mqdefault",
        "default"
    )
    for r ∈ res
        u = "https://img.youtube.com/vi/$(v_id)/$(r).jpg"
        try
            HTTP.head(u)
        catch e
            if e.status == 404
                continue
            else
                @warn("Couldn't get thumbnail.", status = e.status)
                break
            end
        end
        HTTP.open(:GET, u) do stream
            open("$path.jpg", "w") do io
                write(io, stream)
            end
        end
        break
    end
end
# http://img.youtube.com
# http://i.ytimg.com
# http://i3.ytimg.com
# WEBP_URL_PATTERN   /vi_webp/OAOP2JUvxYw/maxresdefault.webp

# ====================================================
# ==================    MetaData    ==================
# ====================================================
function get_metadata(v_id)
    response = get_player_response(v_id)
    j = LazyJSON.value(response)
    haskey(j, "streamingData") || return nothing
    haskey(j["streamingData"], "adaptiveFormats") || return nothing
    cs = Dict{Int,NamedTuple}()
    for d ∈ j["streamingData"]["adaptiveFormats"]
        ps = Pair{Symbol,Any}[]
        if haskey(d, "contentLength")
            push!(ps, :length => d["contentLength"])
        end
        if haskey(d, "url")
            push!(ps, :url => d["url"])
        #=elseif haskey(d, "signatureCipher")
            push!(ps, :cipher => d["signatureCipher"])=#
        else
            return nothing
        end
        push!(cs, d["itag"] => (; ps...))
    end
    return (
        title = j["videoDetails"]["title"],
        owner = j["videoDetails"]["author"],
        owner_id = j["videoDetails"]["channelId"],
        description = j["videoDetails"]["shortDescription"],
        length = j["videoDetails"]["lengthSeconds"],
        date = j["microformat"]["playerMicroformatRenderer"]["publishDate"],
        codecs = cs
    )
end

function get_player_response(v_id)
    u = "https://www.youtube.com/get_video_info?video_id=" * v_id
    msg = HTTP.get(u).body |> String
    response = split(msg, '&') |>
        x -> filter(y->startswith(y, "player_response"), x) |>
        first |>
        x -> split(x, '=', limit=2) |>
        last
    #io = IOBuffer(response)
    #out = IOBuffer()
    # searching player_response
    #=while !eof(io)
        c = read(io, Char)
        if c == '&'
            c_next1 = read(io, Char)
            if c_next1 == 'p'
                c_next2 = read(io, Char)
                c_next2 == 'l' && break
            end
        end
    end=#
    # skip to response body
    #=while !eof(io)
        read(io, Char) == '=' && break
    end=#
    unescape_uri(response)
end

# unescape URI encoded character
#   (forked from HTTP.URIs.unescapeuri)
function unescape_uri(str)
    io = IOBuffer(str)
    out = IOBuffer()
    while !eof(io)
        c = read(io, Char)
        c == '&' && break
        if c == '%'
            c1 = read(io, Char)
            c2 = read(io, Char)
            c_decoded = parse(UInt8, string(c1, c2); base=16)
            write(out, c_decoded)
        # spaces have been replaced with '+'
        elseif c == '+'
            write(out, ' ')
        else
            write(out, c)
        end
    end
    String(take!(out))
end

#=function extract_json()
    TODO
end=#

# ====================================================
# ===================    Cipher    ===================
# ====================================================
function iscipher(codecs_dict)
    keys(codecs_dict) |>
        first |>
        x -> codecs_dict[x].url |>
        isempty
end

# decipher(url) = TODO

# ====================================================
# =============    MetaData for FFMPEG    ============
# ====================================================
# convert metadata to FFMPEG format
#
function info2cmd(info, v_id, opts)
    date = format_date(info.date, opts)
    d = Dict(
        "title" => info.title,
        "comment" => info.description,
        "date" => date,
        "artist" => info.owner,
        "publisher" => "https://www.youtube.com/channel/" * info.owner_id,
        "purl" => "https://www.youtube.com/watch?v=" * v_id
        #"encoding_tool" => ""
    )
    cmd_arr = String[]
    for (k,v) ∈ d
        push!(cmd_arr, "-metadata:g")
        push!(cmd_arr, string(k, "=", v))
    end
    return cmd_arr
end

function format_date(date, opts)
    # youtube default delimiter
    default = '-'
    haskey(opts, :date_format) || return replace(date, default=>'/')
    fmt = opts[:date_format]
    if haskey(fmt, :delim)
        replace(date, default=>fmt[:delim])
    elseif haskey(fmt, :custom)
        # TODO
        # user defined format
        #=ud = fmt[:custom]
        year, month, day = split(date, default)
        if occursin("yyyy", ud)
            date = replace(date, "yyyy"=>year)
        elseif occursin("yy", ud)
            date = replace(date, "yy"=>year)
        elseif occursin("y", ud)
            date = replace(date, "yyyy"=>year)
        end=#
        return date
    else
        replace(date, default=>'/')
    end
end


text:
@gen function propose_to_pose(trace, pose_addr, pose, pos_var, rot_conc)
    {pose_addr} ~ gaussianVMF(pose, pos_var, rot_conc)
end
pose_proposal_move(trace, pose_addr, pose, pose_var, rot_conc) = mh(trace, propose_to_pose, (pose_addr, pose, pose_var, rot_conc))

@gen function gaussian_drift_proposal(trace, pose_addr, pos_var, rot_conc)
    {pose_addr} ~ gaussianVMF(trace[pose_addr], pos_var, rot_conc)
end
drift_move(trace, pose_addr, pos_var, rot_conc) = mh(trace, gaussian_drift_proposal, (pose_addr, pos_var, rot_conc))

@gen function gaussian_drift_scalar(trace, pose_addr, var)
    {pose_addr} ~ normal(trace[pose_addr], var)
end
drift_move(trace, pose_addr, var) = mh(trace, gaussian_drift_scalar, (pose_addr, var))

export propose_to_pose, gaussian_drift_proposal, pose_proposal_move, drift_move

# Proposal to rotate an object's pose by 180 degrees about a specified axis.
@gen function pose_flip_proposal(trace, pose_addr, dimension::Int, rot_conc)
    pose = trace[pose_addr]
    euler = zeros(3)
    euler[dimension] = pi
    flipped_pose = get_c_relative_to_a(pose, Pose(R.RotXYZ(euler...)))
    {pose_addr} ~ gaussianVMF(flipped_pose, 0.0001, rot_conc)
end
pose_flip_move(trace, pose_addr, dimension, rot_conc) = mh(trace, pose_flip_proposal, (pose_addr, dimension, rot_conc))

export pose_flip_move

# Pose involution proposal.
@gen function pose_inv_randomness(trace, pose_addr, modifier)
    bit ~ bernoulli(0.5)
    bit
end
function pose_involution(
    trace,
    fwd_choices,
    fwd_ret,
    proposal_args;
    check = false,
)
    pose_addr, modifier = proposal_args
    bit = fwd_ret
    bwd_choices = choicemap(:bit => !bit)

    if bit
        new_pose = get_c_relative_to_a(trace[pose_addr], modifier)
    else
        new_pose = get_c_relative_to_a(trace[pose_addr], inverse_pose(modifier))
    end

    args = get_args(trace)
    argdiffs = map((_) -> NoChange(), args)
    new_trace, weight, _, _ = Gen.update(trace, args, argdiffs, choicemap(pose_addr => new_pose))
    return (new_trace, bwd_choices, weight)
end
function pose_involution_move(trace, pose_addr, modifier; iterations=1)
    acceptances = false
    for _ in 1:iterations
        trace, acc = mh(
            trace,
            pose_inv_randomness,
            (pose_addr, modifier,),
            pose_involution,
            check = false,
        )
        acceptances = acc || acceptances
    end
    trace, acceptances
end

export pose_involution_move

# Mixture of Gaussian VMF proposal

@gen function propose_pose_mixture(trace, addr, poses, weights, posVar, rotConc)
    {addr} ~ gaussianVMF_weighted_mixture(poses, weights, posVar, rotConc)
end
pose_mixture_move(trace, addr, poses, weights, posVar, rotConc) = mh(trace, propose_pose_mixture, (addr, poses, weights,  posVar, rotConc))

export pose_mixture_move


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

ProjDir = @__DIR__

d_str = "dag{k1 -> k2;k2 -> y;v -> x2;w -> k1;x1 -> v;x1 -> w;x2 -> y;x3 -> w;x3 -> y}"

dag = DAG("test_desc_02", d_str)

@testset "m_separation_03" begin

  @test m_separation(dag, :k1, :v, c=[:x1]) == true
  @test m_separation(dag, :k1, :v, c=[:w]) == true
  @test m_separation(dag, :k1, :x1, c=[:w]) == true
  @test m_separation(dag, :k1, :x2, c=[:v]) == true
  @test m_separation(dag, :k1, :x2, c=[:x1]) == true
  @test m_separation(dag, :k1, :x2, c=[:w]) == true
  @test m_separation(dag, :k1, :x3, c=[:w]) == true
  @test m_separation(dag, :k1, :y, c=[:k2, :x2, :x3]) == true
  @test m_separation(dag, :k1, :y, c=[:k2, :v, :x3]) == true
  @test m_separation(dag, :k1, :y, c=[:k2, :x1, :x3]) == true
  @test m_separation(dag, :k1, :y, c=[:k2, :w]) == true
  @test m_separation(dag, :k2, :v, c=[:x1]) == true
  @test m_separation(dag, :k2, :v, c=[:w]) == true
  @test m_separation(dag, :k2, :v, c=[:k1]) == true
  @test m_separation(dag, :k2, :w, c=[:k1]) == true
  @test m_separation(dag, :k2, :x1, c=[:w]) == true
  @test m_separation(dag, :k2, :x1, c=[:k1]) == true
  @test m_separation(dag, :k2, :x2, c=[:v]) == true
  @test m_separation(dag, :k2, :x2, c=[:x1]) == true
  @test m_separation(dag, :k2, :x2, c=[:w]) == true
  @test m_separation(dag, :k2, :x2, c=[:k1]) == true
  @test m_separation(dag, :k2, :x3, c=[:w]) == true
  @test m_separation(dag, :k2, :x3, c=[:k1]) == true
  @test m_separation(dag, :v, :w; c=:x1) == true
  @test m_separation(dag, :v, :x3) == true
  @test m_separation(dag, :v, :y; c=[:k2, :x2, :x3]) == true
  @test m_separation(dag, :v, :y; c=[:k1, :x2, :x3]) == true
  @test m_separation(dag, :v, :y; c=[:w, :x2, :x3]) == true
  @test m_separation(dag, :v, :y; c=[:x2, :x1]) == true
  @test m_separation(dag, :w, :x2; c=[:v]) == true
  @test m_separation(dag, :w, :x2; c=[:x1]) == true
  @test m_separation(dag, :w, :y; c=[:k2, :x2, :x3]) == true
  @test m_separation(dag, :w, :y; c=[:k2, :v, :x3]) == true
  @test m_separation(dag, :w, :y; c=[:k1, :x2, :x3]) == true
  @test m_separation(dag, :w, :y; c=[:k1, :v, :x3]) == true
  @test m_separation(dag, :w, :y; c=[:k1, :x1, :x3]) == true
  @test m_separation(dag, :w, :y; c=[:k2, :x2, :x3]) == true
  @test m_separation(dag, :x1, :x2; c=[:v]) == true
  @test m_separation(dag, :x1, :x3) == true
  @test m_separation(dag, :x1, :y; c=[:k2, :x2, :x3]) == true
  @test m_separation(dag, :x1, :y; c=[:k1, :x2, :x3]) == true
  @test m_separation(dag, :x1, :y; c=[:w, :x2, :x3]) == true
  @test m_separation(dag, :x1, :y; c=[:k2, :v, :x3]) == true
  @test m_separation(dag, :x1, :y; c=[:k1, :v, :x3]) == true
  @test m_separation(dag, :x1, :y; c=[:v, :w, :x3]) == true
  @test m_separation(dag, :x2, :x3) == true

end


text:
function compute_bounds(tp::TargetParameter,
                        bases::Array{Tuple{MTRBasis, MTRBasis}, 1},
                        assumptions::Dict,
                        dgp::DGP)

    ############################################################################
    # Set up problem
    ############################################################################
    m = Model(() -> Clp.Optimizer())
    set_optimizer_attribute(m, "LogLevel", 0)

    L = length(bases)
    J = [length(basis[d + 1].a) for basis in bases, d in 0:1]
    K = [length(basis[d + 1].b) for basis in bases, d in 0:1]
    @variable(m, θ[ℓ = 1:L, d = 0:1, j = 1:J[ℓ, d + 1], k = 1:K[ℓ, d + 1]])

    ############################################################################
    # Observational equivalence for each model
    ############################################################################
    βₛ = compute_βₛ(dgp)
    Γₛ = compute_Γₛ(bases, dgp)
    @constraint(m, obseq[ℓ = 1:L, s = 1:length(βₛ)],
                sum(θ[ℓ, d, j, k] * Γₛ[ℓ, d + 1][s, j, k]
                    for d = 0:1, j = 1:J[ℓ, d + 1], k = 1:K[ℓ, d + 1]) == βₛ[s])

    ############################################################################
    # Target parameters
    ############################################################################
    Γ⭑ = compute_Γ⭑(tp, bases, dgp)
    @objective(m, Min, sum(θ[ℓ, d, j, k] * Γ⭑[ℓ, d + 1][j, k]
                           for ℓ in 1:L, d in 0:1,
                               j in 1:J[ℓ, d + 1], k in 1:K[ℓ, d + 1]))

    ############################################################################
    # Add assumptions
    ############################################################################
    @variable(m, θz[ℓ = 1:L, d = 0:1,
                    k = 1:K[ℓ, d + 1], z = 1:size(dgp.suppZ, 1)])
    @constraint(m, θz_definition[ℓ = 1:L, d = 0:1,
                                 k = 1:K[ℓ, d + 1], z = 1:size(dgp.suppZ, 1)],
                sum(θ[ℓ, d, j, k] * bases[ℓ][d + 1].a[j](dgp.suppZ[z, :])
                    for j in 1:J[ℓ, d + 1]) == θz[ℓ, d, k, z])

    # Boundedness
    if (haskey(assumptions, :lb))
        @constraint(m, lb[ℓ = 1:L, d = 0:1,
                          k = 1:K[ℓ, d + 1], z = 1:size(dgp.suppZ, 1)],
                    θz[ℓ, d, k, z] >= assumptions[:lb])
    end
    if (haskey(assumptions, :ub))
        @constraint(m, ub[ℓ = 1:L, d = 0:1,
                          k = 1:K[ℓ, d + 1], z = 1:size(dgp.suppZ, 1)],
                    θz[ℓ, d, k, z] <= assumptions[:ub])
    end
    if (haskey(assumptions, :lb_diff) | haskey(assumptions, :ub_diff))
        for ℓ in 1:L @assert K[ℓ, 1] == K[ℓ, 2] end
    end
    if (haskey(assumptions, :lb_diff))
        @constraint(m, lb_diff[ℓ = 1:L,
                          k = 1:K[ℓ, 1], z = 1:size(dgp.suppZ, 1)],
                    θz[ℓ, 1, k, z] - θz[ℓ, 0, k, z] >= assumptions[:lb_diff])
    end
    if (haskey(assumptions, :ub_diff))
        @constraint(m, ub_diff[ℓ = 1:L,
                          k = 1:K[ℓ, 1], z = 1:size(dgp.suppZ, 1)],
                    θz[ℓ, 1, k, z] - θz[ℓ, 0, k, z] <= assumptions[:ub_diff])
    end

    # Mutual consistency
    if (haskey(assumptions, :mutually_consistent) &&
            assumptions[:mutually_consistent] && (L >= 2))
        Γₘ  = compute_Γₘ(bases, dgp)
        @constraint(m, mutually_consistent[ℓ = 1:L, d = 0:1, d′ = 0:1,
                                           z = 1:size(dgp.suppZ, 1)],
                    sum(θ[ℓ, d, j, k] * Γₘ[ℓ, d + 1][d′ + 1, z, j, k] for
                        j = 1:J[ℓ, d + 1], k = 1:K[ℓ, d + 1]) ==
                    sum(θ[1, d, j, k] * Γₘ[1, d + 1][d′ + 1, z, j, k] for
                        j = 1:J[1, d + 1], k = 1:K[1, d + 1]))
    end

    # Decreasing level
    if (haskey(assumptions, :decreasing_level) &&
            length(assumptions[:decreasing_level]) > 0)
        for (ℓ, d) in assumptions[:decreasing_level]
            con = @constraint(m, [k = 2:K[ℓ, d + 1], z = 1:size(dgp.suppZ, 1)],
                              θz[ℓ, d, k - 1, z] >= θz[ℓ, d, k, z])
        end
    end

    # Decreasing difference
    if (haskey(assumptions, :decreasing_difference) &&
            length(assumptions[:decreasing_difference]) > 0)
        for ℓ in assumptions[:decreasing_difference]
            @assert (K[ℓ, 1] == K[ℓ, 2])
            con = @constraint(m, [k = 2:K[ℓ, 1], z = 1:size(dgp.suppZ, 1)],
                              θz[ℓ, 1, k - 1, z] - θz[ℓ, 0, k - 1, z] >=
                              θz[ℓ, 1, k, z] - θz[ℓ, 0, k, z])
        end
    end

    ############################################################################
    # Solve problem
    ############################################################################
    function recover_optimal() # this dance is due to JuMP's container handling
        θjump = value.(θ)
        θarray = [[NaN for j in 1:J[ℓ, d + 1], k in 1:K[ℓ, d + 1]]
                  for ℓ in 1:L, d in 0:1]
        for key in keys(θjump.data)
            θarray[key[1], key[2] + 1][key[3], key[4]] = θjump.data[key]
        end
        return [(MTR(bases[ℓ][1], θarray[ℓ, 1]),
                 MTR(bases[ℓ][2], θarray[ℓ, 2])) for ℓ in 1:length(bases)]
    end

    results = Dict{Symbol, Any}(:lb => +Inf, :ub => -Inf)
    optimize!(m)
    if !(termination_status(m) in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED])
        results[:lb] = objective_value(m)
        results[:mtr_lb] = recover_optimal()
        set_objective_sense(m, MOI.MAX_SENSE)
        optimize!(m)
        results[:ub] = objective_value(m)
        results[:mtr_ub] = recover_optimal()
    end

    return results
end


text:
# v0.12 deprecations
@deprecate Dropout(p, dims) Dropout(p; dims=dims)
@deprecate InstanceNorm(λ, β, γ, μ, σ², ϵ, momentum, active=nothing) InstanceNorm(λ, β, γ, μ, σ², ϵ, momentum, true, true, active, length(β))
@deprecate BatchNorm(λ, β, γ, μ, σ², ϵ, momentum, active=nothing) BatchNorm(λ, β, γ, μ, σ², ϵ, momentum, true, true, active, length(β))
@deprecate GroupNorm(G, λ, β, γ, μ, σ², ϵ, momentum, active=nothing) GroupNorm(G, λ, β, γ, μ, σ², ϵ, momentum, true, true, active, length(β))
@deprecate outdims(f, inputsize) outputsize(f, inputsize)
@deprecate Conv(; weight,  bias, activation=identity, kws...) Conv(weight, bias, activation; kws...)
@deprecate ConvTranspose(; weight, bias, activation=identity, kws...) ConvTranspose(weight, bias, activation; kws...)
@deprecate DepthwiseConv(; weight, bias, activation=identity, kws...) DepthwiseConv(weight, bias, activation; kws...)

function Base.getproperty(a::Dense, s::Symbol)
  if s === :W
    Base.depwarn("field name dense.W is deprecated in favour of dense.weight", :Dense)
    return getfield(a, :weight)
  elseif s === :b
    Base.depwarn("field name dense.b is deprecated in favour of dense.bias", :Dense)
    return getfield(a, :bias)
  end
  return getfield(a, s)
end

function ones(dims...)
  Base.depwarn("Flux.ones(size...) is deprecated, please use Flux.ones32(size...) or Base.ones(Float32, size...)", :ones)
  Base.ones(Float32, dims...)
end
ones(T::Type, dims...) = Base.ones(T, dims...)

function zeros(dims...)
  Base.depwarn("Flux.zeros(size...) is deprecated, please use Flux.zeros32(size...) or Base.zeros(Float32, size...)", :zeros)
  Base.zeros(Float32, dims...)
end
zeros(T::Type, dims...) = Base.zeros(T, dims...)

ones32(::Type, dims...) = throw(ArgumentError("Flux.ones32 is always Float32, use Base.ones to specify the element type"))
zeros32(::Type, dims...) = throw(ArgumentError("Flux.zeros32 is always Float32, use Base.zeros to specify the element type"))


text:
"""
    get_highlights_comments(
        target::String;
        concatenate::Bool=true
    ) -> Tuple{Vector{String}, Vector{String}}

Extract the highlights and comments from a passed PDF or all PDFs found recursively
in the passed directory.

# Arguments
- `target::String`: $(TARGET_PDF_ARGUMENT)

# Keywords
- `concatenate::Bool=true`: $(CONCATENATE_KEYWORD)

# Returns
- `Tuple{Vector{String}, Vector{String}}`: the highlights and comments

# Throws
- Exceptions from: [`get_highlights_comments_pages`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf_dir = joinpath(pathof(PDFHighlights) |> dirname |> dirname, "test", "pdf")
path_to_pdf = joinpath(path_to_pdf_dir, "TestPDF.pdf")

get_highlights_comments(path_to_pdf_dir) ==
get_highlights_comments(path_to_pdf) ==
(
    [
        "Highlight 1",
        "Highlight 2 Highlight 3",
        "Highlight 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ],
    ["Comment 1", "Comment 2 Comment 3", "Comment 4", "", "", "", ""],
)

# output

true
```
"""
function get_highlights_comments(
    target::String;
    concatenate::Bool=true,
)::Tuple{Vector{String}, Vector{String}}

    highlights, comments, _ = get_highlights_comments_pages(target; concatenate)
    return highlights, comments

end


text:
# generated by hand ;(

using CEnum

# HACK: a shim
const hipStream_t = Ptr{Cvoid}

const rocblas_int = Int32
const rocblas_long = Int64
const rocblas_float_complex = ComplexF32
const rocblas_double_complex = ComplexF64
const rocblas_half = Int16
const rocblas_half_complex = ComplexF32 # wtf???
const rocblas_handle = Ptr{Nothing}

const ROCBLASFloat = Union{Float32,Float64}
const ROCBLASReal = Union{Float32,Float64}
const ROCBLASComplex = Union{ComplexF32,ComplexF64}

@cenum(rocblas_operation_t,
    ROCBLAS_OPERATION_NONE = 111,
    ROCBLAS_OPERATION_TRANSPOSE = 112,
    ROCBLAS_OPERATION_CONJUGATE_TRANSPOSE = 113,
)

@cenum(rocblas_fill_t,
    ROCBLAS_FILL_UPPER = 121,
    ROCBLAS_FILL_LOWER = 122,
    ROCBLAS_FILL_FULL = 123,
)

@cenum(rocblas_diagonal_t,
    ROCBLAS_DIAGONAL_NON_UNIT = 131,
    ROCBLAS_DIAGONAL_UNIT = 132,
)

@cenum(rocblas_side_t,
    ROCBLAS_SIDE_LEFT = 141,
    ROCBLAS_SIDE_RIGHT = 142,
    ROCBLAS_SIDE_BOTH = 143,
)

@cenum(rocblas_status_t,
    ROCBLAS_STATUS_SUCCESS = 0,
    ROCBLAS_STATUS_INVALID_HANDLE = 1,
    ROCBLAS_STATUS_NOT_IMPLEMENTED = 2,
    ROCBLAS_STATUS_INVALID_POINTER = 3,
    ROCBLAS_STATUS_INVALID_SIZE = 4,
    ROCBLAS_STATUS_MEMORY_ERROR = 5,
    ROCBLAS_STATUS_INTERNAL_ERROR = 6,
)

@cenum(rocblas_datatype_t,
    ROCBLAS_DATATYPE_F16_R = 150, # 16 bit floating point, real
    ROCBLAS_DATATYPE_F32_R = 151, # 32 bit floating point, real
    ROCBLAS_DATATYPE_F64_R = 152, # 64 bit floating point, real
    ROCBLAS_DATATYPE_F16_C = 153, # 16 bit floating point, complex
    ROCBLAS_DATATYPE_F32_C = 154, # 32 bit floating point, complex
    ROCBLAS_DATATYPE_F64_C = 155, # 64 bit floating point, complex
    ROCBLAS_DATATYPE_I8_R = 160, # 8 bit signed integer, real
    ROCBLAS_DATATYPE_U8_R = 161, # 8 bit unsigned integer, real
    ROCBLAS_DATATYPE_I32_R = 162, # 32 bit signed integer, real
    ROCBLAS_DATATYPE_U32_R = 163, # 32 bit unsigned integer, real
    ROCBLAS_DATATYPE_I8_C = 164, # 8 bit signed integer, complex
    ROCBLAS_DATATYPE_U8_C = 165, # 8 bit unsigned integer, complex
    ROCBLAS_DATATYPE_I32_C = 166, # 32 bit signed integer, complex
    ROCBLAS_DATATYPE_U32_C = 167, # 32 bit unsigned integer, complex
)

@cenum(rocblas_pointer_mode_t,
    ROCBLAS_POINTER_MODE_HOST = 0,
    ROCBLAS_POINTER_MODE_DEVICE = 1,
)

@cenum(rocblas_layer_mode_t,
    ROCBLAS_LAYER_MODE_NONE = 0b0000000000,
    ROCBLAS_LAYER_MODE_LOG_TRACE = 0b0000000001,
    ROCBLAS_LAYER_MODE_LOG_BENCH = 0b0000000010,
    ROCBLAS_LAYER_MODE_LOG_PROFILE = 0b0000000100,
)

@cenum(rocblas_gemm_algo,
    ROCBLAS_GEMM_ALGO_STANDARD = 0b0000000000,
)


text:
using SparseArrays
using Random
using SpatialEcology
using Distances
using Test

@testset "ComMatrix" begin
    Random.seed!(1337)

    datb = sprand(Bool, 12,8,0.9)

    dati = sprand(8,13,0.9).*100
    dati = round.(Int, dati)

    datf = sprand(11,9,0.9)

    cmb = ComMatrix(datb)
    @test cmb isa ComMatrix{Bool}
    cmi = ComMatrix(dati,
        speciesnames = ["sp1", "sp2", "sp3", "sp4", "sp5", "sp6", "sp7", "sp8"],
        sitenames = [:a, :b, :c, :d, :e, :f, :g, :h, :i, :j, :k, :l, :m])
    @test cmi isa ComMatrix{<:Union{Int64, Int32}}
    cmf = ComMatrix(datf, sitecolumns = false)

    @test ComMatrix(Matrix(datf)).occurrences == ComMatrix(datf).occurrences

    @test occupancy(cmb) == [8, 7, 7, 8, 8, 6, 6, 7, 8, 6, 6, 8]
    @test occupancy(cmi) == [11, 12, 11, 13, 10, 13, 10, 12]
    @test occupancy(cmf) == [9, 9, 9, 10, 9, 10, 11, 11, 8]

    @test richness(cmb) == [11, 12, 11, 9, 11, 10, 11, 10]
    @test richness(cmi) == [7, 7, 7, 6, 8, 8, 6, 7, 8, 8, 7, 6, 7]
    @test richness(cmf) == [8, 7, 8, 7, 7, 8, 8, 8, 7, 9, 9]

    @test nsites(cmb) == 8
    @test nsites(cmi) == 13
    @test nsites(cmf) == 11

    @test nspecies(cmb) == 12
    @test nspecies(cmi) == 8
    @test nspecies(cmf) == 9

    @test occurring(cmb) == 1:12
    @test occurring(cmi) == 1:8
    @test occurring(cmf) == 1:9

    @test occupied(cmb) == 1:8
    @test occupied(cmi) == 1:13
    @test occupied(cmf) == 1:11

    @test noccurring(cmb) == 12
    @test noccurring(cmi) == 8
    @test noccurring(cmf) == 9

    @test noccupied(cmb) == 8
    @test noccupied(cmi) == 13
    @test noccupied(cmf) == 11

    @test getspecies(cmb, 3) == [false, true, true, true, true, true, true, true]
    @test getspecies(cmb, "species2") == [true, true, true, true, false, true, true, true]
    @test getspecies(cmi, 2) == [15, 15, 0, 2, 28, 99, 25, 29, 23, 90, 96, 7, 43]
    gcm = getspecies(cmf, 6)
    @test gcm[3] == 0.08724510824225518
    @test gcm isa SubArray

    @test speciesnames(cmi)[8] == "sp8"
    @test speciesnames(cmf)[3] == "species3"

    @test sitenames(cmi)[4] == "d"
    @test sitenames(cmb)[5] == "site5"

    @test length(sitenames(cmf)) == nsites(cmf)
    @test length(sitenames(cmi)) == nsites(cmi)
    @test length(sitenames(cmb)) == nsites(cmb)

    @test length(speciesnames(cmf)) == nspecies(cmf)
    @test length(speciesnames(cmi)) == nspecies(cmi)
    @test length(speciesnames(cmb)) == nspecies(cmb)

    @test sitetotals(cmb) == [11, 12, 11, 9, 11, 10, 11, 10]
    @test sitetotals(cmb) == richness(cmb)
    @test sitetotals(cmi) == [351, 228, 359, 291, 357, 496, 319, 242, 297, 356, 443, 316, 357]
    @test length(sitetotals(cmf)) == nsites(cmf)
    @test sitetotals(cmf)[5] ≈ 2.918329485324914

    @test speciestotals(cmb) == occupancy(cmb)
    @test speciestotals(cmi) == [511, 472, 599, 620, 423, 664, 464, 659]
    @test length(speciestotals(cmf)) == nspecies(cmf)
    @test speciestotals(cmf)[8] ≈ 4.780138474567269

    @test size(cmb) == (12, 8)
    @test size(cmi, 1) == 8

    @test cooccurring(cmb, 1, 3) == [false, true, true, true, true, true, true, true]
    @test cooccurring(cmf, [8, 3]) == [true, true, true, false, true, false, true, true, true, true, true]

    dist = pairwise(BrayCurtis(), view(cmi, sites=1:2))
    @test round.(dist, digits=1) == [0.0 0.3; 0.3 0.0]
  # getindex and setindex are to do

end


text:
"""
    Physics(C0, T0, a, B, Cd, κ, r)

Create data structure for model physics.

See also: [`ClimateModel`](@ref)
"""
mutable struct Physics
    c0::Float64  # initial CO2e concentration, relative to PI
    T0::Float64  # initial temperature, relative to PI

    a::Float64   # logarithmic CO2 forcing coefficient
    B::Float64   # feedback parameter

    Cd::Float64  # deep ocean heat capacity
    κ::Float64   # deep ocean heat uptake rate

    r::Float64   # long-term airborne fraction of CO2e
end


text:
module NURBS

include("Utils.jl")
include("Bases.jl")
include("Objects.jl")

end  # module NURBS


text:
"""
```
measurement{T<:AbstractFloat}(m::Model904{T}, TTT::Matrix{T}, RRR::Matrix{T},
                              CCC::Vector{T})
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::Model904{T},
                                       TTT::Matrix{T},
                                       RRR::Matrix{T},
                                       CCC::Vector{T}) where {T<:AbstractFloat}
    endo     = m.endogenous_states
    endo_new = m.endogenous_states_augmented
    exo      = m.exogenous_shocks
    obs      = m.observables

    _n_observables = n_observables(m)
    _n_states = n_states_augmented(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    ## Output growth - Quarterly!
    ZZ[obs[:obs_gdp], endo[:y_t]]      = 1.0
    ZZ[obs[:obs_gdp], endo_new[:y_t1]] = -1.0
    ZZ[obs[:obs_gdp], endo[:z_t]]      = 1.0
    DD[obs[:obs_gdp]]                  = 100*(exp(m[:z_star])-1)

    ## Hours growth
    ZZ[obs[:obs_hours], endo[:L_t]] = 1.0
    DD[obs[:obs_hours]]             = m[:Lmean]

    ## Labor Share/real wage growth
    ZZ[obs[:obs_wages], endo[:w_t]]      = 1.0
    ZZ[obs[:obs_wages], endo_new[:w_t1]] = -1.0
    ZZ[obs[:obs_wages], endo[:z_t]]      = 1.0
    DD[obs[:obs_wages]]                  = 100*(exp(m[:z_star])-1)

    ## Inflation (GDP Deflator)
    ZZ[obs[:obs_gdpdeflator], endo[:π_t]]  = 1.0
    DD[obs[:obs_gdpdeflator]]              = 100*(m[:π_star]-1)

    ## Nominal interest rate
    ZZ[obs[:obs_nominalrate], endo[:R_t]] = 1.0
    DD[obs[:obs_nominalrate]]             = m[:Rstarn]

    ## Consumption Growth
    ZZ[obs[:obs_consumption], endo[:c_t]]      = 1.0
    ZZ[obs[:obs_consumption], endo_new[:c_t1]] = -1.0
    ZZ[obs[:obs_consumption], endo[:z_t]]      = 1.0
    DD[obs[:obs_consumption]]                  = 100*(exp(m[:z_star])-1)

    ## Investment Growth
    ZZ[obs[:obs_investment], endo[:i_t]]      = 1.0
    ZZ[obs[:obs_investment], endo_new[:i_t1]] = -1.0
    ZZ[obs[:obs_investment], endo[:z_t]]      = 1.0
    DD[obs[:obs_investment]]                  = 100*(exp(m[:z_star])-1)

    ## Spreads
    ZZ[obs[:obs_spread], endo[:ERtil_k_t]] = 1.0
    ZZ[obs[:obs_spread], endo[:R_t]]       = -1.0
    DD[obs[:obs_spread]]                   = 100*log(m[:spr])

    ## 10 yrs infl exp
    TTT10                          = (1/40)*((Matrix{Float64}(I, size(TTT, 1), size(TTT, 1)) - TTT)\(TTT - TTT^41))
    ZZ[obs[:obs_longinflation], :] = TTT10[endo[:π_t], :]
    DD[obs[:obs_longinflation]]    = 100*(m[:π_star]-1)

    QQ[exo[:g_sh], exo[:g_sh]]           = m[:σ_g]^2
    QQ[exo[:b_sh], exo[:b_sh]]           = m[:σ_b]^2
    QQ[exo[:μ_sh], exo[:μ_sh]]           = m[:σ_μ]^2
    QQ[exo[:z_sh], exo[:z_sh]]           = m[:σ_z]^2
    QQ[exo[:λ_f_sh], exo[:λ_f_sh]]       = m[:σ_λ_f]^2
    QQ[exo[:λ_w_sh], exo[:λ_w_sh]]       = m[:σ_λ_w]^2
    QQ[exo[:rm_sh], exo[:rm_sh]]         = m[:σ_r_m]^2
    QQ[exo[:σ_ω_sh], exo[:σ_ω_sh]]       = m[:σ_σ_ω]^2
    QQ[exo[:μ_e_sh], exo[:μ_e_sh]]       = m[:σ_μ_e]^2
    QQ[exo[:γ_sh], exo[:γ_sh]]           = m[:σ_γ]^2
    QQ[exo[:π_star_sh], exo[:π_star_sh]] = m[:σ_π_star]^2

    # These lines set the standard deviations for the anticipated shocks. They
    # are calibrated using the standard deviation of the unanticipated shocks.
    for i = 1:n_anticipated_shocks(m)
        ZZ[obs[Symbol("obs_nominalrate$i")], :]              = ZZ[obs[:obs_nominalrate], :]'*(TTT^i)
        DD[obs[Symbol("obs_nominalrate$i")]]                 = m[:Rstarn]

        if subspec(m) in ["ss9", "realtime1"]
            # Use σ_r_m as SD for all anticipated shocks without normalizing
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[:σ_r_m]^2

        elseif subspec(m) in ["ss10", "ss11", "ss13", "realtime2"]
            # Normalize s.t. sum of variances of anticipated shocks equals
            # the variance of the contemporaneous shock
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[:σ_r_m]^2 / n_anticipated_shocks(m)

        elseif subspec(m) == "realtime3"
            # Normalize s.t. sum of variances of anticipated shocks equals 1/2 the
            # variance of the contemporaneous shock
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[:σ_r_m]^2 / (2*n_anticipated_shocks(m))

        else
            error("Subspec $(subspec(m)) is not defined.")
        end
    end

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0)
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
"""
jInv.ForwardShare

This module contains generic methods for forward problem computations. It handles the
parallelisms and abstractions used in jInv to allow fast and flexible solution of
PDE parameter estimation problems.

Most users will not use ForwardShare directly, but extend its methods by customized
implementations of their problems. Doing so, users can benefit from the full parallel
functionality.

Here is a quick step-by-step introuction on how to define your own forward problem in
a way that is compatible with jInv.

1) Create your own type (documented param that describes a forward problem) as a subtype
   of the generic ForwardProbType

   		type YourParam <: ForwardProbType
      		Mesh
      		Sources
      		Receivers
          Fields
	  	end

2) Import and extend the method 'getData' for your type. Implement it such that it solves
   one forward problem that is stored locally on a worker. The I/O should be as follows

	    function jInv.ForwardShare.getData(sigma,param::YourParam)
        # solve your problem
	    	return dobs, param
	    end

3) Import and extend methods for computing matrix vector products with the sensitivity
   matrix for your problem.

      function jInv.ForwardShare.getSensMatVec(x,sigma,param::YourParam)
        # compute stuff
   			return Sens*x
      end

      function jInv.ForwardShare.getSensTMatVec(x,sigma,param::YourParam)
        # compute stuff
      	return transpose(Sens) * x
      end

	4) clear! method for your type to get rid of unnecessary temp results (e.g., fields, factorizations, preconditioners )

  5) With the previous 4 steps you are already in good shape using ForwardShare to
     reconstruct parameters in parallel or in combination with other physical models.
     For some advanced functionality you will need to extend the methods

        jInv.ForwardShare.getNumberOfData
        jInv.ForwardShare.getSensMatSize
"""
module ForwardShare

	export ForwardProbType
	abstract type ForwardProbType end

	using jInv.Mesh
	using jInv.Utils

	using Distributed
	using SparseArrays
	using LinearAlgebra

	export getSensMatVec
	"""
	Jv  = getSensMatVec(v::Vector,m::Vector,param::ForwardProbType)

	Computes matrix-vector product with the Jacobian.

	"""
	getSensMatVec(v::Vector,m::Vector,param::ForwardProbType) = error("nyi")
	export getSensTMatVec
	"""
	JTv  = getSensMatVec(v::Vector,m::Vector,param::ForwardProbType)

	Computes matrix-vector product with the transpose of Jacobian. Implementation
	depends on forward problem.

	"""
	getSensTMatVec(v::Vector,m::Vector,param::ForwardProbType) = error("nyi")

	"""
	(m,n) = getSensMatSize(pFor)

	Returns size of sensitivity matrix where m is the number of data points
	and n the number of parameters in the model.

	Input

		pFor - forward problem:: Union{ForwardProbType, Array, RemoteChannel}

	This is problem dependent and should be implemented in the respective
	packages.
	"""
	getSensMatSize(pFor::ForwardProbType) = error("getSensMatSize not implemented for forward problems of type $(typeof(pFor))")

	function getSensMatSize(pFor::RemoteChannel)
		if pFor.where != myid()
			return remotecall_fetch(getSensMatSize,pFor.where,pFor)
		else
			return getSensMatSize(fetch(pFor))
		end
	end

	function getSensMatSize(pFor::Array)
		n  = length(pFor)
		sz = [0;0]
		for k=1:n
			(s1,s2) = getSensMatSize(pFor[k])
			sz[2]=s2
			sz[1]+=s1
		end
		return tuple(sz...)
	end

  """
  nd = getNumberOfData(pFor)

  Returns number of data in forward problem
  """
  getNumberOfData(pFor::ForwardProbType) = error("getNumberOfData not implemented for forward problems of type $(typeof(pFor))")

  function getNumberOfData(pFor::RemoteChannel)
		if pFor.where != myid()
			return remotecall_fetch(getNumberOfData,pFor.where,pFor)
		else
			return getNumberOfData(fetch(pFor))
		end
	end

	function getNumberOfData(pFor::Array)
		n  = length(pFor)
		nd = 0
		for k=1:n
			nd += getNumberOfData(pFor[k])
		end
		return nd
	end

	export getSensTMatVec,getSensMatVec, getSensMatSize, getNumberOfData

	Mesh2MeshTypes = Union{RemoteChannel, Future, SparseMatrixCSC, AbstractFloat}

	# # ===== Methods for parallelization =====
	include("getDataParallel.jl")
	include("prepareMesh2Mesh.jl")
	include("interpLocalToGlobal.jl")
	include("getSensMat.jl")
	include("testing.jl")

	import jInv.Utils.clear!
	function clear!(P::ForwardProbType;clearAinv::Bool=true,clearFields::Bool=true, clearMesh::Bool=false, clearSources::Bool=false, clearObs::Bool=false,clearAll::Bool=false)
		if clearAll || clearMesh
			P.M = clear!(P.M)
		end
		if clearAll || clearSources
			P.Sources = clear!(P.Sources)
		end
		if clearAll || clearObs
			P.Obs     = clear!(P.Obs)
		end
		if clearAll || clearFields
		if any(fieldnames(typeof(P)) .== :Fields)
			P.Fields = clear!(P.Fields)
		end

	end

	if clearAll || clearAinv
	if any(fieldnames(typeof(P)) .== :Ainv)
	   clear!(P.Ainv)
	end
	end

end

end


text:
using Test
using WaveProp.Geometry

@testset "Line" begin
    l = ReferenceLine()
    @test ambient_dimension(l) == 1
    @test geometric_dimension(l) == 1
    x = SVector(0.5)
    @test x ∈ l
    x = SVector(1.0)
    @test x ∈ l
    x = SVector(1.1)
    @test !in(x,l)
end
@testset "Triangle" begin
    t = ReferenceTriangle()
    @test ambient_dimension(t) == 2
    @test geometric_dimension(t) == 2
    x = SVector(0.5,0.5)
    @test x ∈ t
    x = SVector(1.0,0.0)
    @test x ∈ t
    x = SVector(1.1,0.0)
    @test !in(x,t)
end
@testset "Tetrahedron" begin
    t = ReferenceTetrahedron()
    @test ambient_dimension(t) == 3
    @test geometric_dimension(t) == 3
    x = SVector(0.5,0.5,0.0)
    @test x ∈ t
    x = SVector(1.0,0.0,0.0) # point on edge
    @test x ∈ t
    x = SVector(1.1,0.0,0.0)
    @test !in(x,t)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
#Code written by Skipper af Rosenborg and Esben Bager
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap

include("ParallelModelGeneration.jl")
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#LOAD DATA
#Skipper's path
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data"

#HPC path
#inputArg = parse(Int64, ARGS[1]) #Should range from 0 to 106
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/"
#mainData = loadCPUData(path)
#mainData = loadConcrete(path)
#mainData = loadHousingData(path)
#mainData = loadDiabetes(path)
mainData = loadWineQuality(path)
mainDataArr = Array{Float64}(mainData)
#mainData = loadIndexDataNoDurLOGReturn(path)

## PIPELINE
dataFile="WineQualityRed"

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results"
path = "/zhome/9f/d/88706/SpecialeCode/Results"
#### MUST CHANGE ####
fileName = path*"/"*dataFile*"/"

colNames = names(mainData)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

function solveStage1(trainingData)
	### STAGE 1 ###
	println("STAGE 1 INITIATED")
	#Define solve
	m = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads = 1))

	#Add binary variables variables
	@variable(m, 0 <= z[1:nCols-1] <= 1, Bin )

	#Calculate highly correlated matix
	HC = cor(trainingData)

	#Define objective function
	@objective(m, Max, sum(z[i] for i=1:length(z)))

	#Define max correlation and constraints
	rho = 0.8
	for i=1:length(z)
		for j=1:length(z)
			if i != j
				if HC[i,j] >= rho
					@constraint(m,z[i]+z[j] <= 1)
				end
			end
		end
	end

	#Get solution status
	status = solve(m);

	#Get objective value
	println("Solved stage 1 with kMax = ", getobjectivevalue(m))
	return kmax = getobjectivevalue(m)
end

function buildAndSolveStage2(standX, standY, curKmax, gamma, warmstartBool, warmStartSol, cuts, fixedZ, HC)
	bigM = 5

	HCPairCounter = 0

	maxTime = 90
	#Define parameters and model
	stage2Model = JuMP.Model(solver = GurobiSolver(TimeLimit = maxTime, OutputFlag = 0, Threads = 1, PreMIQCPForm=0, MIPFocus=1, ImproveStartTime=40))

	#Define variables
	@variable(stage2Model, b[1:bCols]) #Beta values
	@variable(stage2Model, 0 <= z[1:bCols] <= 1, Bin ) #binary z variable
	@variable(stage2Model, v[1:bCols]) # auxiliary variables for abs
	@variable(stage2Model, T) #First objective term
	@variable(stage2Model, G) #Second objective term

	if warmstartBool == true
		for j in find(warmStartSol)
			setvalue(b[j], warmStartSol[j])
			setvalue(v[j], abs(warmStartSol[j]))
			setvalue(z[j], 1)
		end
	end

	#Define objective function (5a)
	@objective(stage2Model, Min, T+G)
	@constraint(stage2Model, soc, norm( [1-T;2*(standX*b-standY)] ) <= 1+T)
	@constraint(stage2Model, gammaConstr, gamma*ones(bCols)'*v <= G) #4bCols

	#Define constraints (5c)
	@constraint(stage2Model, conBigMN, -1*b .<= bigM*z) #from 0 to bCols -1
	@constraint(stage2Model, conBigMP,  1*b .<= bigM*z) #from bCols to 2bCols-1

	#Second objective term
	@constraint(stage2Model, 1*b .<= v) #from 2bCols to 3bCols-1
	@constraint(stage2Model, -1*b .<= v) #from 3bCols to 4bCols-1

	#Define kmax constraint (5d)
	@constraint(stage2Model, kMaxConstr, sum(z[j] for j=1:bCols) <= curKmax) #4bCols+1

	#Constraint 5f - can only select one of a pair of highly correlated features
	rho = 0.8
	for k=1:bCols
		for j=1:bCols
			if k != j
				if HC[k,j] >= rho
					HCPairCounter += 1
					@constraint(stage2Model,z[k]+z[j] <= 1) #from 4bCols+1 to (4bCols+1+HCPairCounter)
				end
			end
		end
	end

	#Constraint (5g) - only one transformation allowed (x, x^2, log(x) or sqrt(x))
	for j=1:(nCols-1)
		@constraint(stage2Model, z[j]+z[j+(nCols-1)]+z[j+2*(nCols-1)]+z[j+3*(nCols-1)] <= 1) #from (4bCols+1+HCPairCounter) to (4bCols+1+HCPairCounter+nCols)
	end

	#Implement cuts
	if !isempty(cuts)
		#nEmptyCuts = 19
		#zOnes = ones(bCols)
		for c = 1:size(cuts)[1]
			@constraint(stage2Model, sum(z[i] for i=find(cuts[c,1:end])) <= countnz(cuts[c,1:end])-1)
			#@constraint(stage2Model, sum(zOnes[i]*z[i] for i=1:bCols) <= bCols) #(4bCols+1+HCPairCounter+nCols) to (4bCols+1+HCPairCounter+nCols+18)
		end
	end

	if !isempty(fixedZ)
		@constraint(stage2Model, sum(z[i] for i=find(fixedZ)) == countnz(fixedZ))
	end

	status = solve(stage2Model)

	if status == :InfeasibleOrUnbounded
		bSolved = zeros(bCols)
		obj = 9999
		objBound = 9999
	else
		bSolved = getvalue(b)
		obj = getobjectivevalue(stage2Model)
		objBound = getobjectivebound(stage2Model)
	end
	return bSolved, obj, objBound
end

function stageThree(best3Beta, X, Y, allCuts)
	bCols = size(X)[2]
	nRows = size(X)[1]
	cuts = Matrix(0, bCols+1)
	rowsPerSample = nRows #All of rows in training data to generate beta estimates, but selected with replacement
	totalSamples = 50 #25 different times we will get a beta estimate
	nBoot = 1000

	#For loop start
	for i = 1:size(best3Beta)[1]
		if signifBoolean[i] == 1 #if the previous solution was already completely significant, skip the bootstrapping
			continue
		end
		bestBeta = best3Beta[i,4:end]
		bestK = best3Beta[i,2]
		bestGamma = best3Beta[i,3]

		xColumns = []
		bSample = Matrix(totalSamples, bCols)

		bZeros = zeros(bCols)
		for j = 1:bCols
			if !isequal(bestBeta[j],0)
				push!(xColumns, j)
			end
		end

		selectedX = X[:,xColumns]
		condNumber = cond(selectedX)
		if condNumber >= 15
			bZeros[xColumns] = 1
			subsetSize = size(xColumns)[1]
			newCut = [bZeros' subsetSize]
			cuts = [cuts; newCut]
			println("A cut based on Condition number = $condNumber has been created from Beta$i")
		end

		# test significance
		bestZ = zeros(bestBeta)
		for l=1:size(bestBeta)[1]
			if bestBeta[l] != 0
				bestZ[l] = 1
			end
		end

		bZeros = zeros(bCols)
		bSample = createBetaDistribution(bSample, X, Y, bestK, totalSamples, rowsPerSample,  bestGamma, allCuts, bestZ) #standX, standY, k, sampleSize, rowsPerSample

		confArray99 = createConfidenceIntervalArray(bSample, nBoot, 0.99)
		confArray95 = createConfidenceIntervalArray(bSample, nBoot, 0.95)
		confArray90 = createConfidenceIntervalArray(bSample, nBoot, 0.90)

		significanceResult = testSignificance(confArray99, confArray95, confArray90, bestBeta)
		significanceResultNONSI = [] # = significanceResult[xColumns]
		subsetSize = size(xColumns)[1]
		for n = 1:size(significanceResult)[1]
			for s = 1:subsetSize
				if significanceResult[n] == 0 && xColumns[s] == n
					push!(significanceResultNONSI,xColumns[s])
					println("Parameter $n is selected, but NOT significant")
				elseif significanceResult[n] > 0 && xColumns[s] == n
					println("Parameter $n is significant with ", significanceResult[n])
				end
			end
		end

		if !isempty(significanceResultNONSI)
			bZeros[significanceResultNONSI] = 1
			subsetSize = size(significanceResultNONSI)[1]
			newCut = [bZeros' subsetSize]
			cuts = [cuts; newCut]
			println("A cut based on parameters being non-significant in Beta$i has been created")
		end
		if isempty(significanceResultNONSI)
			signifBoolean[i] = 1
		end
	end
	return cuts
end

function checkRSquared(standX, standY, bSolved, best3Beta, k, gamma)
	RSquared = getRSquared(standX,standY,bSolved)

	if RSquared > best3Beta[1,1] #largest than 3rd largest RSquared
		if RSquared > best3Beta[2,1] #largest than 2nd largest RSquared
			if RSquared > best3Beta[3,1] #largest than largest RSquared
				best3Beta[1,:] = best3Beta[2,:]
				best3Beta[2,:] = best3Beta[3,:]
				best3Beta[3,1] = RSquared
				best3Beta[3,2] = k
				best3Beta[3,3] = gamma #potentially store actual gamma value
				best3Beta[3,4:end] = bSolved
			else
				best3Beta[1,:] = best3Beta[2,:]
				best3Beta[2,1] = RSquared
				best3Beta[2,2] = k
				best3Beta[2,3] = gamma #potentially store actual gamma value
				best3Beta[2,4:end] = bSolved
			end
		else
			best3Beta[1,1] = RSquared
			best3Beta[1,2] = k
			best3Beta[1,3] = gamma #potentially store actual gamma value
			best3Beta[1,4:end] = bSolved
		end
	end

	return best3Beta
end

bCols = 0
nRows = 0
signifBoolean = 0
best3Beta = 0
run = parse(Int64, ARGS[1])

## SPLIT AND TRANSFORM DATA
#for run=1:10
halfRows = Int64(floor(size(mainDataArr)[1]/2))
train, test = splitDataIn2(mainDataArr, halfRows, size(mainData)[1])
halfRows = Int64(floor(size(test)[1]/2))
test, valid = splitDataIn2(mainDataArr, halfRows, size(test)[1])

# Standardize

standX = zScoreByColumn(expandWithTransformations(train[:,1:end-1]))
standY = zscore(train[:,end])
standXTest = zScoreByColumn(expandWithTransformations(test[:,1:end-1]))
standYTest = zscore(test[:,end])
standXVali = zScoreByColumn(expandWithTransformations(valid[:,1:end-1]))
standYVali = zscore(valid[:,end])

#=
standX = zScoreByColumn(train[:,1:end-1])
standY = zscore(train[:,end])
standXTest = zScoreByColumn(test[:,1:end-1])
standYTest = zscore(test[:,end])
standXVali = zScoreByColumn(valid[:,1:end-1])
standYVali = zscore(valid[:,end])
=#

# Output for Lasso regression in R
writedlm(fileName*string(run)*"XTrain.CSV",standX,",")
writedlm(fileName*string(run)*"YTrain.CSV",standY,",")
writedlm(fileName*string(run)*"XTest.CSV",standXTest,",")
writedlm(fileName*string(run)*"YTest.CSV",standYTest,",")
writedlm(fileName*string(run)*"XVali.CSV",standXVali,",")
writedlm(fileName*string(run)*"YVali.CSV",standYVali,",")

allDataTrain = hcat(standX, standY)
allDataTest = hcat(standXTest, standYTest)
allDataVali = hcat(standXVali, standYVali)
bCols = size(standX)[2]
nRows = size(standX)[1]

trainingData = mainDataArr[:,1:nCols-1]

kmax = solveStage1(trainingData)
#kmax = 40
#Get solution value
#zSolved = getvalue(z)
println("STAGE 1 DONE")

### STAGE 2 ###
println("STAGE 2 INITIATED")
bigM = 5
amountOfGammas = 5

kValue = []
RSquaredValue = []
bSolved = []
bestBeta = []
prevSolutions = zeros(Int64(kmax*amountOfGammas), bCols)

stage2Model = 0
best3Beta = zeros(3,bCols+3)
best3Beta[:,1] = -1e3
signifBoolean = zeros(3)

maxVal = findmax(abs.(standX'*standY))[1]

## RUN LASSO
kZeroCounter = 0
lassoOutput = zeros(300,5+bCols)
tic()
for i = 1:300
	gammaArray = linspace(maxVal, 0, 300)
	lassoOutput[i,5] = gammaArray[i]
	lassoOutput[i,1], lassoOutput[i,2], lassoOutput[i,3], lassoOutput[i,4], lassoOutput[i,6:end] = lassoValidationAndTest(gammaArray[i], standX, standY, standXVali, standYVali, standXTest, standYTest)
	if lassoOutput[i,4] == 0
		kZeroCounter += 1
	end
	if kZeroCounter > 1
		break
	end
	println("Done $i/300")
end
println("LassoTime",run," = ",toc()," seconds.")

### EXTRACT PREDICTION AND ISRS AND DATE FOR LASSO###
f = open(fileName*string(run)*"_Lasso"*dataFile*".csv", "w")
write(f, "Max Correlation,Test R^2,Vali R^2,k,gamma,"*expandedColNamesToString(colNames,1,0)*"\n")
writecsv(f,lassoOutput)
close(f)

#Spaced between 0 and half the SSTO since this would then get SSTO*absSumOfBeta which would force everything to 0
maxVal = findmax(abs.(standX'*standY))[1]
gammaArray = linspace(maxVal, 1, amountOfGammas)

#INITIALISE STORING ARRAYS
bSample = []
allCuts = []
signifBoolean = zeros(3)
best3Beta = zeros(3,bCols+3)
best3Beta[:,1] = -1e3
runCount = 1
curRSquared = -1e3
HC = cor(standX)

#BUILD AND SOLVE AALR MODEL
bSolved = []
warmstart = false
tic()
for i = 1:kmax
	for g = 1:amountOfGammas
		println("\nSolving for k = $i and gamma index = $g\n")
		bSolved, obj, objBound = buildAndSolveStage2(standX, standY, i, gammaArray[g], warmstart, bSolved, [], [], HC)

		#println("Cur index is: ",Int64(g+(i-1)*amountOfGammas))
		prevSolutions[Int64(g+(i-1)*amountOfGammas),:] = bSolved

		best3Beta = checkRSquared(standXVali, standYVali, bSolved, best3Beta, i, gammaArray[g])

		warmstart = true
	 end
end

#println(best3Beta)

#CHECK FOR SIGNIFICANCE (THIS INCLUDES SOLVING THE MODEL WITH FIXED z[] VARIABLES)
cuts = stageThree(best3Beta, standX, standY, [])

#Check for cuts
while !isempty(cuts)
	if runCount > 4
		break
	end
	best3Beta[:,1] = -1e3
	allCuts = vcat(allCuts,cuts[:,1:end-1])

	# Based on the cuts, ONLY resolve the problems that the cut affects.
	# Find the minimum number of variables in over all cuts
	minNumberOfVar = kmax
	for j=1:size(allCuts)[1]
		if countnz(allCuts[1,:]) < kmax
			minNumberOfVar = countnz(allCuts[1,:])
		end
	end

	#Rebuild and solve model with cuts
	for i = 1:kmax
		for g = 1:amountOfGammas
			bSolved = prevSolutions[Int64(g+(i-1)*amountOfGammas),:]
			if i >= minNumberOfVar
				println("\nSolving for k = $i and gamma index = $g\n")
				bSolved, obj, objBound = buildAndSolveStage2(standX, standY, i, gammaArray[g], warmstart, prevSolutions[Int64(g+(i-1)*amountOfGammas),:], allCuts, [], HC)
			end

			#prevSolutions[Int64(i+(i-1)*amountOfGammas),:] = bSolved

			best3Beta = checkRSquared(standXVali, standYVali, bSolved, best3Beta, i, gammaArray[g])

			warmstart = true
		 end
	end

	#Check significance
	cuts = stageThree(best3Beta, standX, standY, allCuts)

	runCount += 1
end
println("AALRTime",run," = ",toc()," seconds")
best3Beta = cat(2,signifBoolean,best3Beta)

rSquared = Array{Float64}(3,1)
#Convert into r^2 in stead

rSquared[1,1] = getRSquared(standXTest, standYTest, best3Beta[1,5:end])
rSquared[2,1] = getRSquared(standXTest, standYTest, best3Beta[2,5:end])
rSquared[3,1] = getRSquared(standXTest, standYTest, best3Beta[3,5:end])

maxCor = Array{Float64}(3,1)
maxCor[1,1] = findMaxCor(standX, best3Beta[1,6:end])
maxCor[2,1] = findMaxCor(standX, best3Beta[2,6:end])
maxCor[3,1] = findMaxCor(standX, best3Beta[3,6:end])

best3Beta = cat(2,rSquared,best3Beta)
best3Beta = cat(2,maxCor,best3Beta)

### EXTRACT PREDICTION AND ISRS AND DATE###
f = open(fileName*string(run)*"_"*dataFile*".csv", "w")
write(f, "Max Correlation,Test R^2,Significant,Vali R^2,k,gamma,"*expandedColNamesToString(colNames,1,0)*"\n")
writecsv(f,best3Beta)
close(f)

println("\nComplete\n")
#end

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/AALR.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/AALR.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
### A Pluto.jl notebook ###
# v0.19.2

using Markdown
using InteractiveUtils

# ╔═╡ 25d21156-d006-11eb-3655-61be76e7db93
begin
import Pkg;
	Pkg.activate("../../../")
end

# ╔═╡ 34aefa9a-e034-4293-b53d-cc1dc349ecc7
let
	using Revise
	using Unfold
	using DataFramesMeta
end

# ╔═╡ c73bf3ed-4af3-4e5a-997e-6399c9b85bbe
using StatsModels,PlutoUI,DataFrames,StatsPlots

# ╔═╡ d1e1de70-7e37-46a4-8ce2-3ca591bb84cf
include(pathof(Unfold)*"../../../test/test_utilities.jl")

# ╔═╡ 2e9ec6e0-03c6-41f6-bd3b-a44eaed5b553
data,evts = loadtestdata("test_case_4a");

# ╔═╡ 0b27e0bb-7aee-4672-8501-6096539c7ad7
evts.continuousA = rand(size(evts,1))

# ╔═╡ 97443eb8-ac67-4f8e-9e8a-eb09176c8516

# ╔═╡ 2a5d5493-b88d-473f-8d81-9e9240b3ffe0
#f  = @formula 0~1+conditionA+continuousA # 1
f  = @formula 0~1+continuousA # 1

# ╔═╡ 5d23b74c-988d-48bc-9c30-c2af0dbabdb4

# ╔═╡ 3f2b7a92-a064-4be6-8a39-cabaa2453777
basisfunction = Unfold.splinebasis(τ=(-1,1),sfreq=20,nsplines=10,name="basisA")

# ╔═╡ 8f32dbf3-8958-4e67-8e25-5f24d48058e1
Unfold.splinebasis

# ╔═╡ 7ef24e5a-63df-4e99-a021-7119d7e8d3bf
fir = firbasis(τ=(-1,1),sfreq=20,name="basisB")

# ╔═╡ 12c617c1-1994-4efe-b8bc-23fd2325c2ca
fir.kernel(2)

# ╔═╡ 6f7bc2f6-28de-40a5-83ee-c15cd6f74f6c
basisfunction.kernel(1)

# ╔═╡ e3d7ac15-e7a2-490b-abbc-aea3d7b4f4c7
size(Unfold.times(basisfunction))

# ╔═╡ 73022b4d-4edd-44e2-8cfa-26ded44a8921
size(Unfold.colnames(basisfunction))

# ╔═╡ d38f9cc9-409f-4c52-b0d2-d09312cc695a
size(basisfunction.colnames)

# ╔═╡ 786d0114-e10e-4e29-b6d7-bbe635c55f08
evts

# ╔═╡ c4227637-81b4-4ea7-a7b7-741860fef8c3
m = fit(UnfoldModel,Dict("eventA"=>(f,fir),"eventB"=>(f,basisfunction)),evts,data,eventcolumn="type");

# ╔═╡ cd19784c-686e-41c5-bd32-cdf793cecf03
res = coeftable(m)

# ╔═╡ 29e875d3-45b9-4a6a-aaa6-c6eaf35244c5
@df res plot(:time,:estimate,group=(:basisname,:coefname))

# ╔═╡ b669f55b-9121-4d66-9c82-f89bbfe0b2df
Unfold.coef(m)

# ╔═╡ 64cf7e78-2677-4551-9491-82c9f560ed61
	res

# ╔═╡ 02b8ff44-c514-492e-bf94-0cbf44431097
x = [1]

# ╔═╡ d964a5ca-7e63-41d5-bf97-526772d259e9
#[kernel(formTerm.basisfunction)(0) for formTerm in getfield.(designmatrix(m).formulas,:rhs)]

# ╔═╡ 1ee404e1-3039-410d-99de-ce713e751280
let
	basisnames = Unfold.get_basis_name(m)
	betaOut = []
	for (k,b) = enumerate(unique(basisnames))

		ix = basisnames .== b

		betaOut = vcat(betaOut,designmatrix(m).formulas[k].rhs.basisfunction.kernel(0)	* m.modelfit.estimate[:,ix]')
	end
	plot(Float64.(betaOut))
end

# ╔═╡ 5466ec0b-d28a-4be7-b686-cb5ad7fea92b
let
	basisnames = Unfold.get_basis_name(m)
	k  = 2
	b = unique(basisnames)[k]

	bf = designmatrix(m).formulas[k].rhs.basisfunction
	ix = basisnames .== b

	hcat(bf.kernel.([0,0])	...) * m.modelfit.estimate[:,ix]'
end

# ╔═╡ 033cf2e0-a3ca-40e3-8f27-e88e017ca8de
size(Unfold.times(basisfunction))

# ╔═╡ a8238bd7-0198-45d2-89bd-d1434c537b2a
	designmatrix(m).formulas[2].rhs.basisfunction.kernel(0)

# ╔═╡ fe5a826f-aede-4e91-b6cd-ee7fca33ce31
zeros()

# ╔═╡ 8b9d8e32-d2ca-49db-986f-851ffabc3a3e
function undoBasis(d,m)
	bname = unique(d.basisname) # the current coefficient basename
	@assert(length(bname)==1)
	@assert(length(unique(d.channel))==1)
	bnames = unique(Unfold.get_basis_name(m)) # all model basenames

	# find out which formula / basisfunction we have
	k = findfirst(bname .== bnames)
	bf = designmatrix(m).formulas[k].rhs.basisfunction

	estimate= bf.kernel.(0) *d.estimate
	@show size(bf.kernel.(0))
	@show size(Unfold.times(bf))
	return DataFrame(:time=>Unfold.times(bf),:estimate=>estimate)

end

# ╔═╡ 350a7f4b-f019-4bf6-8878-2cbb0cf82005
begin
	gd = groupby(res,[:basisname,:channel,:group,:coefname])
	a = combine(gd) do d
		return undoBasis(d,m)
	end
	a
end

# ╔═╡ a876f554-a034-461e-9522-490ca4bab5f8
@df a plot(:time,:estimate,group=(:basisname,:coefname))

# ╔═╡ 1d5b8b05-cdc6-4058-9bc1-bd52e6e3b444
Unfold.get_basis_name(m)

# ╔═╡ 111189ec-be04-4687-91ff-53b540a7b4db
Unfold.get_basis_name(m)

# ╔═╡ 4b572283-14cd-4a75-9ff2-f5a18f40228d
res

# ╔═╡ Cell order:
# ╠═25d21156-d006-11eb-3655-61be76e7db93
# ╠═34aefa9a-e034-4293-b53d-cc1dc349ecc7
# ╠═c73bf3ed-4af3-4e5a-997e-6399c9b85bbe
# ╠═d1e1de70-7e37-46a4-8ce2-3ca591bb84cf
# ╠═2e9ec6e0-03c6-41f6-bd3b-a44eaed5b553
# ╠═0b27e0bb-7aee-4672-8501-6096539c7ad7
# ╠═97443eb8-ac67-4f8e-9e8a-eb09176c8516
# ╠═2a5d5493-b88d-473f-8d81-9e9240b3ffe0
# ╠═5d23b74c-988d-48bc-9c30-c2af0dbabdb4
# ╠═3f2b7a92-a064-4be6-8a39-cabaa2453777
# ╠═8f32dbf3-8958-4e67-8e25-5f24d48058e1
# ╠═7ef24e5a-63df-4e99-a021-7119d7e8d3bf
# ╠═12c617c1-1994-4efe-b8bc-23fd2325c2ca
# ╠═6f7bc2f6-28de-40a5-83ee-c15cd6f74f6c
# ╠═e3d7ac15-e7a2-490b-abbc-aea3d7b4f4c7
# ╠═73022b4d-4edd-44e2-8cfa-26ded44a8921
# ╠═d38f9cc9-409f-4c52-b0d2-d09312cc695a
# ╠═786d0114-e10e-4e29-b6d7-bbe635c55f08
# ╠═c4227637-81b4-4ea7-a7b7-741860fef8c3
# ╠═cd19784c-686e-41c5-bd32-cdf793cecf03
# ╠═29e875d3-45b9-4a6a-aaa6-c6eaf35244c5
# ╠═b669f55b-9121-4d66-9c82-f89bbfe0b2df
# ╠═64cf7e78-2677-4551-9491-82c9f560ed61
# ╠═02b8ff44-c514-492e-bf94-0cbf44431097
# ╠═d964a5ca-7e63-41d5-bf97-526772d259e9
# ╠═1ee404e1-3039-410d-99de-ce713e751280
# ╠═5466ec0b-d28a-4be7-b686-cb5ad7fea92b
# ╠═033cf2e0-a3ca-40e3-8f27-e88e017ca8de
# ╠═a8238bd7-0198-45d2-89bd-d1434c537b2a
# ╠═350a7f4b-f019-4bf6-8878-2cbb0cf82005
# ╠═a876f554-a034-461e-9522-490ca4bab5f8
# ╠═fe5a826f-aede-4e91-b6cd-ee7fca33ce31
# ╠═8b9d8e32-d2ca-49db-986f-851ffabc3a3e
# ╠═1d5b8b05-cdc6-4058-9bc1-bd52e6e3b444
# ╠═111189ec-be04-4687-91ff-53b540a7b4db
# ╠═4b572283-14cd-4a75-9ff2-f5a18f40228d

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
## All this code breaks with JumP 0.19.

function solve(model::PowerOperationModel{T, F}; Solver = nothing) where {T <: AbstractOperationsModel, F <: Array{<:Function}}

    psmodel = model.model # JuMP model
    model_type = JuMP.ProblemTraits(psmodel)

    #==
    function ProblemTraits(m::Model; relaxation=false)
    int = !relaxation & any(c-> !(c == :Cont || c == :Fixed), m.colCat)
    qp = !isempty(m.obj.qvars1)
    qc = !isempty(m.quadconstr)
    nlp = m.nlpdata !== nothing
    soc = !isempty(m.socconstr)
    # will need to change this when we add support for arbitrary variable cones
    sdp = !isempty(m.sdpconstr) || !isempty(m.varCones)
    sos = !isempty(m.sosconstr)
    ProblemTraits(int, !(qp|qc|nlp|soc|sdp|sos), qp, qc, nlp, soc, sdp, sos, soc|sdp)
    end
    =#

    if Solver == nothing

        if (model_type.qc|model_type.qp|model_type.nlp)

            if model_type.int
                @error "The model is a Mixed Integer Non-Linear Problem, please define an appropiate solver manually using the Solver= argument"
            end

            JuMP.setsolver(psmodel, Ipopt.IpoptSolver())
            @warn "The model contains non-linear elements (QP, QC, NLP), by default the solver is Ipopt Solver"

        elseif model_type.lin & model_type.int & !(model_type.qc|model_type.qp|model_type.nlp)

            JuMP.setsolver(psmodel, Cbc.CbcSolver(logLevel = 1))
            @warn "The model is Mixed Integer Linear, the default is Cbc Solver"

        elseif model_type.lin & !(model_type.qc|model_type.qp|model_type.nlp)

            JuMP.setsolver(psmodel, Clp.ClpSolver(SolveType = 5, LogLevel = 4))

            @warn "The model is linear, by default the solver is Clp Solver"

        else

            @error "The PS Model is not a standard LP, MILP, QP, QC or NLP, please define an appropiate solver manually using the Solver argument"

        end

    else
        JuMP.setsolver(psmodel, Solver)
        @warn "A solver has been defined manually, this might break the results output function"

    end

    JuMP.build(psmodel)
    status = solve(psmodel)

    if status != :Optimal
        println(status)
        @error "Problem has no solution."
    end

end


text:
using Test, BufferedFiles, BufferedFiles.Libc

const curdir = pwd()
cd(@__DIR__)
if !isdir("testdir")
    mkdir("testdir")
end
cd("testdir")

@testset "BufferedWriteFile" begin
    for i = 1:1024
        f = open("testfile$i", om"w", i)
        @test f isa BufferedFiles.BufferedWriteFile
        @test position(f) == 0
        @test write(f, UInt8[1, 2, 3, 4, 5, 6, 7, 8]) == 8
        @test position(f) == 8
        @test write(f, 1) == 8
        @test position(f) == 16
        @test write(f, 0x01) == 1
        @test position(f) == 17
        close(f)
    end
end

@testset "BufferedReadFile" begin
    for i = 1:1024
        f = BufferedFiles.open("testfile$i", om"r", i)
        @test f isa BufferedFiles.BufferedReadFile
        @test eof(f) == false
        @test position(f) == 0
        @test read(f, 8) == UInt8[1, 2, 3, 4, 5, 6, 7, 8]
        @test position(f) == 8
        @test eof(f) == false
        @test read(f, Int) == 1
        @test position(f) == 16
        @test eof(f) == false
        @test read(f, UInt8) == 1
        @test position(f) == 17
        @test eof(f) == true
        @test_throws EOFError read(f, UInt8)
        close(f)
    end
end

cd("../")
rm("testdir", recursive = true)
cd(curdir)


text:
"""
    KSparseLMO{T}(K::Int, right_hand_side::T)

LMO for the K-sparse polytope:
```
C = B_1(τK) ∩ B_∞(τ)
```
with `τ` the `right_hand_side` parameter.
The LMO results in a vector with the K largest absolute values
of direction, taking values `-τ sign(x_i)`.
"""
struct KSparseLMO{T} <: LinearMinimizationOracle
    K::Int
    right_hand_side::T
end

function compute_extreme_point(lmo::KSparseLMO{T}, direction; kwargs...) where {T}
    K = min(lmo.K, length(direction))
    K_indices = sortperm(direction[1:K], by=abs, rev=true)
    K_values = direction[K_indices]
    for idx in K+1:length(direction)
        new_val = direction[idx]
        # new greater value: shift everything right
        if abs(new_val) > abs(K_values[1])
            K_values[2:end] .= K_values[1:end-1]
            K_indices[2:end] .= K_indices[1:end-1]
            K_indices[1] = idx
            K_values[1] = new_val
            # new value in the interior
        elseif abs(new_val) > abs(K_values[K])
            # NOTE: not out of bound since unreachable with K=1
            j = K - 1
            while abs(new_val) > abs(K_values[j])
                j -= 1
            end
            K_values[j+1:end] .= K_values[j:end-1]
            K_indices[j+1:end] .= K_indices[j:end-1]
            K_values[j+1] = new_val
            K_indices[j+1] = idx
        end
    end
    v = spzeros(T, length(direction))
    for (idx, val) in zip(K_indices, K_values)
        v[idx] = -lmo.right_hand_side * sign(val)
    end
    return v
end

function convert_mathopt(
    lmo::KSparseLMO{T},
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {T,OT}
    τ = lmo.right_hand_side
    n = dimension
    K = min(lmo.K, n)
    MOI.empty!(optimizer)
    x = MOI.add_variables(optimizer, n)
    tinf = MOI.add_variable(optimizer)
    MOI.add_constraint(optimizer, MOI.VectorOfVariables([tinf; x]), MOI.NormInfinityCone(n + 1))
    MOI.add_constraint(optimizer, tinf, MOI.LessThan(τ))
    t1 = MOI.add_variable(optimizer)
    MOI.add_constraint(optimizer, MOI.VectorOfVariables([t1; x]), MOI.NormOneCone(n + 1))
    MOI.add_constraint(optimizer, t1, MOI.LessThan(τ * K))
    return MathOptLMO(optimizer)
end

"""
    BirkhoffPolytopeLMO

The Birkhoff polytope encodes doubly stochastic matrices.
Its extreme vertices are all permutation matrices of side-dimension `dimension`.
"""
struct BirkhoffPolytopeLMO <: LinearMinimizationOracle end

function compute_extreme_point(
    ::BirkhoffPolytopeLMO,
    direction::AbstractMatrix{T};
    kwargs...,
) where {T}
    n = size(direction, 1)
    n == size(direction, 2) ||
        DimensionMismatch("direction should be square and matching BirkhoffPolytopeLMO dimension")
    res_mat = Hungarian.munkres(direction)
    m = spzeros(Bool, n, n)
    (rows, cols, vals) = SparseArrays.findnz(res_mat)
    @inbounds for i in eachindex(cols)
        m[rows[i], cols[i]] = vals[i] == 2
    end
    m = convert(SparseArrays.SparseMatrixCSC{Float64,Int64}, m)
    return m
end

function compute_extreme_point(
    lmo::BirkhoffPolytopeLMO,
    direction::AbstractVector{T};
    kwargs...,
) where {T}
    nsq = length(direction)
    n = isqrt(nsq)
    return compute_extreme_point(lmo, reshape(direction, n, n); kwargs...)
end

function convert_mathopt(
    ::BirkhoffPolytopeLMO,
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {OT}
    n = dimension
    MOI.empty!(optimizer)
    (x, _) = MOI.add_constrained_variables(optimizer, fill(MOI.Interval(0.0, 1.0), n * n))
    xmat = reshape(x, n, n)
    for idx in 1:n
        # column constraint
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), xmat[:, idx]), 0.0),
            MOI.EqualTo(1.0),
        )
        # row constraint
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), xmat[idx, :]), 0.0),
            MOI.EqualTo(1.0),
        )
    end
    return MathOptLMO(optimizer)
end


text:
"""
    couple_models(aero::LiftingLine, stru::GEBT, flap::LiftingLineFlaps, ctrl::Trim)

Create a coupled model using a lifting line aerodynamic model, a
geometrically exact beam theory model, a lifting line control surface model,
and a controller which maintains trimmed operating conditions.  This model introduces
additional parameters corresponding to the freestream air density ``\\rho_\\infty``,
followed by the external loads ``F_{x,i}, F_{y,i}, F_{z,i}, M_{x,i}, M_{y,i},
M_{z,i}`` or displacements ``u_{x,i}, u_{y,i}, u_{z,i}, \\theta_{x,i},
\\theta_{y,i}, \\theta_{z,i}`` for each node, followed by the constant distributed
loads ``f_{x,i}, f_{y,i}, f_{z,i}, m_{x,i}, m_{y,i}, m_{z,i}`` for each beam
element (excluding aerodynamic loads), followed by the body frame linear and
angular velocities ``u, v, w, p, q, r``, and the control surface deflections
which are not defined by the controller.

**NOTE: When using this model, the local frame for each beam element should be
oriented with the x-axis along the beam's axis, the y-axis forward, and the
z-axis normal to the surface**
"""
function couple_models(aero::LiftingLine, stru::GEBT, flap::LiftingLineFlaps,
    ctrl::Trim)

    return (aero, stru, flap, ctrl)
end

# --- traits --- #

function number_of_additional_parameters(aero::LiftingLine, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}, ctrl::Trim{NC}) where {NA,NF,NG,NC,TF,TG}

    return 1 + 6*length(stru.icol_point) + 6*length(stru.icol_elem) + 6 + NG - NC
end

function coupling_inplaceness(::Type{<:LiftingLine}, ::Type{<:GEBT},
    ::Type{<:LiftingLineFlaps}, ::Type{<:Trim})

    return InPlace()
end

function coupling_rate_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps{NF,NG,TF,TG}}, ::Type{<:Trim}) where {NA,NF,NG,TA,TF,TG}

    aero_model_types = TA.parameters
    flap_model_types = TF.parameters
    if all(islinear.(coupling_rate_jacobian_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Linear()
    else
        return Nonlinear()
    end
end

function coupling_state_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps{NF,NG,TF,TG}}, ::Type{<:Trim}) where {NA,NF,NG,TA,TF,TG}

    aero_model_types = TA.parameters
    flap_model_types = TF.parameters
    if all(islinear.(coupling_state_jacobian_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Linear()
    else
        return Nonlinear()
    end
end

coupling_parameter_jacobian_type(::Type{<:LiftingLine}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps}, ::Type{Trim}) = Nonlinear()

function coupling_time_gradient_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps{NF,NG,TF,TG}}, ::Type{<:Trim}) where {NA,NF,NG,TA,TF,TG}

    aero_model_types = TA.parameters
    flap_model_types = TF.parameters
    if all(isempty.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Empty()
    elseif all(iszero.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Zeros()
    elseif all(isconstant.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Constant()
    elseif all(islinear.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Linear()
    else
        return Nonlinear()
    end
end

# --- methods --- #

function get_coupling_inputs!(y, aero::LiftingLine{NA,TA}, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}, ctrl::Trim{NC}, dx, x, p, t) where {NA,NF,NG,NC,TA,TF,TG}

    # extract number of state variables, inputs, and parameters
    nxa = number_of_states(aero) # number of aerodynamic states
    nxs = number_of_states(stru) # number of structural states
    nxf = number_of_states(flap) # number of control surface states
    nxc = number_of_states(ctrl) # number of control surface states
    nya = number_of_inputs(aero) # number of aerodynamic inputs
    nys = number_of_inputs(stru) # number of structural inputs
    nyf = number_of_inputs(flap) # number of control surface inputs
    nyc = number_of_inputs(ctrl) # number of control surface inputs
    npa = number_of_parameters(aero) # number of aerodynamic parameters
    nps = number_of_parameters(stru) # number of structural parameters
    npf = number_of_parameters(flap) # number of control surface parameters
    npc = number_of_parameters(ctrl) # number of control surface parameters
    npadd = number_of_additional_parameters(aero, stru, flap, ctrl) # number of additional parameters

    # get indices for state variables, inputs, and parameters
    ixa = 1:nxa # indices of aerodynamic states
    iya = 1:nya # indices of aerodynamic inputs
    ipa = 1:npa # indices of aerodynamic parameters
    ixs = nxa + 1 : nxa + nxs # indices of structural states
    iys = nya + 1 : nya + nys # indices of structural inputs
    ips = npa + 1 : npa + nps # indices of structural parameters
    ixf = nxa + nxs + 1 : nxa + nxs + nxf # indices of control surface states
    iyf = nya + nys + 1 : nya + nys + nyf # indices of control surface inputs
    ipf = npa + nps + 1 : npa + nps + npf # indices of control surface parameters
    ixc = nxa + nxs + nxf + 1 : nxa + nxs + nxf + nxc # indices of control surface states
    iyc = nya + nys + nyf + 1 : nya + nys + nyf + nyc # indices of control surface inputs
    ipc = npa + nps + npf + 1 : npa + nps + npf + npc # indices of control surface parameters
    ipadd = npa + nps + npf + npc + 1 : npa + nps + npf + npc + npadd # indices of additional parameters
    ixas = state_indices(aero.models) # indices of aerodynamic states for each section
    iyas = input_indices(aero.models) # indices of aerodynamic inputs for each section
    ipas = parameter_indices(aero.models) # indices of aerodynamic parameters for each section
    ixfs = state_indices(flap.models) # indices of control surface states for each section
    iyfs = input_indices(flap.models) # indices of control surface inputs for each section
    ipfs = parameter_indices(flap.models) # indices of control surface parameters for each section

    # separate state variables, inputs, and parameters
    dxa = view(dx, ixa) # aerodynamic rates
    xa = view(x, ixa) # aerodynamic states
    ya = view(y, iya) # aerodynamic inputs
    pa = view(p, ipa) # aerodynamic parameters
    dxs = view(dx, ixs) # structural rates
    xs = view(x, ixs) # structural states
    ys = view(y, iys) # structural inputs
    ps = view(p, ips) # structural parameters
    dxf = view(dx, ixf) # control surface rates
    xf = view(x, ixf) # control surface states
    yf = view(y, iyf) # control surface inputs
    pf = view(p, ipf) # control surface parameters
    dxc = view(dx, ixc) # controller rates
    xc = view(x, ixc) # controller states
    yc = view(y, iyc) # controller inputs
    pc = view(p, ipc) # controller parameters
    padd = view(p, ipadd) # additional parameters for coupled model
    dxas = view.(Ref(dxa), ixas) # aerodynamic rates for each section
    xas = view.(Ref(xa), ixas) # aerodynamic states for each section
    yas = view.(Ref(ya), iyas) # aerodynamic inputs for each section
    pas = view.(Ref(pa), ipas) # aerodynamic parameters for each section
    dxfs = view.(Ref(dxf), ixfs) # control surface rates for each section
    xfs = view.(Ref(xf), ixfs) # control surface states for each section
    yfs = view.(Ref(yf), iyfs) # control surface inputs for each section
    pfs = view.(Ref(pf), ipfs) # control surface parameters for each section

    # number of points and elements
    npoint = length(stru.icol_point)
    nelem = length(stru.icol_elem)

    # controller state variables
    δc = xc

    # global parameters
    ρ = padd[1]
    ur = padd[1 + 6*npoint + 6*nelem + 1]
    vr = padd[1 + 6*npoint + 6*nelem + 2]
    wr = padd[1 + 6*npoint + 6*nelem + 3]
    pr = padd[1 + 6*npoint + 6*nelem + 4]
    qr = padd[1 + 6*npoint + 6*nelem + 5]
    rr = padd[1 + 6*npoint + 6*nelem + 6]
    δp = view(padd, 1 + 6*npoint + 6*nelem + 6 + 1 : 1 + 6*npoint + 6*nelem + 6 + NG - NC)

    # commanded control deflections
    T = promote_type(eltype(dx), eltype(x), eltype(p))
    dδ = zeros(SVector{NG,T})
    iδu = 1
    for i = 1:NG
        if i in ctrl.state_indices
            dδ = setindex(dδ, δc[iδu], i)
            iδu += 1
        else
            dδ = setindex(dδ, 0, i)
        end
    end

    # commanded control deflections
    T = promote_type(eltype(x), eltype(p))
    δ = zeros(SVector{NG,T})
    iδu = iδp = 1
    for i = 1:NG
        if i in ctrl.state_indices
            δ = setindex(δ, δc[iδu], i)
            iδu += 1
        else
            δ = setindex(δ, δp[iδp], i)
            iδp += 1
        end
    end

    # body linear and angular velocity
    V = SVector(ur, vr, wr)
    Ω = SVector(pr, qr, rr)

    # construct assembly from structural parameters
    assembly = gxbeam_assembly(ps, npoint, nelem, stru.start, stru.stop)

    # initialize total forces and moments
    Ftot = @SVector zeros(3)
    Mtot = @SVector zeros(3)

    # save prescribed point loads/displacements
    for ip = 1:npoint
        yoff = 6*(ip-1)
        poff = 1 + 6*(ip-1)
        y[iys[yoff+1]] = padd[poff+1]
        y[iys[yoff+2]] = padd[poff+2]
        y[iys[yoff+3]] = padd[poff+3]
        y[iys[yoff+4]] = padd[poff+4]
        y[iys[yoff+5]] = padd[poff+5]
        y[iys[yoff+6]] = padd[poff+6]
    end

    # add prescribed loads to total forces and moments
    for ip = 1:npoint
        # point location
        point = assembly.points[ip]
        # prescribed condition identities
        prescribed_forces = SVector{6}(view(stru.displacement, :, ip)) .== false
        # point loads
        poff = 1 + 6*(ip-1)
        F1 = ifelse(prescribed_forces[1], padd[poff+1], zero(padd[poff+1]))
        F2 = ifelse(prescribed_forces[2], padd[poff+2], zero(padd[poff+2]))
        F3 = ifelse(prescribed_forces[3], padd[poff+3], zero(padd[poff+3]))
        M1 = ifelse(prescribed_forces[4], padd[poff+4], zero(padd[poff+4]))
        M2 = ifelse(prescribed_forces[5], padd[poff+5], zero(padd[poff+5]))
        M3 = ifelse(prescribed_forces[6], padd[poff+6], zero(padd[poff+6]))
        Fp = SVector(F1, F2, F3)
        Mp = SVector(M1, M2, M3)
        Ftot += Fp
        Mtot += cross(point, Fp) + Mp
    end

    # calculate aerodynamic inputs, distributed loads, and element properties
    for i = 1:NA
        # models for this section
        section_aero = aero.models[i]
        section_stru = LiftingLineSection()
        section_flap = flap.models[i]
        section_ctrl = LiftingLineSectionControl()
        section_models = (section_aero, section_stru, section_flap, section_ctrl)

        # model dimensions for this section
        Nxai = number_of_states(section_aero)
        Nyai = number_of_inputs(section_aero)
        Npai = number_of_parameters(section_aero)
        Nxsi = number_of_states(section_stru)
        Nysi = number_of_inputs(section_stru)
        Npsi = number_of_parameters(section_stru)
        Nxfi = number_of_states(section_flap)
        Nyfi = number_of_inputs(section_flap)
        Npfi = number_of_parameters(section_flap)
        Nxci = number_of_states(section_ctrl)
        Nyci = number_of_inputs(section_ctrl)
        Npci = number_of_parameters(section_ctrl)

        # local section properties
        icol = stru.icol_elem[i]
        element = assembly.elements[i]
        u_elem = SVector(xs[icol], xs[icol+1], xs[icol+2]) # linear displacement
        θ_elem = SVector(xs[icol+3], xs[icol+4], xs[icol+5]) # angular displacement
        F_elem = SVector(xs[icol+6], xs[icol+7], xs[icol+8]) .* stru.force_scaling # internal force
        M_elem = SVector(xs[icol+9], xs[icol+10], xs[icol+11]) .* stru.force_scaling # internal moment
        P_elem = SVector(xs[icol+12], xs[icol+13], xs[icol+14]) .* stru.mass_scaling # linear momentum
        H_elem = SVector(xs[icol+15], xs[icol+16], xs[icol+17]) .* stru.mass_scaling # angular momentum

        dF_elem = SVector(dxs[icol+6], dxs[icol+7], dxs[icol+8]) .* stru.force_scaling # internal force
        dM_elem = SVector(dxs[icol+9], dxs[icol+10], dxs[icol+11]) .* stru.force_scaling # internal moment
        dP_elem = SVector(dxs[icol+12], dxs[icol+13], dxs[icol+14]) .* stru.mass_scaling # linear momentum
        dH_elem = SVector(dxs[icol+15], dxs[icol+16], dxs[icol+17]) .* stru.mass_scaling # angular momentum

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = GXBeam.rotation_parameter_scaling(θ_elem)
        θ_elem *= scaling # angular displacement (Wiener-Milenkovic parameters)

        # length and location of the element
        ΔL = element.L # element length
        pe = element.x + u_elem # location of element (in body frame)

        # local to body transformation
        CtCab = GXBeam.get_C(θ_elem)'*element.Cab

        # local structural to aerodynamic transformation
        R = @SMatrix [0 -1 0; 1 0 0; 0 0 1]

        # local freestream linear and angular velocities
        vi = -R*(CtCab'*V + GXBeam.element_linear_velocity(element, P_elem, H_elem)) # local linear freestream velocity
        ωi = R*GXBeam.element_angular_velocity(element, P_elem, H_elem) # local angular freestream velocity

        # local freestream linear and angular accelerations
        dvi = -R*GXBeam.element_linear_velocity(element, dP_elem, dH_elem)
        dωi = R*GXBeam.element_angular_velocity(element, dP_elem, dH_elem)

        # NOTE: We assume that the local beam frame y-axis is oriented towards
        # the leading edge and the z-axis is oriented up
        # NOTE: We assume the local aerodynamic frame is oriented with the x-axis
        # in the chordwise direction and the z-axis in the (airfoil) normal
        # direction.

        # section state variable rates
        dxai = SVector{Nxai}(dxas[i]) # aerodynamic state variables
        dxsi = vcat(dvi, dωi) # structural state variables
        dxfi = SVector{Nxfi}(dxfs[i]) # control surface state variables
        dxci = SMatrix{1,NG}(getindex.(flap.gains, i))* dδ # controller state variables
        dxi = vcat(dxai, dxsi, dxfi, dxci)

        # section state variables
        xai = SVector{Nxai}(xas[i]) # aerodynamic state variables
        xsi = vcat(vi, ωi) # structural state variables
        xfi = SVector{Nxfi}(xfs[i]) # control surface state variables
        xci = SMatrix{1,NG}(getindex.(flap.gains, i)) * δ # controller state variables
        xi = vcat(xai, xsi, xfi, xci)

        # section parameters
        pai = SVector{Npai}(pas[i]) # aerodynamic parameters
        psi = SVector(ρ) # structural parameters
        pfi = SVector{Npfi}(pfs[i]) # control surface parameters
        pci = SVector{0,Float64}() # controller parameters
        pi = vcat(pai, psi, pfi, pci)

        # section inputs
        yi = get_coupling_inputs(section_models, dxi, xi, pi, t)

        # separate inputs
        yai = view(yi, 1:Nyai)
        ysi = view(yi, Nyai + 1 : Nyai + Nysi)
        yfi = view(yi, Nyai + Nysi + 1 : Nyai + Nysi + Nyfi)
        yci = view(yi, Nyai + Nysi + Nyfi + 1 : Nyai + Nysi + Nyfi + Nyci)

        # save local inputs
        y[iya[iyas[i]]] .= yai # aerodynamic inputs
        y[iyf[iyfs[i]]] .= yfi # control surface inputs

        # section aerodynamic loads (in body frame)
        fi = CtCab*R'*SVector(yi[Nyai+1], yi[Nyai+2], yi[Nyai+3])
        mi = CtCab*R'*SVector(yi[Nyai+4], yi[Nyai+5], yi[Nyai+6])

        # add constant distributed loads (in body frame)
        poff = 1 + 6*npoint + 6*(i-1)
        fi += SVector(padd[poff+1], padd[poff+2], padd[poff+3])
        mi += SVector(padd[poff+4], padd[poff+5], padd[poff+6])

        # save distributed loads for this element (in body frame)
        yoff = 6*npoint + 6*(i-1)
        y[iys[yoff+1 : yoff+3]] = fi
        y[iys[yoff+4 : yoff+6]] = mi

        # add distributed loads to total forces and moments
        Ftot += ΔL*fi
        Mtot += cross(pe, ΔL*fi) + ΔL*mi

    end

    # save body frame linear/angular velocities
    yoff = 6*npoint + 6*nelem
    y[iys[yoff+1]] = ur
    y[iys[yoff+2]] = vr
    y[iys[yoff+3]] = wr
    y[iys[yoff+4]] = pr
    y[iys[yoff+5]] = qr
    y[iys[yoff+6]] = rr

    # save trim model inputs
    loads = SVector(Ftot[1], Ftot[2], Ftot[3], Mtot[1], Mtot[2], Mtot[3])

    for i = 1:NC
        idx = ctrl.force_indices[i]
        y[iyc[i]] = loads[idx]
    end

    return y
end

# --- performance overloads --- #

# TODO

# --- convenience methods --- #

function set_additional_parameters!(padd, model::LiftingLine, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}, ctrl::Trim{NC}; rho, point_conditions,
    element_loads, u, v, w, p, q, r, delta) where {NF,NG,NC,TF,TG}

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    padd[1] = rho

    for ip = 1:np
        padd[1+6*(ip-1)+1] = point_conditions[6*(ip-1)+1]
        padd[1+6*(ip-1)+2] = point_conditions[6*(ip-1)+2]
        padd[1+6*(ip-1)+3] = point_conditions[6*(ip-1)+3]
        padd[1+6*(ip-1)+4] = point_conditions[6*(ip-1)+4]
        padd[1+6*(ip-1)+5] = point_conditions[6*(ip-1)+5]
        padd[1+6*(ip-1)+6] = point_conditions[6*(ip-1)+6]
    end

    for ie = 1:ne
        padd[1+6*np+6*(ie-1)+1] = element_loads[6*(ie-1)+1]
        padd[1+6*np+6*(ie-1)+2] = element_loads[6*(ie-1)+2]
        padd[1+6*np+6*(ie-1)+3] = element_loads[6*(ie-1)+3]
        padd[1+6*np+6*(ie-1)+4] = element_loads[6*(ie-1)+4]
        padd[1+6*np+6*(ie-1)+5] = element_loads[6*(ie-1)+5]
        padd[1+6*np+6*(ie-1)+6] = element_loads[6*(ie-1)+6]
    end

    padd[1 + 6*np + 6*ne + 1] = u
    padd[1 + 6*np + 6*ne + 2] = v
    padd[1 + 6*np + 6*ne + 3] = w
    padd[1 + 6*np + 6*ne + 4] = p
    padd[1 + 6*np + 6*ne + 5] = q
    padd[1 + 6*np + 6*ne + 6] = r
    padd[1 + 6*np + 6*ne + 6 + 1 : 1 + 6*np + 6*ne + 6 + NG - NC] = delta

    return padd
end

function separate_additional_parameters(model::LiftingLine, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}, ctrl::Trim{NC}, padd) where {NF,NG,NC,TF,TG}

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    rho = padd[1]

    point_conditions = reshape(view(padd, 1 + 1 : 1 + 6*np), 6, np)

    element_loads = reshape(view(padd, 1 + 6*np + 1 : 1 + 6*np + 6*ne), 6, ne)

    u = padd[1 + 6*np + 6*ne + 1]
    v = padd[1 + 6*np + 6*ne + 2]
    w = padd[1 + 6*np + 6*ne + 3]
    p = padd[1 + 6*np + 6*ne + 4]
    q = padd[1 + 6*np + 6*ne + 5]
    r = padd[1 + 6*np + 6*ne + 6]
    delta = view(padd, 1 + 6*np + 6*ne + 6 + 1 : 1 + 6*np + 6*ne + 6 + NG - NC)

    return (rho = rho, point_conditions = point_conditions, element_loads = element_loads,
        u = u, v = v, w = w, p = p, q = q, r = r, delta = delta)
end