text:
# airfoil methods

using Printf: @printf

export afeval, AFType, SimpleAF, AlphaAF, AlphaMachAF, AlphaReAF, AlphaReMachAF
export read_AlphaAF, read_AlphaReAF, read_AlphaMachAF, read_AlphaReMachAF, write_af
export mach_correction, MachCorrection, PrandtlGlauert
export re_correction, ReCorrection, SkinFriction, LaminarSkinFriction, TurbulentSkinFriction
export rotation_correction, RotationCorrection, DuSeligEggers
export tip_correction, TipCorrection, PrandtlTip, PrandtlTipHub
export viterna

# TODO: check for repeat alpha as that will mess up fit.

############# file parsing ##############
"""
A basic airfoil file format.  `nheader` is the number of header lines,
which will be skipped.  For one Reynolds/Mach number.
Additional data like cm is optional but will be ignored.

format:

informational header\n
Re\n
Mach\n
alpha1 cl1 cd1 ...\n
alpha2 cl2 cd2\n
alpha3 cl3 cd3\n
...
"""
function parsefile(filename, radians)

    alpha = Float64[]
    cl = Float64[]
    cd = Float64[]
    info = ""
    Re = 1.0
    Mach = 1.0

    open(filename) do f

        # skip header
        info = readline(f)
        Re = parse(Float64, readline(f))
        Mach = parse(Float64, readline(f))

        for line in eachline(f)
            parts = split(line)
            push!(alpha, parse(Float64, parts[1]))
            push!(cl, parse(Float64, parts[2]))
            push!(cd, parse(Float64, parts[3]))
        end
    end

    if !radians
        alpha *= pi/180
    end

    return info, Re, Mach, alpha, cl, cd
end

function writefile(filename, info, Re, Mach, alpha, cl, cd, radians)
    open(filename, "w") do f
        @printf(f, "%s\n", info)
        @printf(f, "%.17g\n", Re)
        @printf(f, "%.17g\n", Mach)

        factor = 1.0
        if !radians
            factor = 180/pi
        end

        for i = 1:length(alpha)
            @printf(f, "%.17g\t%.17g\t%.17g\n", alpha[i]*factor, cl[i], cd[i])
        end
    end

    return nothing
end
###################################

########### Airfoil Data Types ###############

# -- interface -----
abstract type AFType end

"""
    afeval(af::AFType, alpha, Re, Mach)

Evaluate airfoil aerodynamic performance

**Arguments**
- `af::AFType or Function`: dispatch on AFType or if function call: `cl, cd = af(alpha, Re, Mach)`
- `alpha::Float64`: angle of attack in radians
- `Re::Float64`: Reynolds number
- `Mach::Float64`: Mach number

**Returns**
- `cl::Float64`: lift coefficient
- `cd::Float64`: drag coefficient
"""
function afeval(af::AFType, alpha, Re, Mach)
    @error "need to choose an AFType"
end
# ---------------

# --- Function Type -----
function afeval(af::Function, alpha, Re, Mach)
    return af(alpha, Re, Mach)
end
# ---------------

# --- Simple paramterization ------

"""
    SimpleAF(m, alpha0, clmax, clmin, cd0, cd2)

A simple parameterized lift and drag curve.
- `cl = m (alpha - alpha0)` (capped by clmax/clmin)
- `cd = cd0 + cd2 * cl^2`

**Arguments**
- `m::Float64`: lift curve slope
- `alpha0::Float64`: zero-lift angle of attack
- `clmax::Float64`: maximum lift coefficient
- `clmin::Float64`: minimum lift coefficient
- `cd0::Float64`: zero lift drag
- `cd2::Float64`: quadratic drag term
"""
struct SimpleAF{TF} <: AFType
    m::TF
    alpha0::TF
    clmax::TF
    clmin::TF
    cd0::TF
    cd2::TF
end

function afeval(af::SimpleAF, alpha, Re, Mach)
    cl = af.m*(alpha - af.alpha0)
    cl = min(cl, af.clmax)
    cl = max(cl, af.clmin)

    cd = af.cd0 + af.cd2*cl^2

    return cl, cd
end

# ---------------

# --- Spline Type (just alpha)

"""
    AlphaAF(alpha, cl, cd, info, Re, Mach)
    AlphaAF(alpha, cl, cd, info, Re=0.0, Mach=0.0)
    AlphaAF(alpha, cl, cd, info="CCBlade generated airfoil", Re=0.0, Mach=0.0)
    AlphaAF(filename::String; radians=true)

Airfoil data that varies with angle of attack.  Data is fit with an Akima spline.

**Arguments**
- `alpha::Vector{Float64}`: angles of attack
- `cl::Vector{Float64}`: corresponding lift coefficients
- `cd::Vector{Float64}`: corresponding drag coefficients
- `info::String`: a description of this airfoil data (just informational)
- `Re::Float64`: Reynolds number data was taken at (just informational)
- `Mach::Float64`: Mach number data was taken at (just informational)

or

a file

**Arguments**
- `filename::String`: name/path of file to read in
- `radians::Bool`: true if angle of attack in file is given in radians
"""
struct AlphaAF{TF, TS, TSP} <: AFType
    alpha::Vector{TF}
    cl::Vector{TF}
    cd::Vector{TF}
    info::TS # not used except for info in file
    Re::TF  # not used except for info in file
    Mach::TF # not used except for info in file
    clspline::TSP  # cache these for speed
    cdspline::TSP
end

function AlphaAF(alpha, cl, cd, info, Re, Mach)
    afcl = FLOWMath.Akima(alpha, cl)
    afcd = FLOWMath.Akima(alpha, cd)

    return AlphaAF(alpha, cl, cd, info, Re, Mach, afcl, afcd)
end

AlphaAF(alpha, cl, cd, info) = AlphaAF(alpha, cl, cd, info, 0.0, 0.0)
AlphaAF(alpha, cl, cd) = AlphaAF(alpha, cl, cd, "CCBlade generated airfoil", 0.0, 0.0)

function AlphaAF(filename::String; radians=true)
    info, Re, Mach, alpha, cl, cd = parsefile(filename, radians)
    return AlphaAF(alpha, cl, cd, info, Re, Mach)
end

function afeval(af::AlphaAF, alpha, Re, Mach)
    return af.clspline(alpha), af.cdspline(alpha)
end

"""
    write_af(filename(s), af::AFType; radians=true)

Write airfoil data to file

**Arguments**
- `filename(s)::String or Vector{String} or Matrix{String}`: name/path of file to write to
- `af::AFType`: writing is dispatched based on type (AlphaAF, AlphaReAF, etc.)
- `radians::Bool`: true if you want angle of attack to be written in radians
"""
function write_af(filename, af::AlphaAF; radians=true)
    writefile(filename, af.info, af.Re, af.Mach, af.alpha, af.cl, af.cd, radians)
    return nothing
end
# ---------------

# ----- higher order splines

"""
    AlphaReAF(alpha, Re, cl, cd, info, Mach)
    AlphaReAF(alpha, Re, cl, cd, info)
    AlphaReAF(alpha, Re, cl, cd)
    read_AlphaReAF(filenames::Vector{String}; radians=true)

Airfoil data that varies with angle of attack and Reynolds number.
Data is fit with a recursive Akima spline.

**Arguments**
- `alpha::Vector{Float64}`: angles of attack
- `Re::Vector{Float64}`: Reynolds numbers
- `cl::Matrix{Float64}`: lift coefficients where cl[i, j] corresponds to alpha[i], Re[j]
- `cd::Matrix{Float64}`: drag coefficients where cd[i, j] corresponds to alpha[i], Re[j]
- `info::String`: a description of this airfoil data (just informational)
- `Mach::Float64`: Mach number data was taken at (just informational)

or

filenames with one file per Reynolds number.

**Arguments**
- `filenames::Vector{String}`: name/path of files to read in, each at a different Reynolds number in ascending order
- `radians::Bool`: true if angle of attack in file is given in radians
"""
struct AlphaReAF{TF, TS} <: AFType
    alpha::Vector{TF}
    Re::Vector{TF}
    cl::Matrix{TF}
    cd::Matrix{TF}
    info::TS # not used except for info in file
    Mach::TF # not used except for info in file
end

AlphaReAF(alpha, Re, cl, cd, info) = AlphaReAF(alpha, Re, cl, cd, info, 0.0)
AlphaReAF(alpha, Re, cl, cd) = AlphaReAF(alpha, Re, cl, cd, "CCBlade generated airfoil", 0.0)

function AlphaReAF(filenames::Vector{String}; radians=true)

    info, Re1, Mach, alpha, cl1, cd1 = parsefile(filenames[1], radians)  # assumes common alpha across files, also common info and common Mach
    nalpha = length(alpha)
    ncond = length(filenames)

    cl = Array{Float64}(undef, nalpha, ncond)
    cd = Array{Float64}(undef, nalpha, ncond)
    Re = Array{Float64}(undef, ncond)
    cl[:, 1] = cl1
    cd[:, 1] = cd1
    Re[1] = Re1

    # iterate over remaining files
    for i = 2:ncond
        _, Rei, _, _, cli, cdi = parsefile(filenames[i], radians)
        cl[:, i] = cli
        cd[:, i] = cdi
        Re[i] = Rei
    end

    return AlphaReAF(alpha, Re, cl, cd, info, Mach)
end

function afeval(af::AlphaReAF, alpha, Re, Mach)

    cl = FLOWMath.interp2d(FLOWMath.akima, af.alpha, af.Re, af.cl, [alpha], [Re])[1]
    cd = FLOWMath.interp2d(FLOWMath.akima, af.alpha, af.Re, af.cd, [alpha], [Re])[1]

    return cl, cd
end

function write_af(filenames, af::AlphaReAF; radians=true)
    for i = 1:length(af.Re)
        writefile(filenames[i], af.info, af.Re[i], af.Mach, af.alpha, af.cl[:, i], af.cd[:, i], radians)
    end
    return nothing
end

# ---------------

"""
    AlphaMachAF(alpha, Mach, cl, cd, info, Re)
    AlphaMachAF(alpha, Mach, cl, cd, info)
    AlphaMachAF(alpha, Mach, cl, cd)
    AlphaMachAF(filenames::Vector{String}; radians=true)

Airfoil data that varies with angle of attack and Mach number.
Data is fit with a recursive Akima spline.

**Arguments**
- `alpha::Vector{Float64}`: angles of attack
- `Mach::Vector{Float64}`: Mach numbers
- `cl::Matrix{Float64}`: lift coefficients where cl[i, j] corresponds to alpha[i], Mach[j]
- `cd::Matrix{Float64}`: drag coefficients where cd[i, j] corresponds to alpha[i], Mach[j]
- `info::String`: a description of this airfoil data (just informational)
- `Re::Float64`: Reynolds number data was taken at (just informational)

or

filenames with one file per Mach number.

**Arguments**
- `filenames::Vector{String}`: name/path of files to read in, each at a different Mach number in ascending order
- `radians::Bool`: true if angle of attack in file is given in radians
"""
struct AlphaMachAF{TF, TS} <: AFType
    alpha::Vector{TF}
    Mach::Vector{TF}
    cl::Matrix{TF}
    cd::Matrix{TF}
    info::TS # not used except for info in file
    Re::TF # not used except for info in file
end

AlphaMachAF(alpha, Mach, cl, cd, info) = AlphaMachAF(alpha, Mach, cl, cd, info, 0.0)
AlphaMachAF(alpha, Mach, cl, cd) = AlphaMachAF(alpha, Mach, cl, cd, "CCBlade generated airfoil", 0.0)

function AlphaMachAF(filenames::Vector{String}; radians=true)

    info, Re, Mach1, alpha, cl1, cd1 = parsefile(filenames[1], radians)  # assumes common alpha across files, also common info and common Re
    nalpha = length(alpha)
    ncond = length(filenames)

    cl = Array{Float64}(undef, nalpha, ncond)
    cd = Array{Float64}(undef, nalpha, ncond)
    Mach = Array{Float64}(undef, ncond)
    cl[:, 1] = cl1
    cd[:, 1] = cd1
    Mach[1] = Mach1

    # iterate over remaining files
    for i = 2:ncond
        _, _, Machi, _, cli, cdi = parsefile(filenames[i], radians)
        cl[:, i] = cli
        cd[:, i] = cdi
        Mach[i] = Machi
    end

    return AlphaMachAF(alpha, Mach, cl, cd, info, Re)
end

function afeval(af::AlphaMachAF, alpha, Re, Mach)

    cl = FLOWMath.interp2d(FLOWMath.akima, af.alpha, af.Mach, af.cl, [alpha], [Mach])[1]
    cd = FLOWMath.interp2d(FLOWMath.akima, af.alpha, af.Mach, af.cd, [alpha], [Mach])[1]

    return cl, cd
end

function write_af(filenames, af::AlphaMachAF; radians=true)
    for i = 1:length(af.Mach)
        writefile(filenames[i], af.info, af.Re, af.Mach[i], af.alpha, af.cl[:, i], af.cd[:, i], radians)
    end
    return nothing
end

# ---------------

"""
    AlphaReMachAF(alpha, Re, Mach, cl, cd, info)
    AlphaReMachAF(alpha, Re, Mach, cl, cd)
    AlphaReMachAF(filenames::Matrix{String}; radians=true)

Airfoil data that varies with angle of attack, Reynolds number, and Mach number.
Data is fit with a recursive Akima spline.

**Arguments**
- `alpha::Vector{Float64}`: angles of attack
- `Re::Vector{Float64}`: Reynolds numbers
- `Mach::Vector{Float64}`: Mach numbers
- `cl::Array{Float64}`: lift coefficients where cl[i, j, k] corresponds to alpha[i], Re[j], Mach[k]
- `cd::Array{Float64}`: drag coefficients where cd[i, j, k] corresponds to alpha[i], Re[j], Mach[k]
- `info::String`: a description of this airfoil data (just informational)

or files with one per Re/Mach combination

**Arguments**
- `filenames::Matrix{String}`: name/path of files to read in.  filenames[i, j] corresponds to Re[i] Mach[j] with Reynolds number and Mach number in ascending order.
- `radians::Bool`: true if angle of attack in file is given in radians
"""
struct AlphaReMachAF{TF, TS} <: AFType
    alpha::Vector{TF}
    Re::Vector{TF}
    Mach::Vector{TF}
    cl::Array{TF}
    cd::Array{TF}
    info::TS
end

AlphaReMachAF(alpha, Re, Mach, cl, cd) = AlphaMachAF(alpha, Re, Mach, cl, cd, "CCBlade generated airfoil")

function AlphaReMachAF(filenames::Matrix{String}; radians=true)

    info, Re1, Mach1, alpha, cl1, cd1 = parsefile(filenames[1, 1], radians)  # assumes common alpha and info across files
    nalpha = length(alpha)
    nRe = length(Re1)
    nMach = length(Mach1)

    cl = Array{Float64}(undef, nalpha, nRe, nMach)
    cd = Array{Float64}(undef, nalpha, nRe, nMach)
    Re = Array{Float64}(undef, nRe)
    Mach = Array{Float64}(undef, nMach)

    for j = 1:nMach
        for i = 1:nRe
            _, Rei, Machj, _, clij, cdij = parsefile(filenames[i, j], radians)
            cl[:, i, j] = clij
            cd[:, i, j] = cdij
            Re[i] = Rei
            Mach[j] = Machj  # NOTE: probably should add check to prevent user error here.
        end
    end

    return AlphaReMachAF(alpha, Re, Mach, cl, cd, info)
end

function afeval(af::AlphaReMachAF, alpha, Re, Mach)
    cl = FLOWMath.interp3d(FLOWMath.akima, af.alpha, af.Re, af.Mach, af.cl, [alpha], [Re], [Mach])[1]
    cd = FLOWMath.interp3d(FLOWMath.akima, af.alpha, af.Re, af.Mach, af.cd, [alpha], [Re], [Mach])[1]

    return cl, cd
end

function write_af(filenames, af::AlphaReMachAF; radians=true)
    nre = length(af.Re)
    nm = length(af.Mach)

    for i = 1:nre
        for j = 1:nm
            writefile(filenames[i, j], af.info, af.Re[i], af.Mach[j], af.alpha, af.cl[:, i, j], af.cd[:, i, j], radians)
        end
    end

    return nothing
end

#####################################

########## Mach Corrections ################

abstract type MachCorrection end

struct PrandtlGlauert <: MachCorrection end

"""
    mach_correction(::MachCorrection, cl, cd, Mach)

Mach number correction for lift/drag coefficient

**Arguments**
- `mc::MachCorrection`: used for dispatch
- `cl::Float64`: lift coefficient before correction
- `cd::Float64`: drag coefficient before correction
- `Mach::Float64`: Mach number

**Returns**
- `cl::Float64`: lift coefficient after correction
- `cd::Float64`: drag coefficient after correction
"""
function mach_correction(::MachCorrection, cl, cd, Mach)
    return cl, cd
end

"""
    mach_correction(::PrandtlGlauert, cl, cd, Mach)

Prandtl/Glauert Mach number correction for lift coefficient
"""
function mach_correction(::PrandtlGlauert, cl, cd, Mach)
    beta = sqrt(1 - Mach^2)
    cl /= beta
    return cl, cd
end

########## Reynolds corrections ##############

abstract type ReCorrection end

"""
    re_correction(re::ReCorrection, cl, cd, Re)

Reynolds number correction for lift/drag coefficient

**Arguments**
- `re::ReCorrection`: used for dispatch
- `cl::Float64`: lift coefficient before correction
- `cd::Float64`: drag coefficient before correction
- `Re::Float64`: Reynolds number

**Returns**
- `cl::Float64`: lift coefficient after correction
- `cd::Float64`: drag coefficient after correction
"""
function re_correction(::ReCorrection, cl, cd, Re)
    return cl, cd
end

"""
    SkinFriction(Re0, p)

Skin friction model for a flat plate.
`cd *= (Re0 / Re)^p`

**Arguments**
- `Re0::Float64`: reference Reynolds number (i.e., no corrections at this number)
- `p::Float64`: exponent in flat plate model.  0.5 for laminar (Blasius solution), ~0.2 for fully turbulent (Schlichting empirical fit)
"""
struct SkinFriction{TF} <: ReCorrection
    Re0::TF
    p::TF
end
LaminarSkinFriction(Re0) = SkinFriction(Re0, 0.5)
TurbulentSkinFriction(Re0) = SkinFriction(Re0, 0.2)

"""
    re_correction(sf::SkinFriction, cl, cd, Re)

Skin friction coefficient correction based on flat plat drag increases with Reynolds number.
"""
function re_correction(sf::SkinFriction, cl, cd, Re)
    cd *= (sf.Re0 / Re)^sf.p
    return cl, cd
end

############ Rotation Correction ###################

abstract type RotationCorrection end

"""
    rotation_correction(rc::RotationCorrection, cl, cd, cr, rR, tsr, alpha, phi=alpha, alpha_max_corr=30*pi/180)

Rotation correction (3D stall delay).

**Arguments**
- `rc::RotationCorrection`: used for dispatch
- `cl::Float64`: lift coefficient before correction
- `cd::Float64`: drag coefficient before correction
- `cr::Float64`: local chord / local radius
- `rR::Float64`: local radius / tip radius
- `tsr::Float64`: local tip speed ratio (Omega r / Vinf)
- `alpha::Float64`: local angle of attack
- `phi::Float64`: local inflow angles (defaults to angle of attack is precomputing since it is only known for on-the-fly computations)
- `alpha_max_corr::Float64`: angle of attack for maximum correction (tapers off to zero by 90 degrees)

**Returns**
- `cl::Float64`: lift coefficient after correction
- `cd::Float64`: drag coefficient after correction
"""
function rotation_correction(rc::RotationCorrection, cl, cd, cr, rR, tsr, alpha, phi=alpha, alpha_max_corr=30*pi/180)
    return cl, cd
end

"""
    DuSeligEggers(a, b, d, m, alpha0)
    DuSeligEggers(a=1.0, b=1.0, d=1.0, m=2*pi, alpha0=0.0)  # uses defaults

DuSelig correction for lift an Eggers correction for drag.

**Arguments**
- `a, b, d::Float64`: parameters in Du-Selig paper.  Normally just 1.0 for each.
- `m::Float64`: lift curve slope.  Defaults to 2 pi for zero argument version.
- `alpha0::Float64`: zero-lift angle of attack.  Defaults to 0 for zero argument version.
"""
struct DuSeligEggers{TF} <: RotationCorrection
    a::TF
    b::TF
    d::TF
    m::TF
    alpha0::TF
end

DuSeligEggers() = DuSeligEggers(1.0, 1.0, 1.0, 2*pi, 0.0)

function rotation_correction(du::DuSeligEggers, cl, cd, cr, rR, tsr, alpha, phi=alpha, alpha_max_corr=30*pi/180)

    # Du-Selig correction for lift
    Lambda = tsr / sqrt(1 + tsr^2)
    expon = du.d / (Lambda * rR)
    fcl = 1.0/du.m*(1.6*cr/0.1267*(du.a-cr^expon)/(du.b+cr^expon)-1)

    # linear lift line
    cl_linear = du.m*(alpha - du.alpha0)

    # adjustment for max correction
    amax = atan(1/0.12) - 5*pi/180  # account for singularity in Eggers (not pi/2)
    if abs(real(alpha)) >= amax
        adj = 0.0
    elseif abs(real(alpha)) > alpha_max_corr
        adj = ((amax-abs_cs_safe(alpha))/(amax-alpha_max_corr))^2
    else
        adj = 1.0
    end

    # increment in cl
    deltacl = adj*fcl*(cl_linear - cl)
    cl += deltacl

    # Eggers correction for drag
    deltacd = deltacl * (sin(phi) - 0.12*cos(phi))/(cos(phi) + 0.12*sin(phi))  # note that we can actually use phi instead of alpha as is done in airfoilprep.py b/c this is done at each iteration
    cd += deltacd

    return cl, cd
end

# ------ Extract Linear Portion of Curve -----------

function linearliftcoeff(af::T, Re, Mach) where T <: Union{AlphaAF, AlphaReAF, AlphaMachAF, AlphaReMachAF}

    clwrap(alpha) = afeval(af, alpha, Re, Mach)[1]
    alpha = range(-2*pi/180, 5*pi/180, length=20)
    cl = clwrap.(alpha)

    return linearliftregression(alpha, cl)
end

function linearliftcoeff(alpha, cl, alphamin=-2*pi/180, alphamax=5*pi/180)

    idxmin = findfirst(alpha .> alphamin)
    idxmax = findfirst(alpha .> alphamax)
    alphasub = @view alpha[idxmin:idxmax]
    clsub = @view cl[idxmin:idxmax]
    return linearliftregression(alphasub, clsub)
end

function linearliftregression(alphasub, clsub)
    # linear regression
    coeff = [alphasub ones(length(alphasub))] \ clsub
    m = coeff[1]  # lift curve slope
    alpha0 = -coeff[2] / m  # zero-lift angle of attack

    return m, alpha0
end

############## tip correction ################

abstract type TipCorrection end

"""
    tip_correction(::TipCorrection, r, Rhub, Rtip, phi, B)

Tip corrections for 3D flow.

**Arguments**
- `tc::TipCorrection`: used for dispatch
- `r::Float64`: local radius
- `Rhub::Float64`: hub radius
- `Rtip::Float64`: tip radius
- `phi::Float64`: inflow angle
- `B::Integer`: number of blades

**Returns**
- `F::Float64`: tip loss factor to multiple against loads.
"""
function tip_correction(::TipCorrection, cl, cd, Mach)
    return 1.0
end

"""
    PrandtlTip()

Standard Prandtl tip loss correction.
"""
struct PrandtlTip <: TipCorrection end

"""
    PrandtlTipHub()

Standard Prandtl tip loss correction plus hub loss correction of same form.
"""
struct PrandtlTipHub <: TipCorrection end

function tip_correction(::PrandtlTip, r, Rhub, Rtip, phi, B)

    asphi = abs_cs_safe(sin(phi))
    factortip = B/2.0*(Rtip/r - 1)/asphi
    F = 2.0/pi*acos(exp(-factortip))

    return F
end

function tip_correction(::PrandtlTipHub, r, Rhub, Rtip, phi, B)

    # Prandtl's tip and hub loss factor
    asphi = abs_cs_safe(sin(phi))
    factortip = B/2.0*(Rtip/r - 1)/asphi
    Ftip = 2.0/pi*acos(exp(-factortip))
    factorhub = B/2.0*(r/Rhub - 1)/asphi
    Fhub = 2.0/pi*acos(exp(-factorhub))
    F = Ftip * Fhub

    return F
end

############## Extrapolation ##################

"""
    viterna(alpha, cl, cd, cr75, nalpha=50)

Viterna extrapolation.  Follows Viterna paper and somewhat follows NREL version of AirfoilPrep, but with some modifications for better robustness and smoothness.

**Arguments**
- `alpha::Vector{Float64}`: angles of attack
- `cl::Vector{Float64}`: correspnding lift coefficients
- `cd::Vector{Float64}`: correspnding drag coefficients
- `cr75::Float64`: chord/Rtip at 75% Rtip
- `nalpha::Int64`: number of discrete points (angles of attack) to include in extrapolation

**Returns**
- `alpha::Vector{Float64}`: angle of attack from -pi to pi
- `cl::Vector{Float64}`: correspnding extrapolated lift coefficients
- `cd::Vector{Float64}`: correspnding extrapolated drag coefficients
"""
function viterna(alpha, cl, cd, cr75, nalpha=50)

    # estimate cdmax
    AR = 1.0 / cr75
    cdmaxAR = 1.11 + 0.018*AR
    cdmax = max(maximum(cd), cdmaxAR)

    # find clmax
    i_ps = argmax(cl)  # positive stall
    cl_ps = cl[i_ps]
    cd_ps = cd[i_ps]
    a_ps = alpha[i_ps]

    # and clmin
    i_bs = alpha .< a_ps  # before stall
    i_ns = argmin(cl[i_bs])  # negative stall
    cl_ns = cl[i_bs][i_ns]
    cd_ns = cd[i_bs][i_ns]
    a_ns = alpha[i_bs][i_ns]

    # coefficients in method
    B1pos = cdmax
    A1pos = B1pos/2.0 * ones(nalpha)
    sa = sin(a_ps)
    ca = cos(a_ps)
    A2pos = (cl_ps - cdmax*sa*ca)*sa/ca^2
    B2pos = (cd_ps - cdmax*sa^2)/ca * ones(nalpha)

    B1neg = cdmax
    A1neg = B1neg/2.0
    sa = sin(a_ns)
    ca = cos(a_ns)
    A2neg = (cl_ns - cdmax*sa*ca)*sa/ca^2 * ones(nalpha)
    B2neg = (cd_ns - cdmax*sa^2)/ca * ones(nalpha)

    # angles of attack to extrapolate to
    apos = range(alpha[end], pi, length=nalpha+1)
    apos = apos[2:end]  # don't duplicate point
    aneg = range(-pi, alpha[1], length=nalpha+1)
    aneg = aneg[1:end-1]  # don't duplicate point

    # high aoa adjustments
    adjpos = ones(nalpha)
    idx = findall(apos .>= pi/2)
    adjpos[idx] .= -0.7
    A1pos[idx] .*= -1
    B2pos[idx] .*= -1

    # idx = findall(aneg .<= -alpha[end])

    adjneg = ones(nalpha)
    idx = findall(aneg .<= -pi/2)
    adjneg[idx] .= 0.7
    A2neg[idx] .*= -1
    B2neg[idx] .*= -1

    # extrapolate
    clpos = @. adjpos * (A1pos*sin(2*apos) + A2pos*cos(apos)^2/sin(apos))
    cdpos = @. B1pos*sin(apos)^2 + B2pos*cos(apos)
    clneg = @. adjneg * (A1neg*sin(2*aneg) + A2neg*cos(aneg)^2/sin(aneg))
    cdneg = @. B1neg*sin(aneg)^2 + B2neg*cos(aneg)

    # # override region between -alpha_high and alpha_low (if it exists)
    # idx = findall(-alpha[end] .<= aneg .<= alpha[1])
    # @. clneg[idx] = -cl[end]*0.7 + (aneg[idx]+alpha[end])/(alpha[1]+alpha[end])*(cl[1]+cl[end]*0.7)
    # @. cdneg[idx] = cd[1] + (aneg[idx]-alpha[1])/(-alpha[end]-alpha[1])*(cd[end]-cd[1])

    # override with linear variation at ends
    idx = findall(apos .>= pi-a_ps)
    @. clpos[idx] = (apos[idx] - pi)/a_ps*cl_ps*0.7
    idx = findall(aneg .<= -pi-a_ns)
    @. clneg[idx] = (aneg[idx] + pi)/a_ns*cl_ns*0.7

    # concatenate
    alphafull = [aneg; alpha; apos]
    clfull = [clneg; cl; clpos]
    cdfull = [cdneg; cd; cdpos]

    # don't allow negative drag
    cdfull = max.(cdfull, 0.0001)
    return alphafull, clfull, cdfull
end


text:
module TestHighlevel

using Test
using ONNXRunTime
const OX = ONNXRunTime
using ONNXRunTime: juliatype

@testset "high level" begin
    @testset "increment2x3.onnx" begin
        path = OX.testdatapath("increment2x3.onnx")
        model = OX.load_inference(path, execution_provider=:cpu)
        @test OX.input_names(model) == ["input"]
        @test OX.output_names(model) == ["output"]
        input = randn(Float32, 2,3)
        #= this works             =# model(Dict("input" => randn(Float32, 2,3)), ["output"])
        @test_throws OX.ArgumentError model(Dict("nonsense" => input), ["output"])
        @test_throws OX.ArgumentError model(Dict("input" => input), ["nonsense"])
        @test_throws OX.OrtException  model(Dict("input" => input), String[])
        @test_throws OX.ArgumentError model(Dict("input" => input, "unused"=>input), ["output"])
        @test_throws OX.ArgumentError model(Dict("input" => input, "unused"=>input), ["output"])
        @test_throws OX.OrtException model(Dict("input" => randn(Float32, 3,2)), ["output"])
        @test_throws Exception       model(Dict("input" => randn(Int, 2,3)    ), ["output"])
        @test_throws OX.OrtException model(Dict("input" => randn(Float64, 2,3)), ["output"])
        y = model(Dict("input" => input), ["output"])
        @test y == Dict("output" => input .+ 1f0)
        y = model(Dict("input" => input))
        @test y == Dict("output" => input .+ 1f0)
    end
    @testset "adder1x2x3.onnx" begin
        path = OX.testdatapath("adder1x2x3.onnx")
        model = OX.load_inference(path)
        @test OX.input_names(model) == ["x", "y"]
        @test OX.output_names(model) == ["sum"]
        x = randn(Float32, 1,2,3)
        y = randn(Float32, 1,2,3)
        d = model(Dict("x" => x, "y"=>y))
        @test d isa AbstractDict
        @test d == Dict("sum" => x+y)
    end
    @testset "diagonal1x2x3x4.onnx" begin
        path = OX.testdatapath("diagonal1x2x3x4.onnx")
        model = OX.load_inference(path)
        @test OX.input_names(model) == ["in"]
        @test OX.output_names(model) == ["out1", "out2"]
        x = randn(Float64, 1,2,3,4)
        d = model(Dict("in" => x))
        @test d == Dict("out1" => x, "out2" => x)
    end
    @testset "swap_x_.onnx" begin
        path = OX.testdatapath("swap_x_.onnx")
        model = OX.load_inference(path)
        @test OX.input_names(model)  == ["in1", "in2"]
        @test OX.output_names(model) == ["out1", "out2"]
        in1 = randn(Float32, 2,3)
        in2 = randn(Float32, 4,5)
        res = model((;in1, in2))
        @test keys(res) === (:out1, :out2)
        @test res isa NamedTuple
        @test res.out1 == in2
        @test res.out2 == in1
    end
    @testset "getindex_12.onnx" begin
        path = OX.testdatapath("getindex_12.onnx")
        model = OX.load_inference(path)
        inputs = (input=collect(reshape(1f0:20, 4,5)),)
        out = model(inputs).output
        @test inputs.input[2,3] == only(out)
    end
    @testset "copy2d.onnx" begin
        path = OX.testdatapath("copy2d.onnx")
        model = OX.load_inference(path)
        inputs = (input=randn(Float32,3,4),)
        out = model(inputs).output
        @test inputs.input == out
    end
    @testset "matmul.onnx" begin
        path = OX.testdatapath("matmul.onnx")
        model = OX.load_inference(path)
        inputs = (
                  input1 = randn(Float32, 2,3),
                  input2 = randn(Float32, 3,4),
        )
        out = model(inputs).output
        @test out ≈ inputs.input1 * inputs.input2
    end
    @testset "xyz_3x4x5.onnx" begin
        path = OX.testdatapath("xyz_3x4x5.onnx")
        model = OX.load_inference(path)
        inputs = (input=randn(Float32,4,10),)
        out = model(inputs)
        @test out.identity == inputs.input
        @test size(out.X) == size(out.Y) == size(out.Z) == (3,4,5)
        for _ in 1:20
            ix = rand(axes(out.X,1))
            iy = rand(axes(out.X,1))
            iz = rand(axes(out.X,1))
            @test out.X[ix,iy,iz] == ix-1
            @test out.Y[ix,iy,iz] == iy-1
            @test out.Z[ix,iy,iz] == iz-1
        end
    end
end

end#module

ext: jl


text:
using PoissonFE
using Test

@testset "PoissonFE.jl" begin
    y = [1.0, 0.0, 0.0, 0.0, 7.0, 1.0, 0.0, 1.0, 0.0, 0.0, 6.0,
    2.0, 3.0, 0.0, 1.0, 6.0, 0.0, 7.0, 0.0, 21.0]
    id = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
    X = [0.398106 0.496961; -0.612026 -0.224875; 0.34112 -1.11714;
    -1.12936 -0.394995; 1.43302 1.54983; 1.9804 -0.743514; -0.367221 -2.33171;
    -1.04413 0.812245; 0.56972 -0.501311; -0.135055 -0.510887;
    2.40162 -1.21536; -0.03924 -0.0225586; 0.689739 0.701239;
    0.0280022 -0.587482; -0.743273 -0.606728; 0.188792 1.09664;
    -1.80496 -0.24751; 1.46555 -0.159902; 0.153253 -0.625778;
    2.17261 0.900435]
    m = PoissonFEModel(y, X, id)
    @test isa(m, PoissonFEModel)
    coefs, se = fit(m)
    @test isapprox(coefs, [0.9244995462032737, 0.8693710731458504]; rtol=0.00001)
    @test isapprox(se, [0.013507707199824406, 0.12757995666403202]; rtol=0.00001)
end


text:
eval(Base,:(have_color = true))
using ASTInterpreter
using Base.Test

# Properly handle `call`
function CallTest()
    UnitRange{Int64}(2,2)
end

interp = enter(CallTest, Environment())
@test ASTInterpreter.finish!(interp; print_step = false, recursive=true) == CallTest()

# Properly handle :meta annotations
function MetaTest()
    @Base._pure_meta
    0
end

interp = enter(MetaTest, Environment())
@test ASTInterpreter.finish!(interp; print_step = false, recursive=true) == MetaTest()

# Test Vararg handling
function VATest(x...)
    x
end
callVA() = VATest()

interp = enter(callVA, Environment())
@test ASTInterpreter.finish!(interp; print_step = false, recursive=true) == callVA()

# Test Loops
function LoopTest()
    x = Int[]
    for i = 1:2
        push!(x, i)
    end
    x
end

interp = ASTInterpreter.enter_call_expr(nothing,:($(LoopTest)()))
@test ASTInterpreter.finish!(interp; print_step = false, recursive=true) == LoopTest()

# Test Loops
function ContinueTest()
    x = Int[]
    for i = 1:3
        if true
            push!(x, i)
            continue
        end
        error("Fell through")
    end
    x
end

interp = ASTInterpreter.enter_call_expr(nothing,:($(ContinueTest)()))
@test ASTInterpreter.finish!(interp; print_step = false, recursive=true) == ContinueTest()

#foo() = 1+1
function foo(n)
    x = n+1
    ((BigInt[1 1; 1 0])^x)[2,1]
end

interp = ASTInterpreter.enter_call_expr(nothing,:($(foo)(20)))
#@test ASTInterpreter.finish!(interp; print_step = false, recursive=true) == foo(20)


text:
using Sundials, Test
using DiffEqProblemLibrary
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

prob = prob_ode_linear
dt = 1 // 2^(4)
saveat = float(collect(0:dt:1))
sol = solve(prob, CVODE_BDF())
sol = solve(prob, CVODE_Adams())
@test sol.errors[:l2] < 1e-3
sol = solve(prob, CVODE_Adams(), reltol = 1e-5)
@test sol.errors[:l2] < 1e-5
sol = solve(prob, CVODE_Adams(), reltol = 1e-5, abstol = [1e-7])
sol = solve(prob, CVODE_Adams(), saveat = saveat)

@test sol.t == saveat

sol = solve(prob, CVODE_Adams(), saveat = dt)

@test sol.t == saveat

sol = solve(prob, CVODE_Adams(), saveat = saveat, save_everystep = true)

@test sol.t != saveat
@test intersect(sol.t, saveat) == saveat

sol = solve(prob, CVODE_Adams(), saveat = saveat, save_everystep = true, save_start = false)

@test sol.t[1] != 0

sol = solve(prob, CVODE_Adams(), tstops = [0.2, 0.5, 0.7])
@test all(t ∈ sol.t for t in [0.2, 0.5, 0.7])

prob = prob_ode_2Dlinear
sol = solve(prob, CVODE_BDF())
sol = solve(prob, CVODE_Adams())
sol = solve(prob, CVODE_Adams(), saveat = saveat)

@test sol.t == saveat

sol = solve(prob, CVODE_Adams(), saveat = [prob.tspan[2]])

@test sol.t == [prob.tspan[2]]

sol = solve(prob, CVODE_Adams(), saveat = saveat, save_everystep = false)

@test sol.t == saveat

sol = solve(prob, CVODE_Adams(), tstops = [0.9])

@test 0.9 ∈ sol.t

# Test the other function conversions
k = (du, u, p, t) -> du[1] = u[1]
prob = ODEProblem(k, [1.0], (0.0, 1.0))
sol = solve(prob, CVODE_BDF())
h = (u, p, t) -> u
u0 = [
    1.0 2.0
    3.0 2.0
]
prob = ODEProblem(h, u0, (0.0, 1.0))
sol = solve(prob, CVODE_BDF())

# Test Algorithm Choices
sol1 = solve(prob, CVODE_BDF(method = :Functional))
sol2 = solve(prob, CVODE_BDF(linear_solver = :Band, jac_upper = 3, jac_lower = 3))
sol3 = solve(prob, CVODE_BDF(linear_solver = :Diagonal))
sol4 = solve(prob, CVODE_BDF(linear_solver = :GMRES))
sol5 = solve(prob, CVODE_BDF(linear_solver = :FGMRES))
sol6 = solve(prob, CVODE_BDF(linear_solver = :PCG))
sol7 = solve(prob, CVODE_BDF(linear_solver = :BCG))
sol8 = solve(prob, CVODE_BDF(linear_solver = :TFQMR))
sol9 = solve(prob, CVODE_BDF(linear_solver = :Dense))
#sol9 = solve(prob,CVODE_BDF(linear_solver=:KLU)) # Requires Jacobian
sol10 = solve(prob, CVODE_BDF(linear_solver = :LapackDense))
sol11 = solve(prob, CVODE_BDF(linear_solver = :LapackBand, jac_upper = 3, jac_lower = 3))

@test isapprox(sol1[end], sol2[end], rtol = 1e-3)
@test isapprox(sol1[end], sol3[end], rtol = 1e-3)
@test isapprox(sol1[end], sol4[end], rtol = 1e-3)
@test isapprox(sol1[end], sol5[end], rtol = 1e-3)
@test isapprox(sol1[end], sol6[end], rtol = 1e-3)
@test isapprox(sol1[end], sol7[end], rtol = 1e-3)
@test isapprox(sol1[end], sol8[end], rtol = 1e-3)
#@test isapprox(sol1[end],sol9[end],rtol=1e-3)

# Test identity preconditioner
global prec_used = false
global psetup_used = false
prec = (z, r, p, t, y, fy, gamma, delta, lr) -> (global prec_used = true; z .= r)
psetup =
    (p, t, u, du, jok, jcurPtr, gamma) -> (global psetup_used = true; jcurPtr[] = false)
sol4 = solve(prob, CVODE_BDF(linear_solver = :GMRES, prec_side = 3, prec = prec))
@test isapprox(sol1[end], sol4[end], rtol = 1e-3)
@test prec_used
sol4 = solve(
    prob,
    CVODE_BDF(linear_solver = :GMRES, prec_side = 3, prec = prec, psetup = psetup),
)
@test isapprox(sol1[end], sol4[end], rtol = 1e-3)
@test psetup_used

# Backwards
prob = deepcopy(prob_ode_2Dlinear)
prob2 = ODEProblem(prob.f, prob.u0, (1.0, 0.0), 1.01)
sol = solve(prob2, CVODE_BDF())
@test maximum(diff(sol.t)) < 0 # Make sure all go negative

number_test(u, p, t) = -u^2 + (p[1] + t + p[2]) * u + p[2]
u0 = 0.0;
tspan = (0.0, 10)
prob = ODEProblem(number_test, u0, tspan, (2.0, 0.01))
sol = solve(prob, CVODE_BDF())

sol = solve(prob2, CVODE_BDF(), maxiters = 1)
@test length(sol.t) == 2


text:
#!/usr/bin/env julia

function part1(elems)
    [(e ÷ 3 - 2) for e in elems] |> sum
end

function helper(elem)
    elem = (elem ÷ 3) - 2
    ans = elem
    while ans >= 0
        ans = (ans ÷ 3) - 2
        if ans > 0
            elem += ans;
        end # if
    end
    elem
end # function

function part2(elems)
    [helper(e) for e in elems] |> sum
end

function main()
    if isempty(ARGS)
        println(stderr, "Usage:\n\tjulia main.jl <INPUT_FILE>\nOR\n\tmain.jl <INPUT_FILE>")
        exit(1)
    else
        elems = [parse(Int64, e) for e in read(ARGS[1], String) |> split]
    end

    println("Part 1: $(part1(elems))")
    println("Part 2: $(part2(elems))")
end

main()

max_stars_repo_name: GrayJack/Advent-Of-Code

max_issues_repo_name: GrayJack/Advent-Of-Code

max_forks_repo_name: GrayJack/Advent-Of-Code


text:
cond1 = make_condition_dict(num_nets=750,mcts_d=[8],net_type=:random,t_bnds=[0.7],seed=5,n_bnds=[8,120],role=:trusted)
cond2 = make_condition_dict(num_nets=250,mcts_d=[3],net_type=:random,t_bnds=[0.7],seed=4,n_bnds=[8,120],role=:candidate)
cond3 = make_condition_dict(num_nets=250,mcts_d=[1],net_type=:random,t_bnds=[0.7],seed=3,n_bnds=[8,120],role=:candidate)

exp_dict = Dict(:name=>"n_vary", :conditions=>Dict(1=>cond1, 2=>cond2, 3=> cond3))

fldr = "logs"
for x in keys(exp_dict[:conditions])
    fname = make_useful_fname(exp_dict,x)
    exp_dict[:conditions][x][:fname] = fname
    exp_dict[:conditions][x][:fldr] = fldr
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using DataFrames, Gadfly

set_default_plot_size(21cm, 8cm)

D = DataFrame(x=1:10, y=[0.11, 0.2, 0.1, 0.15, 0.3, 0.45, 0.37, 0.81, 0.62, 0.5])
palettef = Scale.lab_gradient("darkgreen","orange", "blue")
p1 = plot(D, x=:x, y=:y, color=:x, Geom.point,
    Scale.color_continuous(colormap=palettef, minvalue=0, maxvalue=10),
    Guide.title("Scale.color_continuous, Theme(alphas=[0.5])"),
    Theme(alphas=[0.5], continuous_highlight_color=identity,
        point_size=2mm)
)
p2 = plot(D, x=:x, y=:y, alpha=:x, Geom.point,
    Scale.alpha_continuous(minvalue=0, maxvalue=10),
    Guide.title("Scale.alpha_continuous, Theme(default_color=\"blue\")"),
    Theme(default_color="blue", discrete_highlight_color=c->"gray",
        point_size=2mm)
)
hstack(p1, p2)


text:
module OdsReaders

using PyCall
using Printf

export openods, readsheet, getsheets
export searchods, searchsheet

const py_pip = PyNULL()
const py_ods = PyNULL()

function __init__()
    u = pyimport("importlib.util")
    if u[:find_spec]("pyexcel_ods") == nothing
        py_pip = pyimport("pip._internal")
        py_pip[:main](["install","pyexcel-ods"])
        #@pyimport pip._internal as p
    end
    copy!(py_ods, pyimport("pyexcel_ods"))
end

"""
   openods function

   Open ods file and return all sheet data
   ods file is for OpenDocument SpreadSheet (e.g. LibreOffice Calc)
   if no file, return nothing

"""
function openods(filename::AbstractString)
    if !isfile(filename)
        error("$filename : file not found")
    else
        return py_ods.get_data(filename)
    end
end

"""
    readsheet function

    open ods file, set a target sheet name and args
    return sheet data for a target sheet

"""
function readsheet(filename::AbstractString, shname::AbstractString; args...)
    data = openods(filename)
    if haskey(data, shname)
        return conv_args(data[shname]; args...)
    else
        @printf "%s : sheet not found\n" shname
        return error()
    end
end

"""
    convert args function
    if header flag is false, first column of data is skipped
    if header flag is true, first column of data is not skipped

"""
function conv_args(sheet; header=true)
    if !isa(header, Bool)
        error("true or false is a valid argument for header")
    else
        if header
            return sheet
        else
            (r, c) = size(sheet)
            return sheet[2:r,:]
        end
    end
end

"""
    getsheet name function

    return sheet name lists of data

"""
function getsheets(filename::AbstractString)
    data = openods(filename)
    return [ k for (k,v) in data]
end

"""
    search string function in ods file (All Sheets)

    retrun search string dictionary of data
    key is sheet name, value is list of tuple(row, col)

"""
function searchods(filename::AbstractString, target)
    D = Dict() # Initilization
    for (k,v) in openods(filename)
        D = merge(D, searchdata(k, v, target))
    end
    D
end

"""
    search string function in a sheet

    retrun search string dictionary of data
    key is sheet name, value is list of tuple(row, col)
"""
function searchsheet(filename::AbstractString, shname::AbstractString, target)
   d = readsheet(filename, shname)
   searchdata(shname, d, target)
end

"""
    internal function of searching row and col value
"""
function searchdata(key, val, m)
    D = Dict()
    if length(size(val)) == 1
        # lack of data matrix
        for i in 1:length(val)
            o = findall(x->x==m, val[i])
            if length(o) != 0
                get!(D, key, [])
                append!(D[key], [ (i,j) for j in o ])
            end
        end
    else
        # full data matrix
        o = findall(x->x==m, val)
        if length(o) != 0
            D[key] = o
        end
    end
    D
end

end # module


text:
fname = :TRKsparse
c = Combi(hessian_sparse,PDataK,solve_modelKTR,preprocessKTR,decreaseKTR,TparamsKTR())
include("Template.jl")


text:
"""
    K, Cl = h2synthesize(P::ExtendedStateSpace, γ = nothing)

Synthesize H₂-optimal controller K and calculate the closed-loop transfer function from `w` to `z`.
Ref: Cha. 14.5 in Robust and Optimal Control.

If `γ = nothing`, use the formulas for H₂ in Ch 14.5. If γ is a large value, the H∞ formulas are used. As γ → ∞, these two are equivalent. The h∞ formulas do a coordinate transfromation that handles slightly more general systems so if you run into an error, it might be worth trying setting γ to something large, e.g., 1000.
"""
function h2synthesize(P::ExtendedStateSpace, γ = nothing)

    if γ === nothing
        X2, Y2, F2, L2 = _solvematrixequations2(P)
        Â2 = P.A + P.B2*F2 + L2*P.C2
        K = ss(Â2, -L2, F2, 0)
        return K, lft(ss(P), K)
    end

    P̄, Ltrans12, Rtrans12, Ltrans21, Rtrans21 = _transformp2pbar(P)
    X2, Y2, F2, L2 = _solvematrixequations(P̄, γ)

    _checkfeasibility(
        X2,
        Y2,
        γ,
        1e-3,
        10;
        verbose = false,
    ) || throw(DomainError("Not feasible for γ = $γ"))

    K = _synthesizecontroller(
        P̄,
        X2,
        Y2,
        F2,
        L2,
        γ,
        Ltrans12,
        Rtrans12,
        Ltrans21,
        Rtrans21,
    )
    Cl = lft(ss(P), K)
    K, Cl
end

function _solvematrixequations2(P::ExtendedStateSpace)

    A = P.A
    B1 = P.B1
    B2 = P.B2
    C1 = P.C1
    C2 = P.C2
    D11 = P.D11
    D12 = P.D12
    D21 = P.D21
    D22 = P.D22

    P1 = size(C1, 1)
    P2 = size(C2, 1)
    M1 = size(B1, 2)
    M2 = size(B2, 2)

    HX = [A zeros(size(A)); -C1'*C1 -A'] - [B2; -C1'*D12] * [D12'*C1 B2']
    HY = [A' zeros(size(A)); -B1*B1' -A] - [C2'; -B1*D21'] * [D21*B1' C2]

    # Solve matrix equations
    X2 = _solvehamiltonianare(HX)
    Y2 = _solvehamiltonianare(HY)

    F2 = -(B2'X2 + D12'C1)

    L2 = -(B1 * D21' + Y2 * C2')

    return X2, Y2, F2, L2
end


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash(show_undo_redo=true)

app.layout = html_div() do
 dcc_input(id="a"),
 html_div(id="b")
end

callback!(app,
    Output("b","children"),
    Input("a","value")
    ) do inputValue
    return inputValue
end

run_server(app)


text:
struct LibSixelDecoder <: AbstractSixelDecoder
    # (Experimental) internal fields

    # We need `allocator` to constructor libsixel objects, however, users of this
    # Julia package is not expected to use this field as it really should just live
    # in the C world.
    allocator::SixelAllocator

    function LibSixelDecoder(allocator=SixelAllocator())
        new(allocator)
    end
end

function (dec::LibSixelDecoder)(bytes::Vector{UInt8}; transpose=false)
    pixels = Ref{Ptr{Cuchar}}()
    palette = Ref{Ptr{Cuchar}}()
    pwidth = Ref(Cint(0))
    pheight = Ref(Cint(0))
    ncolors = Ref(Cint(0))

    status = C.sixel_decode_raw(
        bytes, length(bytes),
        pixels, pwidth, pheight,
        palette, ncolors,
        dec.allocator.ptr
    )
    check_status(status)

    palette_length = 256
    ncolors[] <= palette_length || throw(ArgumentError("palette `ncolors` should be less than $(palette_length), instead it is $(ncolors[])"))

    index = unsafe_wrap(Matrix{Cuchar}, pixels[], (pwidth[], pheight[]); own=false)

    # libsixel declares it to be ARGB but it's actually RGB{N0f8}
    PT = RGB{N0f8}
    pvalues = convert(Ptr{PT}, palette[])
    # TODO: we can acutally further compress the data with `values` of length `ncolors`
    #       Does it worth it? Maybe, but not very much, though.
    values = unsafe_wrap(Vector{PT}, pvalues, (palette_length, ); own=false)
    values = OffsetArray(values, OffsetArrays.Origin(0)) # again, libsixel assumes 0-based indexing

    return index, values
end


text:
_eltype(::Type{T}) where {T <: AbstractArray} = eltype(T)
_eltype(::Type{T}) where {T} = T
for (i, field) in enumerate((:vertextype, :facetype, :normaltype,
                       :texturecoordinatetype, :colortype))
    @eval begin
        $field(t::Type{T}) where {T <: HomogenousMesh} = _eltype(fieldtype(t, $i))
        $field(mesh::HomogenousMesh) = $field(typeof(mesh))
    end
end

hasvertices(msh) = vertextype(msh) != Nothing
hasfaces(msh) = facetype(msh) != Nothing
hasnormals(msh) = normaltype(msh) != Nothing
hastexturecoordinates(msh) = texturecoordinatetype(msh) != Nothing
hascolors(msh) = colortype(msh) != Nothing

vertices(msh) = msh.vertices
faces(msh) = msh.faces
normals(msh) = msh.normals
texturecoordinates(msh) = msh.texturecoordinates
colors(msh) = msh.color

# Bad, bad name! But it's a little tricky to filter out faces and verts from the attributes, after get_attribute
function attributes_noVF(m::T) where T <: AbstractMesh
    fielddict = Dict{Symbol, Any}(map(fieldnames(T)[3:end]) do field
        field => getfield(m, field)
    end)
    return filter(fielddict) do (key, val)
        val != nothing && val != Nothing[]
    end
end
#Gets all non Nothing attributes from a mesh in form of a Dict fieldname => value
function attributes(m::AbstractMesh)
    filter(((key,val),) -> (val != nothing && val != Nothing[]), all_attributes(m))
end
#Gets all non Nothing attributes types from a mesh type fieldname => ValueType
function attributes(m::Type{M}) where M <: HMesh
    filter(((key,val),) -> (val != Nothing && val != Vector{Nothing}), all_attributes(M))
end

function all_attributes(m::Type{M}) where M <: HMesh
    Dict{Symbol, Any}(map(field -> (field => fieldtype(M, field)), fieldnames(M)))
end
function all_attributes(m::M) where M <: HMesh
    Dict{Symbol, Any}(map(field -> (field => getfield(m, field)),  fieldnames(M)))
end

# Needed to not get into an stack overflow
convert(::Type{M}, mesh::Union{AbstractGeometry, AbstractMesh}) where {M <: AbstractMesh} = M(mesh)
convert(::Type{T}, mesh::T) where T <: AbstractMesh = mesh
# (::Type{HM1}){HM1 <: AbstractMesh}(mesh::HM1) = mesh

isvoid(::Type{T}) where {T} = false
isvoid(::Type{Nothing}) = true
isvoid(::Type{Vector{T}}) where {T} = isvoid(T)

@generated function (::Type{HM1})(primitive::HM2) where {HM1 <: HomogenousMesh, HM2 <: HomogenousMesh}
    fnames = fieldnames(HM1)
    expr = Expr(:call, HM1)
    for i in 1:fieldcount(HM1)
        field = fnames[i]
        target_type = fieldtype(HM1, i)
        source_type = fieldtype(HM2, i)
        if !isconcretetype(fieldtype(HM1, i))  # target is not defined
            push!(expr.args, :(getfield(primitive, $(QuoteNode(field)))))
        elseif !isvoid(target_type) && isvoid(source_type) # target not there yet, maybe we can decompose though (e.g. normals)
            push!(expr.args, :(decompose($(HM1.parameters[i]), primitive)))
        elseif isvoid(target_type)
            push!(expr.args, :($(target_type())))
        else
            push!(expr.args, :(convert($target_type, getfield(primitive, $(QuoteNode(field))))))
        end
    end
    expr
end

#Should be:
#function call{M <: HMesh, VT <: Point, FT <: Face}(::Type{M}, vertices::Vector{VT}, faces::Vector{FT})
# Haven't gotten around to implement the types correctly with abstract types in FixedSizeArrays
function (::Type{M})(
        vertices::AbstractVector{Point{3, VT}}, faces::AbstractVector{FT}
    ) where {M <: HMesh, VT, FT <: Face}
    msh = PlainMesh{VT, FT}(vertices = vertices, faces = faces)
    convert(M, msh)
end
get_default(x::Union{Type, TypeVar}) = nothing
get_default(x::Type{X}) where {X <: AbstractArray} = Nothing[]

"""
generic constructor for abstract HomogenousMesh, infering the types from the keywords (which have to match the field names)
some problems with the dispatch forced me to use this method name... need to further investigate this
"""
function homogenousmesh(attribs::Dict{Symbol, Any})
    newfields = map(fieldnames(HMesh)) do name
        get(attribs, name, get_default(fieldtype(HMesh, name)))
    end
    HomogenousMesh(newfields...)
end

"""
Creates a mesh from keyword arguments, which have to match the field types of the given concrete mesh
"""
(::Type{M})(; kw_args...) where {M <: HMesh} = M(Dict{Symbol, Any}(kw_args))

"""
Creates a new mesh from a dict of `fieldname => value` and converts the types to the given meshtype
"""
function (::Type{M})(attribs::Dict{Symbol, Any}) where M <: HMesh
    newfields = map(fieldnames(HomogenousMesh)) do field
        target_type = fieldtype(M, field)
        default = fieldtype(HomogenousMesh, field) <: AbstractVector ? Nothing[] : nothing
        get(attribs, field, default)
    end
    convert(M, HomogenousMesh(newfields...))
end

"""
Creates a new mesh from an old one, with changed attributes given by the keyword arguments
"""
function (::Type{M})(mesh::AbstractMesh, attributes::Dict{Symbol, Any}) where M <: HMesh
    newfields = map(fieldnames(HomogenousMesh)) do field
        get(attributes, field, getfield(mesh, field))
    end
    HomogenousMesh(newfields...)
end
"""
Creates a new mesh from an old one, with a new constant attribute (like a color)
"""
function (::Type{HM})(mesh::AbstractMesh, constattrib::ConstAttrib) where {HM <: HMesh, ConstAttrib}
    result = Dict{Symbol, Any}()
    for (field, target_type) in zip(fieldnames(HM), HM.parameters)
        if target_type <: ConstAttrib
            result[field] = constattrib
        elseif target_type != Nothing
            result[field] = decompose(target_type, mesh)
        end
    end
    HM(result)
end
function add_attribute(m::AbstractMesh, attribute)
    attribs = attributes(m) # get all attribute values as a Dict fieldname => value
    attribs[:color] = attribute # color will probably be renamed to attribute. not sure yet...
    homogenousmesh(attribs)
end

"""
Creates a new mesh from a tuple of a geometry type and a constant attribute
"""
function (::Type{HM})(x::Tuple{P, ConstAttrib}) where {HM <: HMesh, ConstAttrib, P<:AbstractGeometry}
    any, const_attribute = x
    add_attribute(HM(any), const_attribute)
end

merge(m::AbstractVector{M}) where {M <: AbstractMesh} = merge(m...)

"""
Merges an arbitrary mesh. This function probably doesn't work for all types of meshes
parameters:
`m1` first mesh
`meshes...` other meshes
"""
function merge(m1::M, meshes::M...) where M <: AbstractMesh
    v = copy(m1.vertices)
    f = copy(m1.faces)
    attribs = deepcopy(attributes_noVF(m1))
    for mesh in meshes
        append!(f, mesh.faces .+ length(v))
        append!(v, mesh.vertices)
        for (v1, v2) in zip(values(attribs), values(attributes_noVF(mesh)))
            append!(v1, v2)
        end
    end
    attribs[:vertices]  = v
    attribs[:faces]     = f
    return M(attribs)
end

"""
A mesh with one constant attribute can be merged as an attribute mesh.
Possible attributes are FSArrays
"""
function merge(
        m1::HMesh{_1, _2, _3, _4, ConstAttrib, _5, _6},
        meshes::HMesh{_1, _2, _3, _4, ConstAttrib, _5, _6}...
    ) where {_1, _2, _3, _4, ConstAttrib <: Colorant, _5, _6}
    vertices     = copy(m1.vertices)
    faces        = copy(m1.faces)
    attribs      = Dict{Symbol, Any}(filter(((k,v),) -> k != :color, attributes_noVF(m1)))
    attribs      = Dict{Symbol, Any}([(k, copy(v)) for (k,v) in attribs])
    color_attrib = RGBA{U8}[RGBA{U8}(m1.color)]
    index        = Float32[length(color_attrib)-1 for i=1:length(m1.vertices)]
    for mesh in meshes
        append!(faces, mesh.faces .+ length(vertices))
        append!(vertices, mesh.vertices)
        attribsb = attributes_noVF(mesh)
        for (k,v) in attribsb
            k != :color && append!(attribs[k], copy(v))
        end
        push!(color_attrib, mesh.color)
        append!(index, Float32[length(color_attrib)-1 for i=1:length(mesh.vertices)])
    end
    delete!(attribs, :color)
    attribs[:vertices]      = vertices
    attribs[:faces]         = faces
    attribs[:attributes]    = color_attrib
    attribs[:attribute_id]  = index
    return HMesh{_1, _2, _3, _4, Nothing, typeof(color_attrib), eltype(index)}(attribs)
end

# Fast but slightly ugly way to implement mesh multiplication
# This should probably go into FixedSizeArrays.jl, Vector{FSA} * FSA
struct MeshMulFunctor{T}
    matrix::Mat{4,4,T}
end
(m::MeshMulFunctor{T})(vert) where {T} = Vec{3, T}(m.matrix*Vec{4, T}(vert..., 1))
function *(m::Mat{4,4,T}, mesh::AbstractMesh) where T
    msh = deepcopy(mesh)
    map!(MeshMulFunctor(m), msh.vertices, msh.vertices)
    msh
end


text:
# Scaling for numerical stability - may not be needed
scale(x) = x / 10000

get_reference(S::TensorMeta, node::XNode) = S.reference_map[node]
graph(S::TensorMeta) = S.graph
Profiler.users(S::TensorMeta) = S.users

# Accessor methods
rb(I::ILPHolder) = I.read_bandwidth
wb(I::ILPHolder) = I.write_bandwidth
rba(I::ILPHolder) = I.read_bandwidth_async
wba(I::ILPHolder) = I.write_bandwidth_async

# Common Methods

# Length check because ngraph compilation is not 100 % consistent and can sometimes have
# a few more nodes than it began with ...
limit(F::Frame, args...) = limit(F.modeltype, args...)
limit(S::ILPHolder, i) = i > length(S.dram_limits) ? maxlimit(S) : S.dram_limits[i]
maxlimit(S::ILPHolder) = maximum(S.dram_limits)

predict(F::Frame) = objective_value(F.model)
descriptor(F::Frame, tensor::XTensor) = F.modeltype.descriptors[tensor]

#####
##### Entry Point
#####

# For procedurally building JuMP expressions for the ILP model
const AffExprType = typeof(AffExpr())

function create_model(modeltype::ILPHolder, profile_data::FunctionData)
    preprocess!(modeltype, profile_data)

    # Start with an empty model that we will progressively build.
    model = Model(AutoTM.getoptimizer())
    frame = Frame(modeltype, model, profile_data)

    # Going deep into JuMP here - the idea is to build the objective as a bunch of aff exprs
    # and eventually combine all of them together.
    model[:node_times] = Dict{String, AffExprType}()
    model[:tensor_async] = Dict{String, AffExprType}()
    model[:tensor_sync] = AffExprType()

    add_tensors!(frame)
    add_nodes!(frame)
    add_constraints!(frame)

    # Default objective is to just sum all of the node times
    objective_expr = model[:tensor_sync]
    for (node_name, node_times) in model[:node_times]
        # Check to see if there are overlapping async transfers.
        #
        # If so, take the max of the sum of the overlapping transfers and the node time.
        _async = get(model[:tensor_async], node_name, nothing)
        if isnothing(_async)
            add_to_expression!(objective_expr, node_times)
        else
            var = @variable(model, lower_bound = 0.0)
            @constraint(model, var >= node_times)
            @constraint(model, var >= _async)
            add_to_expression!(objective_expr, var)
        end
    end
    # Quick optimization to remove zero terms
    drop_zeros!(objective_expr)
    @objective(frame.model, Min, objective_expr)

    return frame
end

#####
##### Adding Tensors
#####

isfixed(frame::Frame, tensor) = descriptor(frame, tensor).isfixed
#isfixed(frame::Frame, tensor) = false
free_tensors(frame::Frame) = [t for t in tensors(frame.profile_data) if !isfixed(frame, t)]
fixed_tensors(frame::Frame) = [t for t in tensors(frame.profile_data) if isfixed(frame, t)]

function add_tensors!(frame::Frame)
    data = frame.profile_data
    modeltype = frame.modeltype

    # Create variables for the tensors and add flow constraints to the to the tensor graphs
    @variable(frame.model,
        tensor_graphs[
            tensor = free_tensors(frame),
            e = edges(graph(descriptor(frame, tensor)))
        ],
        Bin
    )

    @showprogress 1 "Creating Flow Formulation " for tensor in free_tensors(frame)
        g = graph(descriptor(frame, tensor))
        # Iterate through nodes in the graph - generating constraints based on the type
        # of node.
        for v in vertices(g)
            # Set flow coming out of the source node
            if getmeta(g, v).location == LOC_SOURCE
                @constraint(frame.model,
                    sum(tensor_graphs[tensor, e] for e in outedges(g, v)) == 1
                )

            # Set flow going into the sink node
            elseif getmeta(g, v).location == LOC_SINK
                @constraint(frame.model,
                    sum(tensor_graphs[tensor, e] for e in inedges(g, v)) == 1
                )

            # All other ops must conserve flow
            else
                oe = collect(outedges(g, v))
                ie = collect(inedges(g, v))
                @constraint(frame.model,
                    sum(tensor_graphs[tensor, e] for e in oe) -
                    sum(tensor_graphs[tensor, e] for e in ie) == 0
                )
            end
        end
    end

    #####
    ##### Add objective penalty for moving data
    #####

    add_movement_formulations!(frame)

    #####
    ##### Create variables to determine if a tensor is in DRAM.
    #####

    @variable(frame.model,
        tensor_in_dram[
            tensor = free_tensors(frame),
            user = nGraph.name.(users(descriptor(frame, tensor)))
        ],
        Bin
    )

    @variable(frame.model,
        tensor_in_dram_post[
            tensor = free_tensors(frame),
            user = nGraph.name.(users(descriptor(frame, tensor)))
        ],
        Bin
    )

    # A tensor in DRAM is live if any of its incoming edges are used.
    for tensor in free_tensors(frame)
        desc = descriptor(frame, tensor)
        g = graph(desc)

        # Create a container for the critical edge from LOC_SOURCE to LOC_SINK if this is
        # an argument tensor
        skip_edge = []
        if isarg(tensor)
            # Find the source to sink edge.
            # Will error if this edge doesn't exist, so serves as a form of error checking
            edge = find_edge(
                (g, e) -> getmeta(g, src(e)).location == LOC_SOURCE &&
                          getmeta(g, dst(e)).location == LOC_SINK,
                g
            )
            push!(skip_edge, edge)
        end

        for user in users(desc)
            # Get the DRAM for this op.
            verts = filter(
                v -> isdram(getmeta(g,v).location) && getmeta(g,v).op == user,
                vertices(g)
            )

            # Map `inedges` to `vertex_iter` and iterats over all those edges
            _iter = vflatten(Iterators.flatten(inedges.(Ref(g), verts)), skip_edge)
            for e in _iter
                @constraint(
                    frame.model,
                    tensor_in_dram[tensor, nGraph.name(user)] >= tensor_graphs[tensor, e]
                )
            end

            # If all incoming edges are not taken, tensor MUST not be in DRAM.
            @constraint(frame.model,
                sum(tensor_graphs[tensor, e] for e in _iter) >=
                    tensor_in_dram[tensor, nGraph.name(user)]
            )

            # Similary, set the post DRAM constraints
            _edges = filter(
                e -> (isdram(getmeta(g, src(e)).location) &&
                    # Need to check "LOC_SINK" for the static case
                    (
                        isdram(getmeta(g, dst(e)).location) ||
                        getmeta(g, dst(e)).location == LOC_SINK
                    ) && getmeta(g, src(e)).op == user),
                collect(edges(g))
            )
            _iter = vflatten(_edges, skip_edge)
            for e in _iter
                @constraint(
                    frame.model,
                    tensor_in_dram_post[tensor, nGraph.name(user)] >= tensor_graphs[tensor, e]
                )
            end

            @constraint(frame.model,
                sum(tensor_graphs[tensor, e] for e in _iter) >=
                    tensor_in_dram_post[tensor, nGraph.name(user)]
            )
        end
    end
    return nothing
end

# Filter on edge type, sort by parent index to get the edges in execution order.
_find_edges(g, edgetype) = sort(
    filter(e -> getmeta(g, e).edgetype == edgetype, collect(edges(g))),
    by = src
)

# Formulation specific move node stuff
add_movement_formulations!(frame::Frame{ILPHolder{Static}}) = nothing
function add_movement_formulations!(frame::Frame)
    # Unpack variables
    data = frame.profile_data
    modeltype = frame.modeltype

    tensor_sync_expr = frame.model[:tensor_sync]
    tensor_async_dict = frame.model[:tensor_async]
    tensor_graphs = frame.model[:tensor_graphs]

    # A tensor is written to dram if:
    # - It was not created into PMEM
    # - Any edge from DRAM to PMEM is taken
    #
    # NOTE: We only pay the write cost once.
    @variable(frame.model, tensor_write[tensor = free_tensors(frame)], Bin)

    # Add objective terms for all read ops
    for tensor in free_tensors(frame)
        # Some unpacking
        g = graph(descriptor(frame, tensor))
        bytes = sizeof(tensor)

        # Take the ceiling of all these to ensure there's always a cost to moving.
        read_cost        = scale(bytes / rb(modeltype))
        write_cost       = scale(bytes / wb(modeltype))
        read_cost_async  = scale(bytes / rba(modeltype))
        write_cost_async = scale(bytes / wba(modeltype))

        # Collect Edges according to type.
        sync_reads = _find_edges(g, EDGE_SYNC_READ)
        sync_writes = _find_edges(g, EDGE_SYNC_WRITE)
        async_reads = _find_edges(g, EDGE_ASYNC_READ)
        async_writes = _find_edges(g, EDGE_ASYNC_WRITE)

        #####
        ##### Constraints for sync write variable
        #####

        # No Sync write if any async write
        for e in sync_writes
            @constraint(
                frame.model,
                tensor_write[tensor] >= tensor_graphs[tensor, e]
            )
        end

        # No sync write if no edges taken
        @constraint(
            frame.model,
            tensor_write[tensor] <= sum(tensor_graphs[tensor, e] for e in sync_writes)
        )

        #####
        ##### Constraints on async write variables
        #####

        # Assign each edge to the kernel it overlaps with.
        #
        # Just go overkill and grab all the edges, even though we end up only using a subset.
        kernels = Dict(e => getmeta(g,src(e)).op for e in edges(g))

        # Create read variables expressions
        for e in async_reads
            _expr = get!(tensor_async_dict, nGraph.name(kernels[e]), AffExprType())
            move_var = tensor_graphs[tensor, e]
            add_to_expression!(_expr, read_cost_async, move_var)
            dict_push!(frame.modeltype.async_move_vars, kernels[e], move_var)
        end

        # Create write variables.
        for e in async_writes
            _expr = get!(tensor_async_dict, nGraph.name(kernels[e]), AffExprType())
            move_var = tensor_graphs[tensor, e]
            add_to_expression!(_expr, write_cost_async, move_var)
            dict_push!(frame.modeltype.async_move_vars, kernels[e], move_var)
        end

        #####
        ##### Finally, add all the synchonous move costs.
        #####
        for e in sync_reads
            add_to_expression!(tensor_sync_expr, read_cost, tensor_graphs[tensor, e])
        end
        add_to_expression!(tensor_sync_expr, write_cost, tensor_write[tensor])
    end
    return nothing
end

# There's an issue when trying to reference whether or not a tensor is in DRAM.
#
# If we're on an op where the tensor is used, we have to look at the inputs to a
# graph verted with LOC_DRAM or LOC_PREAD to see if the tensor was fetched or already
# lived in dram.
#
# If we're on an op where a tensor is LIVE but not READ, we need to check the outgoing
# edge of the correct DRAM -> DRAM node to see if the tensor just lives around in DRAM.
function get_tensor_in_dram(F::Frame, tensor::XTensor, node::XNode)
    # First - check if the node is fixed. If it is fixed, then we can returna constant
    # depending on its location.
    if isfixed(F, tensor)
        location = collect(locations(tensor))
        if length(location) != 1
            @show location
            @show tensor
            @show node
            error()
        end
        if first(location) == DRAM
            return 1
        else
            return 0
        end
    end

    # Otherwise, look through the generated model to find the variable that will indicate
    # if this tensor is in DRAM or not.
    desc = descriptor(F, tensor)
    if in(node, users(desc))
        return F.model[:tensor_in_dram][tensor, nGraph.name(node)]
    else
        return F.model[:tensor_in_dram_post][tensor, nGraph.name(get_reference(desc, node))]
    end
end

function add_nodes!(F::Frame)
    data = F.profile_data

    # Create decision variables for all nodes that have a choice of backend algorithm.
    select_nodes = filter(x -> hasprofile(x) && can_select_algo(x), nodes(data))
    if !isempty(select_nodes)
        @info "Creating Algorithms Variables"
        @variable(
            F.model,
            algo_var[
                node = select_nodes,
                enum = enums(gettime(node))
            ],
            Bin
        )

        # Constrain so only one algorithm may be selected.
        for node in select_nodes
            @constraint(
                F.model,
                sum(algo_var[node, e] for e in enums(gettime(node))) == 1
            )
        end
    end

    for node in nodes(data)
        # We don't profile all ops, so perform a quick check to see if this is an op
        # the we have profile information for. If not, there's nothing to do as far as the
        # ILP model is concerned.
        hasprofile(node) || continue

        # The GPU path of this code will just return an all DRAM config - which will be
        # useful for generating the constraint that all kernel IO for the GPU case must
        # reside in GPU DRAM.
        #
        # The CPU path will yield a bunch of DRAM/PMEM combinations
        configs = configs_for(node)

        # Create a variable for each config.
        vars = @variable(F.model, [config = configs], Bin)

        for config in configs
            # Create an expression for the input and output locations
            expr = AffExprType()
            iter = Iterators.flatten((
                zip(config.inputs, inputs(node)),
                zip(config.outputs, outputs(node))
            ))

            for (location, tensor) in iter
                # use `jump_tensor` because it's really a JuMP variable that is returned
                # by this call.
                jump_tensor = get_tensor_in_dram(F, tensor, node)
                if location == DRAM
                    # This branch handles the jump_tensor being either an Int or a
                    # VariableRef just fine.
                    add_to_expression!(expr, jump_tensor)
                    @constraint(F.model, vars[config] <= jump_tensor)
                else
                    # The return value from `get_tensor_in_dram` can be a literal value,
                    # so make sure we handle this case.
                    if isa(jump_tensor, Int )
                        add_to_expression!(expr, jump_tensor)
                    else
                        add_to_expression!(expr, 1)
                        add_to_expression!(expr, -1, jump_tensor)
                    end
                    @constraint(F.model, vars[config] <= 1 - jump_tensor)
                end
            end

            @constraint(
                F.model,
                vars[config] + length(config.inputs) + length(config.outputs) >= 1 + expr
            )
        end

        # Add a valid contraint to help the solver
        @constraint(F.model, sum(vars[config] for config in configs) == 1)

        # Create an expression for this node's expected running time.
        node_times = AffExprType()
        for config in configs
            # If we can select the algorithm for this node, we need to generate some more
            # variables to "AND" the config with the algorithm selection to ensure that
            # we only get a single algorithm out at the end.
            #
            # If there are not multiple algorithms, then we don't have to worry about it.
            if can_select_algo(node, config)
                v = @variable(F.model, [enum = enums(gettime(node, config))], Bin)
                for enum in enums(gettime(node, config))
                    @constraint(F.model, v[enum] <= algo_var[node, enum])
                    @constraint(F.model, v[enum] <= vars[config])
                    @constraint(F.model, v[enum] + 1 >= vars[config] + algo_var[node, enum])

                    coeff = scale(Profiler.timeat(gettime(node, config), enum))
                    add_to_expression!(node_times, coeff, v[enum])
                end
            else
                coeff = scale(gettime(node, config))
                add_to_expression!(node_times, coeff, vars[config])
            end
        end
        F.model[:node_times][nGraph.name(node)] = node_times
    end
    return nothing
end

# Allocations in ngraph happen on 4096 bytes boundaries. For better accuracty, round
# up to the nearest multiple of 4096 before figuring out the number of bytes.
#
# Take the floor to introduce more zeros into the ILP formulation. This shouldn't really
# make much of a difference.
tensor_size(t::XTensor) = tensor_size(sizeof(t))
tensor_size(sz) = ceil(Int, ceil(Int, sz / 4096) * 4096 / 1E6)

function add_constraints!(F::Frame)
    # Unpack some variables
    data = F.profile_data

    for (index, tensors) in enumerate(live_tensors(data))
        node = nodes(data)[index]
        hasprofile(node) || continue
        F.modeltype.node_to_limit_index[nGraph.name(node)] = index

        # Add DRAM constraint for the workspace
        if can_select_algo(node)
            v = F.model[:algo_var]
            algo_expr = @expression(
                F.model,
                sum(
                    tensor_size(Profiler.bytesat(gettime(node), e)) *
                    v[node, e] for e in enums(gettime(node))
                )
            )
        else
            # If we can't select the algorithm, just create an empty expression that will
            # be hapily optimized away when we create the size constraint.
            algo_expr = AffExprType()
        end

        if !isempty(tensors)
            @constraint(F.model,
                algo_expr + sum(tensor_size(t) * get_tensor_in_dram(F, t, node)
                    for t in tensors
                    if !iszero(tensor_size(t))) <= limit(F, index)
            )
        end
    end

    return nothing
end

#####
##### Queries
#####

function find_shortcut_edge(frame, tensor)
    g = descriptor(frame, tensor).graph
    edge = find_edge(
        (g, e) -> getmeta(g, src(e)).location == LOC_SOURCE && getmeta(g, dst(e)).location == LOC_SINK,
        g
    )
    return edge
end

# Check if a given tensor is:
# 1. A function argument (nGraph function input or output)
# 2. Assigned to DRAM for the lifetime of the function
function islocalarg(frame, tensor::XTensor)
    # Determine if it's an argument
    isarg(tensor) || return false

    # Next - figure out if it was fixed in a memory pool.
    if isfixed(frame, tensor)
        location = first(locations(tensor))
        return location == DRAM
    end

    # Find the critical edge from LOC_SOURCE to LOC_SINK
    # This tensor is assigned to the local memory if this edge is taken (has a solved value
    # of one)
    edge = find_shortcut_edge(frame, tensor)

    return approx_one(frame.model[:tensor_graphs][tensor, edge])
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


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

using Markdown
using InteractiveUtils

# ╔═╡ 4d38c306-5ced-11ec-32c8-354a64d62d72
using Pkg

# ╔═╡ feec0e1a-0b9c-4e18-9131-bdaa4ba8981d
Pkg.develop(path="../.")

# ╔═╡ 94a14d53-eda6-4840-a389-51c3df60f017
using Plots, FourierTools, NDTools

# ╔═╡ 80c0e146-0ddf-42f5-a4e2-b528b8e73b40
x = 2π .* fftpos(1,5,3)

# ╔═╡ 75810e0b-442f-46d1-a484-f78d638cf576
f(x) = sin(x)

# ╔═╡ 9903e739-58eb-403e-a903-94a3b4815ab9
collect(fftpos(1, 5, 3))

# ╔═╡ 97dd506f-46f8-46a8-bbd8-a6bd6b58ef79
collect(fftpos(1, 101, 51))

# ╔═╡ 4735a86c-0ecc-45d9-91c7-26eed0f51804
collect(fftpos(1, 100, 51))

# ╔═╡ e0ecb758-1825-46c8-9d72-2baa39c1eaf6
md"
The even interpolation somehow extrapolates (!) only the left tail. the right tail
"

# ╔═╡ 82b62bf7-587a-4787-81ae-46fc328bf87a
plot(fftpos(1, 20, 11), circshift(resample(circshift(f.(x), -2), 20), 10))

# ╔═╡ 7475a27b-4423-49db-9b69-0ec15c7c2a89
plot(fftpos(1, 21, 11), circshift(resample(circshift(f.(x), -2), 21), 10))

# ╔═╡ Cell order:
# ╠═4d38c306-5ced-11ec-32c8-354a64d62d72
# ╠═feec0e1a-0b9c-4e18-9131-bdaa4ba8981d
# ╠═94a14d53-eda6-4840-a389-51c3df60f017
# ╠═80c0e146-0ddf-42f5-a4e2-b528b8e73b40
# ╠═75810e0b-442f-46d1-a484-f78d638cf576
# ╟─9903e739-58eb-403e-a903-94a3b4815ab9
# ╠═97dd506f-46f8-46a8-bbd8-a6bd6b58ef79
# ╠═4735a86c-0ecc-45d9-91c7-26eed0f51804
# ╟─e0ecb758-1825-46c8-9d72-2baa39c1eaf6
# ╠═82b62bf7-587a-4787-81ae-46fc328bf87a
# ╠═7475a27b-4423-49db-9b69-0ec15c7c2a89


text:
using Test
using DynamicPolynomials
using MultivariatePolynomials
using LazySets: Hyperrectangle, low, high
using CarlemanLinearization: generate_monomials, build_matrix, quadratic_matrix_form, kron_pow, lift_vector

@testset "Kronecker power (symbolic)" begin
    @polyvar x[1:2]
    y = kron_pow(x, 2)
    @test findall(==(x[1]*x[2]), y) == [2, 3]
    @test findall(==(x[1]^2), y) == [1]
    @test findall(==(x[2]^2), y) == [4]
    @test findall(==(x[2]^3), y) == Int[]
end

@testset "Conversion from polynomial to matrix representation" begin
    vars = @polyvar x y
    dx = 3x + y^2
    dy = x - y - 2.2 * x * y + x^2
    f = [dx, dy]
    F1, F2 = quadratic_matrix_form(f, vars)
    @test F1 == [3.0 0; 1 -1]
    @test F2 == [0.0 0 0 1; 1 -2.2 0 0]
end

@testset "Generating all commutative monomials" begin
    for (n, D) in [(3, 0), (3, 1), (3, 2), (5, 3), (5, 10)]
        monomials = generate_monomials(n, D)
        # No reprtitions
        @test length(monomials) == length(Set(monomials))
        # Correct size
        @test length(monomials) == binomial(n + D, n)
        # Nonnegativity
        @test all([sum(m) == sum(map(abs, m)) for m in monomials])
        # ordering
        @test all([sum(monomials[i - 1]) <= sum(monomials[i]) for i in 2:length(monomials)])
        # not exceeding D
        @test sum(monomials[end]) == D
    end
end

@testset "Compressed matrices" begin
    F1 = [3.0 0; 1 -1]
    F2 = [0.0 0 0 1; 1 -2.2 0 0]
    result = build_matrix(F1, F2, 1; compress = true)
    @test result == F1

    result = build_matrix(F1, F2, 2; compress = true)
    @test result == [
        3. 0. 0. 0. 1.;
        1. -1. 1. -2.2 0;
        0 0 6. 0 0;
        0 0 1. -1. 0;
        0 0 0 2. -2.
    ]
end

@testset "Lifting vectors" begin
    x0 = Hyperrectangle(low=[0, 1, -1], high=[1, 2, 3])
    N = 2
    lifted = lift_vector(x0, N)
    println(zip(low(lifted), high(lifted)))
    sides = Set(zip(low(lifted), high(lifted)))
    @test sides == Set{Tuple{Rational, Rational}}([(0, 1), (1, 2), (-1, 3), (0, 1), (1, 4), (0, 9), (0, 2), (-1, 3), (-2, 6)])
end


text:
const eidx = 30

```
    ui_advanced(img::AstroImage, indx::Int = 1)

Opens an advanced GUI panel to visualize AstroImages.
```
function ui_advanced(img::AstroImage, indx::Int = 1)              # 3 X eidx grid
    win = Window("Image", 700,700, border_width = 5)
    g = Grid()
    set_gtk_property!(g, :column_homogeneous, true)
    set_gtk_property!(g, :column_spacing, 15)
    push!(win, g)
    img_labels = Array{Tuple{Tuple{Real,Real},String},1}(undef,0)

    # Sets toolbar widget
    toolbar = Toolbar()
    headerbtn = get_button("Header")
    closebtn = get_button("Close")
    labelbtn = get_button("Show Labels")
    savebtn = get_button("Save")
    resetbtn = get_button("Reset")
    map(w->push!(toolbar,w),[headerbtn, labelbtn, savebtn,resetbtn, closebtn])
    g[1:3,1] = toolbar

    # Create canvas to display image
    c = canvas(UserUnit)
    g[1:2,2:eidx-3] = c
    set_gtk_property!(c, :expand, true)

    if(length(img) > 1)
        g[3,6] = get_label("AstroImage index-> 1:"* string.(length(img)))
        g[3,7] = GtkEntry()
        g[3,8] = GtkButton("Change Index")
    end

    # Image Label widgets
    g[3,9] = get_label("Pixel Coordinate X")
    g[3,10] = GtkEntry()
    g[3,11]= get_label("Pixel Coordinate Y")
    g[3,12]= GtkEntry()
    g[3,13]= get_label("Label")
    g[3,14]= GtkEntry()
    g[3,15]= GtkButton("Add Label")

    # Sets labels showing informations
    g[1,eidx-2] = get_label("Physical Coordinate X =")
    g[2,eidx-2] = get_label("Physical Coordinate Y =")
    g[1,eidx-1] = get_label("World Coordinate X =")
    g[2,eidx-1] = get_label("World Coordinate Y =")
    g[1:2,eidx] = get_label("Value = ")
    g[3,eidx] = get_label("Tip: Ctrl+Scroll to zoom")

    fv,zr,surf = form_image(img, indx)

    # interactive signals
    init_zoom_scroll(c,zr)                     # Makes canvas scrollable
    close_window(closebtn, win)                # Handles close button
    draw_canvas(c,zr,surf)                     # Draws canvas on window

    # Mapping mouse motion
    map(c.mouse.motion) do btn
        xu, yu = btn.position.x.val, btn.position.y.val
        world_coord = pix_to_world(img.wcs[indx], [xu, yu])
        GAccessor.text(g[1,eidx-2] ,"Pixel Coordinate X = $(round(xu, digits = 2))")
        GAccessor.text(g[2,eidx-2] ,"Pixel Coordinate Y = $(round(yu, digits = 2))")
        GAccessor.text(g[1,eidx-1] ,"World Coordinate X = $(round(world_coord[1], digits = 2))"*"°")
        GAccessor.text(g[2,eidx-1] ,"World Coordinate Y = $(round(world_coord[2], digits = 2))"*"°")
        if xu >= 1 && yu >= 1
            GAccessor.text(g[1,eidx] ,"Value = $(img.data[indx][Int(round(xu)),Int(round(yu))])")
        else
            GAccessor.text(g[1,eidx] , "Value = ")
        end
    end

    # Handles label adding
    signal_connect(g[3,15], :clicked) do btn
        x = parse(Float64, get_gtk_property(g[3,10], :text, String))
        y = parse(Float64, get_gtk_property(g[3,12], :text, String))
        label = get_gtk_property(g[3,14], :text, String)
        push!(img_labels, ((x,y), label))
        set_gtk_property!(g[3,10], :text, "")
        set_gtk_property!(g[3,12], :text, "")
        set_gtk_property!(g[3,14], :text, "")
    end

    # Handles index changing
    if(length(img) > 1)
        signal_connect(g[3,8], :clicked) do btn
            tmp = parse(Int, get_gtk_property(g[3,7], :text, String))
            if tmp <= length(img)
                indx = tmp
                set_gtk_property!(g[3,7], :text, "")
                fv,zr,surf = form_image(img, indx)
                draw_canvas(c,zr,surf)
                map(c.mouse.motion) do btn
                    xu, yu = btn.position.x.val, btn.position.y.val
                    world_coord = pix_to_world(img.wcs[indx], [xu, yu])
                    GAccessor.text(g[1,eidx-2] ,"Pixel Coordinate X = $(round(xu, digits = 2))")
                    GAccessor.text(g[2,eidx-2] ,"Pixel Coordinate Y = $(round(yu, digits = 2))")
                    GAccessor.text(g[1,eidx-1] ,"World Coordinate X = $(round(world_coord[1], digits = 2))"*"°")
                    GAccessor.text(g[2,eidx-1] ,"World Coordinate Y = $(round(world_coord[2], digits = 2))"*"°")
                    if xu >= 1 && yu >= 1
                        GAccessor.text(g[1,eidx] ,"Value = $(img.data[indx][Int(round(xu)),Int(round(yu))])")
                    else
                        GAccessor.text(g[1,eidx] , "Value = ")
                    end
                end
                init_zoom_scroll(c,zr)
            else
                set_gtk_property!(g[3,7], :text, "Set index within bound")
            end
        end
    end

    # Handles header showing
    signal_connect(headerbtn, :clicked) do widget
        showall(header_window(img.wcs[indx]))
    end

    # Handles Label showing
    signal_connect(labelbtn, :clicked) do widget
        showall(label_window(img_labels))
    end

    # Handles reset button
    signal_connect(resetbtn, :clicked) do widget
        img_labels = Array{Tuple{Tuple{Real,Real},String},1}(undef,0)
    end

    showall(win)
end


text:
# The paper:
#
#  Yi Sun, F. Gomez, T. Schaul, J. Schmidhuber. A Linear Time Natural Evolution
#  Strategy for Non-Separable Functions. In Proceedings of the Genetic and
#  Evolutionary Computation Conference (GECCO), Amsterdam, 2013.
#
# proposes a variant of a Natural Evolution Strategy which can successfully
# solve very high-dimensional problems. In particular they show state-of-the-art
# performance on the 512-dimensional Rosenbrock function.
#
# Here we try our DE implementations on this problem to see what kind of
# performance we get.

using BlackBoxOptim

function fitness_for_opt(problemname, numDimensions, populationSize, numSteps, method)

  p = BlackBoxOptim.example_problems[problemname]
  problem = instantiate(p, numDimensions)

  println("\n$(problem.name), n = $(numdims(problem)), optimizer = $(string(method))")

  result = bboptimize(problem; Method = method,
    NumDimensions = numDimensions, PopulationSize = populationSize,
    MaxSteps = numSteps
  )

  best_fitness(result)
end

problemname = "Rosenbrock"

# Call once to compile everything
fitness_for_opt(problemname, 2, 10, 100, :adaptive_de_rand_1_bin)

for (dims, steps) in [(128, 1e6), (256, 2e6), (512, 4e6)]
  println("Num dims = ", dims, ", Num steps = ", steps)
  fitness = fitness_for_opt(problemname, dims, 50, steps, :adaptive_de_rand_1_bin)
end