text:
#=
Author: Andrew Ning

A general blade element momentum (BEM) method for propellers/fans and turbines.

Some unique features:
- a simple yet very robust solution method
- allows for non-ideal conditions like reversed flow, or free rotation
- allows arbitrary inflow
- convenience methods for common wind turbine inflow scenarios

=#

module CCBlade

import FLOWMath

export Rotor, Section, OperatingPoint, Outputs
export af_from_files, af_from_data
export simple_op, windturbine_op
export solve, thrusttorque, nondim

# --------- structs -------------

"""
    Rotor(Rhub, Tip, B, turbine, pitch, precone)

Scalar parameters defining the rotor.

**Arguments**
- `Rhub::Float64`: hub radius (along blade length)
- `Rtip::Float64`: tip radius (along blade length)
- `B::Int64`: number of blades
- `turbine::Bool`: true if turbine, false if propeller
- `pitch::Float64`: pitch angle (rad).  defined same direction as twist.
- `precone::Float64`: precone angle
"""
mutable struct Rotor{TF, TI, TB}

    Rhub::TF
    Rtip::TF
    B::TI
    turbine::TB
    pitch::TF
    precone::TF

end

# convenience constructor for no precone and no pitch
Rotor(Rhub, Rtip, B, turbine) = Rotor(Rhub, Rtip, B, turbine, zero(Rhub), zero(Rhub))
Rotor(Rhub, Rtip, B, turbine, pitch) = Rotor(Rhub, Rtip, B, turbine, pitch, zero(Rhub))

"""
    Section(r, chord, theta, af)

Define sectional properties for one station along rotor

**Arguments**
- `r::Float64`: radial location along blade (`Rhub < r < Rtip`)
- `chord::Float64`: corresponding local chord length
- `theta::Float64`: corresponding twist angle (radians)
- `af::function`: a function of the form: `cl, cd = af(alpha, Re, Mach)`
"""
mutable struct Section{TF, TAF}

    r::TF
    chord::TF
    theta::TF
    af::TAF

end

# make rotor broadcastable as a single entity
Base.Broadcast.broadcastable(r::Rotor) = Ref(r)

# convenience function to set fields within an array of structs
function Base.setproperty!(obj::Array{Rotor{TF, TI, TB}, N}, sym::Symbol, x) where {TF, TI, TB, N}
    setfield!.(obj, sym, x)
end

# convenience function to access fields within an array of structs
function Base.getproperty(obj::Array{Section{TF, TAF}, N}, sym::Symbol) where {TF, TAF, N}
    return getfield.(obj, sym)
end

# convenience function to set fields within an array of structs
function Base.setproperty!(obj::Array{Section{TF, TAF}, N}, sym::Symbol, x) where {TF, TAF, N}
    setfield!.(obj, sym, x)
end

"""
    OperatingPoint(Vx, Vy, rho, mu=1.0, asound=1.0)

Operation point for a rotor.
The x direction is the axial direction, and y direction is the tangential direction in the rotor plane.
See Documentation for more detail on coordinate systems.
Vx and Vy vary radially at same locations as `r` in the rotor definition.

**Arguments**
- `Vx::Float64`: velocity in x-direction along blade
- `Vy::Float64`: velocity in y-direction along blade
- `rho::Float64`: fluid density
- `mu::Float64`: fluid dynamic viscosity (unused if Re not included in airfoil data)
- `asound::Float64`: fluid speed of sound (unused if Mach not included in airfoil data)
"""
mutable struct OperatingPoint{TF, TF2}
    Vx::TF
    Vy::TF
    rho::TF2  # different type to accomodate ReverseDiff
    mu::TF2
    asound::TF2
end

# convenience constructor when Re and Mach are not used.
OperatingPoint(Vx, Vy, rho) = OperatingPoint(Vx, Vy, rho, one(rho), one(rho))

# convenience function to access fields within an array of structs
function Base.getproperty(obj::Array{OperatingPoint{TF, TF2}, N}, sym::Symbol) where {TF, TF2, N}
    return getfield.(obj, sym)
end

# convenience function to set fields within an array of structs
function Base.setproperty!(obj::Array{OperatingPoint{TF, TF2}, N}, sym::Symbol, x) where {TF, TF2, N}
    setfield!.(obj, sym, x)
end

"""
    Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)

Outputs from the BEM solver along the radius.

**Arguments**
- `Np::Vector{Float64}`: normal force per unit length
- `Tp::Vector{Float64}`: tangential force per unit length
- `a::Vector{Float64}`: axial induction factor
- `ap::Vector{Float64}`: tangential induction factor
- `u::Vector{Float64}`: axial induced velocity
- `v::Vector{Float64}`: tangential induced velocity
- `phi::Vector{Float64}`: inflow angle
- `alpha::Vector{Float64}`: angle of attack
- `W::Vector{Float64}`: inflow velocity
- `cl::Vector{Float64}`: lift coefficient
- `cd::Vector{Float64}`: drag coefficient
- `cn::Vector{Float64}`: normal force coefficient
- `ct::Vector{Float64}`: tangential force coefficient
- `F::Vector{Float64}`: hub/tip loss correction
- `G::Vector{Float64}`: effective hub/tip loss correction for induced velocities: `u = Vx * a * G, v = Vy * ap * G`
"""
struct Outputs{TF}

    Np::TF
    Tp::TF
    a::TF
    ap::TF
    u::TF
    v::TF
    phi::TF
    alpha::TF
    W::TF
    cl::TF
    cd::TF
    cn::TF
    ct::TF
    F::TF
    G::TF

end

# convenience constructor to initialize
Outputs() = Outputs(0.0, 0.0, 0.0, 0.0, 0.1, 0.1, -5.0*pi/180.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

# convenience function to access fields within an array of structs
function Base.getproperty(obj::Array{Outputs{TF}, N}, sym::Symbol) where {TF, N}
    return getfield.(obj, sym)
end

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

# ----------- airfoil ---------------

"""
    parse_af_file(filename, header=1)

Private function. Read an airfoil file.  For one Reynolds/Mach number.
Additional data like cm is optional but will be ignored.
alpha should be in degrees

format:

header\n
alpha1 cl1 cd1 ...\n
alpha2 cl2 cd2\n
alpha3 cl3 cd3\n
...

Returns arrays for alpha (in radians), cl, cd
"""
function parse_af_file(filename, header=1)

    alpha = Float64[]
    cl = Float64[]
    cd = Float64[]

    open(filename) do f

        # skip header
        for i = 1:header
            readline(f)
        end

        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

    return alpha*pi/180, cl, cd

end

"""
    af_from_files(filenames; Re=[], Mach=[])

Read airfoil file(s) and return a function of the form `cl, cd = func(alpha, Re, M)`

If filenames is just one file, then Re and Mach are ignored (just aoa variation).
`af_file_files("somefile.dat")`

If filenames is a vector then there is variation with either Re or Mach (both not both).
`af_file_files(["f1.dat", "f2.dat", "f3.dat"], Mach=[0.5, 0.6, 0.7])`

Filenames can be a matrix for variation with both.
`af_file_files(filematrix, Re=[3e6, 5e6], Mach=[0.5, 0.6, 0.7])`
where `filematrix[i, j]` correspnds to `Re[i]`, `Mach[j]`

Uses the `af_from_data` function.
"""
function af_from_files(filenames; Re=[], Mach=[])

    # evalutae first file to determine size
    if isa(filenames, Array)
        sz = size(filenames)
    else
        sz = 0
    end
    nRe = length(Re)
    nMach = length(Mach)

    if sz == 0  # just one file
        alpha, cl, cd = parse_af_file(filenames)

    elseif length(sz) == 1  # list of files
        alpha, cl1, cd1 = parse_af_file(filenames[1])
        nalpha = length(alpha)
        ncond = length(filenames)

        cl = zeros(nalpha, ncond)
        cd = zeros(nalpha, ncond)
        cl[:, 1] = cl1
        cd[:, 1] = cd1

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

    else  # matrix of filenames
        alpha, cl1, cd1 = parse_af_file(filenames[1, 1])
        nalpha = length(alpha)

        cl = zeros(nalpha, nRe, nMach)
        cd = zeros(nalpha, nRe, nMach)

        for j = 1:nMach
            for i = 1:nRe
                _, clij, cdij = parse_af_file(filenames[i, j])
                cl[:, i, j] = clij
                cd[:, i, j] = cdij
            end
        end

    end

    return af_from_data(alpha, Re, Mach, cl, cd)
end

"""private function
1d interpolation vs alpha
"""
function afalpha(alpha, Re, Mach, cl, cd)

    nRe = length(Re)
    nMach = length(Mach)

    # squeeze out singleton dimensions if necessary
    if nRe == 1 && nMach == 1  # could be cases with 1, 0, but this would be inconcistent input and I'm not going to bother handlingn it.
        cl = cl[:, 1, 1]
        cd = cd[:, 1, 1]
    elseif nRe == 0 && nMach == 0
        cl = cl[:, 1]
        cd = cd[:, 1]
    end

    afcl = FLOWMath.Akima(alpha, cl)
    afcd = FLOWMath.Akima(alpha, cd)

    afeval(alpha_pt, Re_pt, M_pt) = afcl(alpha_pt)[1],
                                    afcd(alpha_pt)[1]
    return afeval
end

"""private function
2d interpolation vs alpha, Re
"""
function afalphaRe(alpha, Re, Mach, cl, cd)

    nRe = length(Re)
    nMach = length(Mach)

    # squeeze out singleton dimensions if necessary
    if nMach == 1
        cl = cl[:, :, 1]
        cd = cd[:, :, 1]
    end

    afeval(alpha_pt, Re_pt, M_pt) = FLOWMath.interp2d(FLOWMath.akima, alpha, Re, cl, alpha_pt, Re_pt)[1],
                                    FLOWMath.interp2d(FLOWMath.akima, alpha, Re, cd, alpha_pt, Re_pt)[1]
    return afeval

end

"""private function
2d interpolation vs alpha, Mach
"""
function afalphaMach(alpha, Re, Mach, cl, cd)

    nRe = length(Re)
    nMach = length(Mach)

    # squeeze out singleton dimensions if necessary
    if nRe == 1
        cl = cl[:, 1, :]
        cd = cd[:, 1, :]
    end

    afeval(alpha_pt, Re_pt, M_pt) = FLOWMath.interp2d(FLOWMath.akima, alpha, Mach, cl, alpha_pt, M_pt)[1],
                                    FLOWMath.interp2d(FLOWMath.akima, alpha, Mach, cd, alpha_pt, M_pt)[1]
    return afeval
end

"""private function
3d interpolation vs alpha, Re, Mach
"""
function afalphaReMach(alpha, Re, Mach, cl, cd)

    nRe = length(Re)
    nMach = length(Mach)

    afeval(alpha_pt, Re_pt, M_pt) = FLOWMath.interp3d(FLOWMath.akima, alpha, Re, Mach, cl, alpha_pt, Re_pt, M_pt)[1],
                                    FLOWMath.interp3d(FLOWMath.akima, alpha, Re, Mach, cd, alpha_pt, Re_pt, M_pt)[1]
    return afeval
end

"""
Create an airfoil function directly from alpha, cl, and cd arrays.
The function of the form `cl, cd = func(alpha, Re, M)`
alpha should be in radians.  Uses an akima spline.  `af_from_files` calls this function.

`cl[i, j, k]` corresponds to `alpha[i]`, `Re[j]`, `Mach[k]`

If `Mach=[]`
`cl[i, j]` corresponds to `alpha[i]`, `Re[j]`
`size(cl) = (length(alpha), length(Re))`
But you can use a singleton dimension for the constant Mach if desired.
`size(cl) = (length(alpha), length(Re), 1)`
The above also applies for `Re=[]` where variation is with alpha and Mach.

There is also a convenience method for vector data with just aoa variation
`af_from_data(alpha, cl, cd)` which just corresponds to `af_from_data(alpha, Re=[], Mach=[], cl, cd)`
"""
function af_from_data(alpha, Re, Mach, cl, cd)

    nRe = length(Re)
    nMach = length(Mach)

    if nRe <= 1 && nMach <= 1
        return afalpha(alpha, Re, Mach, cl, cd)
    elseif nMach <= 1
        return afalphaRe(alpha, Re, Mach, cl, cd)
    elseif nRe <= 1
        return afalphaMach(alpha, Re, Mach, cl, cd)
    else
        return afalphaReMach(alpha, Re, Mach, cl, cd)
    end
end

# convenience wrappers
af_from_data(alpha, cl, cd) = af_from_data(alpha, Re=[], Mach=[], cl, cd)

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

# ------------ BEM core ------------------

"""
(private) residual function
"""
function residual(phi, rotor, section, op)

    # unpack inputs
    r = section.r
    chord = section.chord
    theta = section.theta
    af = section.af
    Rhub = rotor.Rhub
    Rtip = rotor.Rtip
    B = rotor.B
    turbine = rotor.turbine
    pitch = rotor.pitch
    Vx = op.Vx
    Vy = op.Vy
    rho = op.rho

    # constants
    sigma_p = B*chord/(2.0*pi*r)
    sphi = sin(phi)
    cphi = cos(phi)

    # angle of attack
    alpha = phi - (theta + pitch)

    # Reynolds number
    W0 = sqrt(Vx^2 + Vy^2)  # ignoring induction, which is generally a very minor difference and only affects Reynolds/Mach number
    Re = rho * W0 * chord / op.mu

    # Mach number
    Mach = W0/op.asound  # also ignoring induction

    # airfoil cl/cd
    if turbine
        cl, cd = af(alpha, Re, Mach)
    else
        cl, cd = af(-alpha, Re, Mach)
        cl *= -1
    end

    # resolve into normal and tangential forces
    cn = cl*cphi + cd*sphi
    ct = cl*sphi - cd*cphi

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

    # sec parameters
    k = cn*sigma_p/(4.0*F*sphi*sphi)
    kp = ct*sigma_p/(4.0*F*sphi*cphi)

    # # parameters used in Vx=0 and Vy=0 cases
    k0 = cn*sigma_p/(4.0*F*sphi*cphi)
    k0p = ct*sigma_p/(4.0*F*sphi*sphi)

    # --- solve for induced velocities ------
    if isapprox(Vx, 0.0, atol=1e-6)

        u = sign(phi)*k0*Vy
        v = 0.0
        a = 0.0
        ap = 0.0
        R = sin(phi)^2 + sign(phi)*cn*sigma_p/(4.0*F)

    elseif isapprox(Vy, 0.0, atol=1e-6)

        u = 0.0
        v = k0p*abs(Vx)
        a = 0.0
        ap = 0.0
        R = sign(Vx)*4*F*sphi*cphi - ct*sigma_p

    else

        if phi < 0
            k *= -1
        end

        if isapprox(k, -1.0, atol=1e-6)  # state corresopnds to Vx=0, return any nonzero residual
            return 1.0
        end

        if k <= 2.0/3  # momentum region
            a = k/(1 + k)

        else  # empirical region
            g1 = 2.0*F*k - (10.0/9-F)
            g2 = 2.0*F*k - (4.0/3-F)*F
            g3 = 2.0*F*k - (25.0/9-2*F)

            if isapprox(g3, 0.0, atol=1e-6)  # avoid singularity
                a = 1.0 - 1.0/(2.0*sqrt(g2))
            else
                a = (g1 - sqrt(g2)) / g3
            end
        end

        u = a * Vx

        # -------- tangential induction ----------
        if Vx < 0
            kp *= -1
        end

        if isapprox(kp, 1.0, atol=1e-6)  # state corresopnds to Vy=0, return any nonzero residual
            return 1.0
        end

        ap = kp/(1 - kp)
        v = ap * Vy

        # ------- residual function -------------
        R = sin(phi)/(1 - a) - Vx/Vy*cos(phi)/(1 + ap)
    end

    # ------- loads ---------
    W = sqrt((Vx - u)^2 + (Vy + v)^2)
    Np = cn*0.5*rho*W^2*chord
    Tp = ct*0.5*rho*W^2*chord

    # The BEM methodology applies hub/tip losses to the loads rather than to the velocities.
    # This is the most common way to implement a BEM, but it means that the raw velocities are misleading
    # as they do not contain any hub/tip loss corrections.
    # To fix this we compute the effective hub/tip losses that would produce the same thrust/torque.
    # In other words:
    # CT = 4 a (1 - a) F = 4 a G (1 - a G)\n
    # This is solved for G, then multiplied against the wake velocities.

    if isapprox(Vx, 0.0, atol=1e-6)
        G = sqrt(F)
    elseif isapprox(Vy, 0.0, atol=1e-6)
        G = F
    else
        G = (1.0 - sqrt(1.0 - 4*a*(1.0 - a)*F))/(2*a)
    end
    u *= G
    v *= G

    if turbine
        return R, Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)
    else
        return R, Outputs(-Np, -Tp, -a, -ap, -u, -v, phi, -alpha, W, -cl, cd, -cn, -ct, F, G)
    end

end

"""
(private) Find a bracket for the root closest to xmin by subdividing
interval (xmin, xmax) into n intervals.

Returns found, xl, xu.
If found = true a bracket was found between (xl, xu)
"""
function firstbracket(f, xmin, xmax, n, backwardsearch=false)

    xvec = range(xmin, xmax, length=n)
    if backwardsearch  # start from xmax and work backwards
        xvec = reverse(xvec)
    end

    fprev = f(xvec[1])
    for i = 2:n
        fnext = f(xvec[i])
        if fprev*fnext < 0  # bracket found
            if backwardsearch
                return true, xvec[i], xvec[i-1]
            else
                return true, xvec[i-1], xvec[i]
            end
        end
        fprev = fnext
    end

    return false, 0.0, 0.0

end

"""
    solve(rotor, section, op)

Solve the BEM equations for given rotor geometry and operating point.

**Arguments**
- `rotor::Rotor`: rotor properties
- `section::Section`: section properties
- `op::OperatingPoint`: operating point

**Returns**
- `outputs::Outputs`: BEM output data including loads, induction factors, etc.
"""
function solve(rotor, section, op)

    # error handling
    if typeof(section) <: Vector
        error("You passed in an vector for section, but this funciton does not accept an vector.\nProbably you intended to use broadcasting (notice the dot): solve.(Ref(rotor), sections, ops)")
    end

    # parameters
    npts = 20  # number of discretization points to find bracket in residual solve

    # unpack
    Vx = op.Vx
    Vy = op.Vy
    theta = section.theta + rotor.pitch

    # ---- determine quadrants based on case -----
    Vx_is_zero = isapprox(Vx, 0.0, atol=1e-6)
    Vy_is_zero = isapprox(Vy, 0.0, atol=1e-6)

    # quadrants
    epsilon = 1e-6
    q1 = [epsilon, pi/2]
    q2 = [-pi/2, -epsilon]
    q3 = [pi/2, pi-epsilon]
    q4 = [-pi+epsilon, -pi/2]

    if Vx_is_zero && Vy_is_zero
        return Outputs()

    elseif Vx_is_zero

        startfrom90 = false  # start bracket search from 90 deg instead of 0 deg.

        if Vy > 0 && theta > 0
            order = (q1, q2)
        elseif Vy > 0 && theta < 0
            order = (q2, q1)
        elseif Vy < 0 && theta > 0
            order = (q3, q4)
        else  # Vy < 0 && theta < 0
            order = (q4, q3)
        end

    elseif Vy_is_zero

        startfrom90 = true  # start bracket search from 90 deg

        if Vx > 0 && abs(theta) < pi/2
            order = (q1, q3)
        elseif Vx < 0 && abs(theta) < pi/2
            order = (q2, q4)
        elseif Vx > 0 && abs(theta) > pi/2
            order = (q3, q1)
        else  # Vx < 0 && abs(theta) > pi/2
            order = (q4, q2)
        end

    else  # normal case

    # for i = 1:nr

        startfrom90 = false

        if Vx > 0 && Vy > 0
            order = (q1, q2, q3, q4)
        elseif Vx < 0 && Vy > 0
            order = (q2, q1, q4, q3)
        elseif Vx > 0 && Vy < 0
            order = (q3, q4, q1, q2)
        else  # Vx[i] < 0 && Vy[i] < 0
            order = (q4, q3, q2, q1)
        end

    end

    # ----- solve residual function ------

    # # wrapper to residual function to accomodate format required by fzero
    R(phi) = residual(phi, rotor, section, op)[1]

    success = false
    for j = 1:length(order)  # quadrant orders.  In most cases it should find root in first quadrant searched.
        phimin, phimax = order[j]

        # check to see if it would be faster to reverse the bracket search direction
        backwardsearch = false
        if !startfrom90
            if phimin == -pi/2 || phimax == -pi/2  # q2 or q4
                backwardsearch = true
            end
        else
            if phimax == pi/2  # q1
                backwardsearch = true
            end
        end

        # force to dual numbers if necessary
        phimin = phimin*one(section.chord)
        phimax = phimax*one(section.chord)

        # find bracket
        success, phiL, phiU = firstbracket(R, phimin, phimax, npts, backwardsearch)

        # once bracket is found, solve root finding problem and compute loads
        if success

            phistar = FLOWMath.brent(R, phiL, phiU)
            _, outputs = residual(phistar, rotor, section, op)
            return outputs
        end
    end

    # it shouldn't get to this point.  if it does it means no solution was found
    # it will return empty outputs
    # alternatively, one could increase npts and try again

    return Outputs()
end

# ------------ inflow ------------------

"""
    simple_op(Vinf, Omega, r, rho, mu=1.0, asound=1.0, precone=0.0)

Uniform inflow through rotor.  Returns an Inflow object.

**Arguments**
- `Vinf::Float`: freestream speed (m/s)
- `Omega::Float`: rotation speed (rad/s)
- `r::Float`: radial location where inflow is computed (m)
- `rho::Float`: air density (kg/m^3)
- `mu::Float`: air viscosity (Pa * s)
- `asounnd::Float`: air speed of sound (m/s)
- `precone::Float`: precone angle (rad)
"""
function simple_op(Vinf, Omega, r, rho; mu=one(rho), asound=one(rho), precone=zero(Vinf))
    # TODO: change this to keyword args in #master

    # error handling
    if typeof(r) <: Vector
        error("You passed in an vector for r, but this function does not accept an vector.\nProbably you intended to use broadcasting")
    end

    Vx = Vinf * cos(precone)
    Vy = Omega * r * cos(precone)

    return OperatingPoint(Vx, Vy, rho, mu, asound)

end

"""
    windturbine_op(Vhub, Omega, pitch, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho, mu=1.0, asound=1.0)

Compute relative wind velocity components along blade accounting for inflow conditions
and orientation of turbine.  See Documentation for angle definitions.

**Arguments**
- `Vhub::Float64`: freestream speed at hub (m/s)
- `Omega::Float64`: rotation speed (rad/s)
- `r::Float64`: radial location where inflow is computed (m)
- `precone::Float64`: precone angle (rad)
- `yaw::Float64`: yaw angle (rad)
- `tilt::Float64`: tilt angle (rad)
- `azimuth::Float64`: azimuth angle to evaluate at (rad)
- `hubHt::Float64`: hub height (m) - used for shear
- `shearExp::Float64`: power law shear exponent
- `rho::Float64`: air density (kg/m^3)
- `mu::Float64`: air viscosity (Pa * s)
- `asound::Float64`: air speed of sound (m/s)
"""
function windturbine_op(Vhub, Omega, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho, mu=1.0, asound=1.0)

    sy = sin(yaw)
    cy = cos(yaw)
    st = sin(tilt)
    ct = cos(tilt)
    sa = sin(azimuth)
    ca = cos(azimuth)
    sc = sin(precone)
    cc = cos(precone)

    # coordinate in azimuthal coordinate system
    x_az = -r*sin(precone)
    z_az = r*cos(precone)
    y_az = 0.0  # could omit (the more general case allows for presweep so this is nonzero)

    # get section heights in wind-aligned coordinate system
    heightFromHub = (y_az*sa + z_az*ca)*ct - x_az*st

    # velocity with shear
    V = Vhub*(1 + heightFromHub/hubHt)^shearExp

    # transform wind to blade c.s.
    Vwind_x = V * ((cy*st*ca + sy*sa)*sc + cy*ct*cc)
    Vwind_y = V * (cy*st*sa - sy*ca)

    # wind from rotation to blade c.s.
    Vrot_x = -Omega*y_az*sc
    Vrot_y = Omega*z_az

    # total velocity
    Vx = Vwind_x + Vrot_x
    Vy = Vwind_y + Vrot_y

    # operating point
    return OperatingPoint(Vx, Vy, rho, mu, asound)

end

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

# -------- convenience methods ------------

"""
    thrusttorque(rotor, sections, outputs::Vector{Outputs{TF}}) where TF

integrate the thrust/torque across the blade,
including 0 loads at hub/tip, using a trapezoidal rule.

**Arguments**
- `rotor::Rotor`: rotor object
- `sections::Vector{Section}`: rotor object
- `outputs::Vector{Outputs}`: output data along blade

**Returns**
- `T::Float64`: thrust (along x-dir see Documentation).
- `Q::Float64`: torque (along x-dir see Documentation).
"""
# function thrusttorque(rotor, sections, outputs)
function thrusttorque(rotor, sections, outputs::Vector{Outputs{TF}}) where TF

    # add hub/tip for complete integration.  loads go to zero at hub/tip.
    rfull = [rotor.Rhub; sections.r; rotor.Rtip]
    Npfull = [0.0; outputs.Np; 0.0]
    Tpfull = [0.0; outputs.Tp; 0.0]

    # integrate Thrust and Torque (trapezoidal)
    thrust = Npfull*cos(rotor.precone)
    torque = Tpfull.*rfull*cos(rotor.precone)

    T = rotor.B * FLOWMath.trapz(rfull, thrust)
    Q = rotor.B * FLOWMath.trapz(rfull, torque)

    return T, Q
end

"""
    thrusttorque(rotor, sections, outputs::Array{Outputs{TF}, 2}) where TF

Integrate the thrust/torque across the blade given an array of output data.
Generally used for azimuthal averaging of thrust/torque.
`outputs[i, j]` corresponds to `sections[i], azimuth[j]`.  Integrates across azimuth
"""
function thrusttorque(rotor, sections, outputs::Matrix{Outputs{TF}}) where TF

    T = 0.0
    Q = 0.0
    nr, naz = size(outputs)

    for j = 1:naz
        Tsub, Qsub = thrusttorque(rotor, sections, outputs[:, j])
        T += Tsub / naz
        Q += Qsub / naz
    end

    return T, Q
end

"""
    nondim(T, Q, Vhub, Omega, rho, rotor)

Nondimensionalize the outputs.

**Arguments**
- `T::Float64`: thrust (N)
- `Q::Float64`: torque (N-m)
- `Vhub::Float64`: hub speed used in turbine normalization (m/s)
- `Omega::Float64`: rotation speed used in propeller normalization (rad/s)
- `rho::Float64`: air density (kg/m^3)
- `rotor::Rotor`: rotor object

**Returns**

if windturbine
- `CP::Float64`: power coefficient
- `CT::Float64`: thrust coefficient
- `CQ::Float64`: torque coefficient

if propeller
- `eff::Float64`: efficiency
- `CT::Float64`: thrust coefficient
- `CQ::Float64`: torque coefficient
"""
function nondim(T, Q, Vhub, Omega, rho, rotor)

    P = Q * Omega
    Rp = rotor.Rtip*cos(rotor.precone)

    if rotor.turbine  # wind turbine normalizations

        q = 0.5 * rho * Vhub^2
        A = pi * Rp^2

        CP = P / (q * A * Vhub)
        CT = T / (q * A)
        CQ = Q / (q * Rp * A)

        return CP, CT, CQ

    else  # propeller

        n = Omega/(2*pi)
        Dp = 2*Rp

        if T < 0
            eff = 0.0  # creating drag not thrust
        else
            eff = T*Vhub/P
        end
        CT = T / (rho * n^2 * Dp^4)
        CQ = Q / (rho * n^2 * Dp^5)

        return eff, CT, CQ

    # elseif rotortype == "helicopter"

    #     A = pi * Rp^2

    #     CT = T / (rho * A * (Omega*Rp)^2)
    #     CP = P / (rho * A * (Omega*Rp)^3)
    #     FM = CT^(3/2)/(sqrt(2)*CP)

    #     return FM, CT, CP
    end

end

end  # module


text:
include("../src/ekp.jl") # for the ekpobj
include("../src/truth.jl") #for the covariance structure
include("../src/transforms.jl") # for transform_prior_to_real, transform_real_to_prior
include("../src/units.jl") # for apply_units_transform

using Plots
using Plots.PlotMeasures
using ColorSchemes
using Distributions
using Random
using LinearAlgebra
using JLD2
using Statistics

using .EKP

"""
Plotting script for the Parameter Ensemble and the GP,

produces two plots,

Requires:
--------
- EKP object
- parameter grid defining the bounds of the plotting region
- gp prediction of mean on the parameter grid
- gp prediction of covariance on the parameter grid
- transformations and inverse transformations real->prior space (if required)
- units transform to add units (if required)
- truthobject (the covariances)
"""

######################################################################33
function main()

    @load "ekp.jld2" ekpobj
    @load "truth.jld2" truthobj

    @load "mean_mesh.jld2" mean_mesh
    @load "var_mesh.jld2" var_mesh
    @load "param1vals.jld2" param1vals
    @load "param2vals.jld2" param1vals

    #Plotting parameters
    min_ekp_it = 1
    max_ekp_it = 6

    #contour plot levels
    n_levels = 30

    yt = ekpobj.g_t #the inflated sample used in ekp

    # the truth
    @load truthdir*"truth.jld2" truthobj
    gcm_truth_cov = cov(truthobj.sample)
    obs_truth_cov = truthobj.cov
    inflation_cov = obs_truth_cov - 2*gcm_truth_cov

    u = ekpobj.u[min_ekp_it:max_ekp_it]#it x [enssize x param]
    u_flat = cat(u...,dims=1)#[(itxens) x param]
    ens_size = ekpobj.J

    # meshes, see emulate_sample.jl, "emulator_uncertainty_graphs(...)"
    # note we do not account for the additional terms from prior in this objective

    # transform y into SVD coordinates first.
    SVD = svd(gcm_truth_cov) #svd.U * svd.S * svd.Vt (can also get Vt)
    Dinv = Diagonal(1.0 ./ sqrt.(SVD.S)) #diagonal matrix of 1/eigenvalues
    y_transform = Dinv*SVD.Vt*yt
    inflation_cov_transform = Dinv*SVD.Vt*inflation_cov*SVD.V*Dinv

    objective = zeros((length(param1vals),length(param2vals)))
    for i = 1:length(param1vals)
        for j = 1:length(param2vals)
            diff = mean_mesh[:,i,j] - y_transform
            data_cov = Diagonal(var_mesh[:,i]) + inflation_cov_transform #Gamma + Delta
            gp_fidelity = 0.5*log(det(data_cov))
            gp_objective_fixed_tau[i] = 0.5*diff'*inv(data_cov)*diff + gp_fidelity
        end
    end
    minobj = minimum(objective)
    maxobj = maximum(objective)
    println("minimum of objective", minobj)
    println("maximum of objective", maxobj)
    param1_bd = [minimum(param1vals), maximum(param1vals)]
    param2_bd = [minimum(param2vals), maximum(param2vals)]

    #plot - transformed
    gr(size = (500,500))
    objlevels = collect(minobj:(maxobj - minobj) / (n_levels - 1):maxobj)

    #Plot the objective function
    plot(param1vals,
         param2vals,
         objective',
         seriestype=:contour,
         levels=objlevels,
         seriescolor=cgrad(:tempo), #white = small obj., dark = big obj.
         fill=true,
         legend=false,
         grid=false,
         xlims=param1_bd,
         ylims=param2_bd,
         dpi=300,
         colorbar=:right,
         clims=[minobj,maxobj],
         framestyle=:box,
         top_margin=25px,
         left_margin=50px,
         bottom_margin=50px)

    #plot the EKP points over the top
    circ = Shape(Plots.partialcircle(0, 2π))

    #get EKP points in domain
    param1_in_bd = (u_flat[:,1] .> param1_bd[1]) .& (u_flat[:,1].<param1_bd[2])
    param2_in_bd = (u_flat[:,2] .> param2_bd[1]) .& (u_flat[:,2].<param2_bd[2])
    u_in_bd = u_flat[(param1_in_bd) .& (param2_in_bd),:]

    plot!(u_in_bd[:,1],
          u_in_bd[:,2],
          seriestype=:scatter,
          markershape=circ,
          markercolor=:black,
          markersize=5,
          msw=0)

    xlabel!("transformed param1")
    ylabel!("transformed param2")
    savefig(outdir*"trainpts_over_objective_transformed.pdf")

    #untransformed 'raw' coordinates
    param1vals_raw,param2vals_raw = transform_prior_to_real([param1vals,param2vals])
    param1vals_raw,param2vals_raw = apply_units_transform([param1vals_raw,param2vals_raw])

    u_raw_in_bd = transform_prior_to_real(u_in_bd)
    u_raw_in_bd = apply_units_transform(u_raw_in_bd)

    param1_bd_raw = [minimum(param1vals_raw), maximum(param1vals_raw)]
    param2_bd_raw = [minimum(param2vals_raw), maximum(param2vals_raw)]

    #plot in untransformed coordinates

    #Plot the objective function

    plot(param1vals_raw,
         param2vals_raw,
         objective_raw', #may need transpose here
         #seriestype=:heatmap,
         seriestype=:contour,
         levels=objlevels,
         seriescolor=cgrad(:tempo), #white = small obj., dark = big obj.
         fill=true,
         legend=false,
         grid=false,
         xlims=param1_bd_raw,
         ylims=param2_bd_raw,
         dpi=300,
         colorbar=:right,
         clims=[minobj,maxobj],
         framestyle=:box,
         top_margin=25px,
         left_margin=50px,
         bottom_margin=50px)

    plot!(u_raw_in_bd[:,1],
          u_raw_in_bd[:,2],
          seriestype=:scatter,
          markershape=circ,
          markercolor=:black,
          markersize=5,
          msw=0)

    xlabel!("param1")
    ylabel!("param2")
    savefig("trainpts_over_objective_untransformed.pdf")

end

main()

ext: jl


text:
@testset "missing_data" begin

    @testset "nothing actually missing" begin
        # Construct a toy data set with no missing data.
        Y = randn(2, 5)
        y = ColVecs(Y)

        # Copy the toy data set, but change the type so that it _might_ have missing data.
        Y_missing = Matrix{Union{Missing, Float64}}(undef, size(Y))
        copy!(Y_missing, Y)
        y_missing = ColVecs(Y_missing)

        # Construct a random projection.
        U, s, _ = svd(randn(2, 1))
        S = Diagonal(s)

        # Construct noise model.
        σ² = 0.93
        D = Diagonal(abs.(randn(size(S, 1))) .+ 1e-2)

        # Ensure that both projection operations produce the same result.
        Yproj, ΣT = OILMMs.project(S, U, y, σ², D)
        Yproj_missing, ΣT_missing = OILMMs.project(S, U, y_missing, σ², D)
        @test Yproj ≈ Yproj_missing
        @test ΣT ≈ ΣT_missing

        # Ensure that the regularisation terms agree.
        reg = OILMMs.regulariser(S, U, σ², y_missing)
        @test reg ≈ OILMMs.regulariser(S, U, σ², y)

        @testset "OILMMs.project AD" begin

            # Perform forwards-pass and construct pullback.
            (Yproj_ad, ΣT_ad), project_pb = Zygote.pullback(
                OILMMs.project, S, U, y_missing, σ², D,
            )

            # Ensure that the forwards-pass is consistent with usual evaluation.
            @test Yproj_ad ≈ Yproj_missing
            @test ΣT_ad ≈ ΣT_missing

            # Estimate / evaluate cotangent of inputs.
            ΔYproj = randn(size(Yproj_ad))
            ΔΣT_proj = randn(size(ΣT_ad))
            Δout = (ΔYproj, ΔΣT_proj)
            dX_fd = FiniteDifferences.j′vp(
                central_fdm(5, 1),
                (S, U, σ², D) -> OILMMs.project(S, U, y_missing, σ², D),
                Δout, S, U, σ², D,
            )
            dX_ad = project_pb(Δout)

            # Check for (approximate) equality beteen AD and finite differencing.
            @test dX_fd[1] ≈ dX_ad[1]
            @test dX_fd[2] ≈ dX_ad[2]
            @test dX_fd[3] ≈ dX_ad[4]
            @test dX_fd[4] ≈ dX_ad[5]
        end

        @testset "OILMMs.regulariser AD" begin

            # Perform forwards-pass and construct pullback.
            reg_ad, regulariser_pd = Zygote.pullback(
                OILMMs.regulariser, S, U, σ², y_missing,
            )

            # Ensure that the forwards-pass is consistent with usual evaluation.
            @test reg_ad ≈ reg

            # Estimate / evaluate cotangent of inputs.
            Δreg = randn()
            dX_fd = FiniteDifferences.j′vp(
                central_fdm(5, 1),
                (S, U, σ²) -> OILMMs.regulariser(S, U, σ², y_missing),
                Δreg, S, U, σ²,
            )
            dX_ad = regulariser_pd(Δreg)

            # Check for (approximate) equality beteen AD and finite differencing.
            @test dX_fd[1] ≈ dX_ad[1]
            @test dX_fd[2] ≈ dX_ad[2]
            @test dX_fd[3] ≈ dX_ad[3]
        end
    end

    @testset "identity projection" begin

        # Construct a toy data set and make some bits of it missing.
        Y = Matrix{Union{Missing, Float64}}(undef, 3, 4)
        Y .= randn(3, 4)
        Y[1, 2] = missing
        Y[2, 1] = missing
        y = ColVecs(Y)

        # Compute the output of the project. The missings are just zeroed.
        Yproj = deepcopy(Y)
        Yproj[1, 2] = 0.0
        Yproj[2, 1] = 0.0

        # Construct the identity transformation.
        U, s, _ = svd(Matrix{Float64}(I, 3, 3))
        S = Diagonal(s)

        # Construct noise model.
        σ² = 0.93
        D = Diagonal(abs.(randn(size(Y, 1))) .+ 1e-2)

        # Ensure that the projection operation recovers the zeroed data.
        @test OILMMs.project(S, U, y, σ², D)[1] == Yproj
    end
end


text:
module TestIsLess
using Test
using CategoricalArrays

pool = CategoricalPool([1, 2, 3])

v1 = CategoricalValue(pool, 1)
v2 = CategoricalValue(pool, 2)
v3 = CategoricalValue(pool, 3)

@testset "values from unordered CategoricalPool" begin
    @test isordered(v1) === false
    @test isordered(v2) === false
    @test isordered(v3) === false

    @test_throws ArgumentError v1 < v1
    @test_throws ArgumentError v1 < v2
    @test_throws ArgumentError v1 < v3
    @test_throws ArgumentError v2 < v1
    @test_throws ArgumentError v2 < v2
    @test_throws ArgumentError v2 < v3
    @test_throws ArgumentError v3 < v1
    @test_throws ArgumentError v3 < v2
    @test_throws ArgumentError v3 < v3

    @test_throws ArgumentError v1 <= v1
    @test_throws ArgumentError v1 <= v2
    @test_throws ArgumentError v1 <= v3
    @test_throws ArgumentError v2 <= v1
    @test_throws ArgumentError v2 <= v2
    @test_throws ArgumentError v2 <= v3
    @test_throws ArgumentError v3 <= v1
    @test_throws ArgumentError v3 <= v2
    @test_throws ArgumentError v3 <= v3

    @test_throws ArgumentError v1 > v1
    @test_throws ArgumentError v1 > v2
    @test_throws ArgumentError v1 > v3
    @test_throws ArgumentError v2 > v1
    @test_throws ArgumentError v2 > v2
    @test_throws ArgumentError v2 > v3
    @test_throws ArgumentError v3 > v1
    @test_throws ArgumentError v3 > v2
    @test_throws ArgumentError v3 > v3

    @test_throws ArgumentError v1 >= v1
    @test_throws ArgumentError v1 >= v2
    @test_throws ArgumentError v1 >= v3
    @test_throws ArgumentError v2 >= v1
    @test_throws ArgumentError v2 >= v2
    @test_throws ArgumentError v2 >= v3
    @test_throws ArgumentError v3 >= v1
    @test_throws ArgumentError v3 >= v2
    @test_throws ArgumentError v3 >= v3

    @test isless(v1, v1) === false
    @test isless(v1, v2) === true
    @test isless(v1, v3) === true
    @test isless(v2, v1) === false
    @test isless(v2, v2) === false
    @test isless(v2, v3) === true
    @test isless(v3, v1) === false
    @test isless(v3, v2) === false
    @test isless(v3, v3) === false

    @testset "comparison with values of different types" begin
        @test_throws ArgumentError isless(v1, 1)
        @test_throws ArgumentError isless(v1, 2)
        @test_throws ArgumentError isless(v1, 10)
        @test_throws ArgumentError isless(v1, "a")
        @test_throws ArgumentError isless(1, v1)
        @test_throws ArgumentError isless("a", v1)
        @test_throws ArgumentError v1 < 1
        @test_throws ArgumentError v1 < 10
        @test_throws ArgumentError v1 < 2
        @test_throws ArgumentError v1 < "a"
        @test_throws ArgumentError v1 <= 1
        @test_throws ArgumentError v1 <= 2
        @test_throws ArgumentError v1 <= "a"
        @test_throws ArgumentError v1 > 1
        @test_throws ArgumentError v1 > 2
        @test_throws ArgumentError v1 > "a"
        @test_throws ArgumentError v1 >= 1
        @test_throws ArgumentError v1 >= 2
        @test_throws ArgumentError v1 >= "a"
    end

    @testset "comparison with missing" begin
        @test isless(v1, missing)
        @test !isless(missing, v1)
        @test ismissing(v1 < missing)
        @test ismissing(v1 <= missing)
        @test ismissing(v1 > missing)
        @test ismissing(v1 >= missing)
        @test ismissing(missing < v1)
        @test ismissing(missing <= v1)
        @test ismissing(missing > v1)
        @test ismissing(missing >= v1)
    end
end

@testset "values from ordered CategoricalPool" begin
    @test ordered!(pool, true) === pool
    @test isordered(pool) === true
    @test isordered(v1) === true
    @test isordered(v2) === true
    @test isordered(v3) === true

    @test (v1 < v1) === false
    @test (v1 < v2) === true
    @test (v1 < v3) === true
    @test (v2 < v1) === false
    @test (v2 < v2) === false
    @test (v2 < v3) === true
    @test (v3 < v1) === false
    @test (v3 < v2) === false
    @test (v3 < v3) === false

    @test (v1 <= v1) === true
    @test (v1 <= v2) === true
    @test (v1 <= v3) === true
    @test (v2 <= v1) === false
    @test (v2 <= v2) === true
    @test (v2 <= v3) === true
    @test (v3 <= v1) === false
    @test (v3 <= v2) === false
    @test (v3 <= v3) === true

    @test (v1 > v1) === false
    @test (v1 > v2) === false
    @test (v1 > v3) === false
    @test (v2 > v1) === true
    @test (v2 > v2) === false
    @test (v2 > v3) === false
    @test (v3 > v1) === true
    @test (v3 > v2) === true
    @test (v3 > v3) === false

    @test (v1 >= v1) === true
    @test (v1 >= v2) === false
    @test (v1 >= v3) === false
    @test (v2 >= v1) === true
    @test (v2 >= v2) === true
    @test (v2 >= v3) === false
    @test (v3 >= v1) === true
    @test (v3 >= v2) === true
    @test (v3 >= v3) === true

    @test isless(v1, v1) === false
    @test isless(v1, v2) === true
    @test isless(v1, v3) === true
    @test isless(v2, v1) === false
    @test isless(v2, v2) === false
    @test isless(v2, v3) === true
    @test isless(v3, v1) === false
    @test isless(v3, v2) === false
    @test isless(v3, v3) === false

    @testset "comparison with values of different types" begin
        @test_throws ArgumentError isless(v1, 1)
        @test_throws ArgumentError isless(v1, 2)
        @test_throws ArgumentError isless(v1, 10)
        @test_throws ArgumentError isless(v1, "a")
        @test_throws ArgumentError isless(1, v1)
        @test_throws ArgumentError isless("a", v1)
        @test_throws ArgumentError (v1 < 1)
        @test_throws ArgumentError (v1 < 2)
        @test_throws ArgumentError v1 < 10
        @test_throws ArgumentError v1 < "a"
        @test_throws ArgumentError (v1 <= 1)
        @test_throws ArgumentError (v1 <= 2)
        @test_throws ArgumentError v1 <= "a"
        @test_throws ArgumentError (v1 > 1)
        @test_throws ArgumentError (v1 > 2)
        @test_throws ArgumentError v1 > "a"
        @test_throws ArgumentError (v1 >= 1)
        @test_throws ArgumentError (v1 >= 2)
        @test_throws ArgumentError v1 >= "a"
    end

    @testset "comparison with values from different pools" begin
        poolb = copy(pool)
        v1b = CategoricalValue(poolb, 1)
        v2b = CategoricalValue(poolb, 2)
        v3b = CategoricalValue(poolb, 3)

        @test (v1 < v1b) === false
        @test (v1 < v2b) === true
        @test (v1 < v3b) === true
        @test (v2 < v1b) === false
        @test (v2 < v2b) === false
        @test (v2 < v3b) === true
        @test (v3 < v1b) === false
        @test (v3 < v2b) === false
        @test (v3 < v3b) === false

        @test isless(v1, v1b) === false
        @test isless(v1, v2b) === true
        @test isless(v1, v3b) === true
        @test isless(v2, v1b) === false
        @test isless(v2, v2b) === false
        @test isless(v2, v3b) === true
        @test isless(v3, v1b) === false
        @test isless(v3, v2b) === false
        @test isless(v3, v3b) === false

        poolc = copy(poolb)
        @test poolb == poolc # To set subsetof field
        v1c = CategoricalValue(poolc, 1)
        push!(poolc, 4)

        @test_throws ArgumentError v1b < v1c
        @test_throws ArgumentError isless(v1b, v1c)

        poolc = copy(poolb)
        @test poolb == poolc # To set subsetof field
        v1c = CategoricalValue(poolc, 1)
        push!(poolc, 4)

        @test_throws ArgumentError v1b < v1c
        @test_throws ArgumentError isless(v1b, v1c)

        poolc = copy(poolb)
        @test poolb == poolc # To set subsetof field
        v1c = CategoricalValue(poolc, 1)
        levels!(poolc, collect(1:4))

        @test_throws ArgumentError v1b < v1c
        @test_throws ArgumentError isless(v1b, v1c)
    end

    @testset "comparison with missing" begin
        @test isless(v1, missing)
        @test !isless(missing, v1)
        @test ismissing(v1 < missing)
        @test ismissing(v1 <= missing)
        @test ismissing(v1 > missing)
        @test ismissing(v1 >= missing)
        @test ismissing(missing < v1)
        @test ismissing(missing <= v1)
        @test ismissing(missing > v1)
        @test ismissing(missing >= v1)
    end
end

end


text:
using Random

mutable struct Experience{O,A}
    observation::O
    action::A
    reward::Float64
    value::Float64
end

mutable struct FixedSizeBuffer{O,A} <: AbstractBuffer
    experiences::Vector{Experience{O,A}}
    buffersize::Int
    tail::Int
    len::Int
    sample_order::Vector{Int}
    last_sample::Int

    function FixedSizeBuffer{O,A}(buffersize::Int, rng::MersenneTwister) where {O,A}
        v = Vector{Experience{O,A}}(undef, buffersize)

        new(v, buffersize, 1, 0, randperm(rng, buffersize), 1)
    end
end

function add!(buffer::FixedSizeBuffer{O,A}, obs::O, action::A, reward::Float64) where {O,A}
    buffer.experiences[buffer.tail] = Experience{O,A}(obs, action, reward, 0.0)
    buffer.tail = mod1(buffer.tail + 1, buffer.buffersize)
    buffer.len = buffer.len == buffer.buffersize ? buffer.buffersize : buffer.len + 1

    nothing
end

Base.length(buffer::FixedSizeBuffer) = buffer.len

# Should i conform to the AbstractArray interface?
get(buffer::FixedSizeBuffer, i::Int) = buffer.experiences[mod1(i, buffer.len)]

# Preciso de um mecanismo para iterar pelas últimas n experiências (para atualizar seus valores)
struct BufferIterator{O,A}
    buffer::FixedSizeBuffer{O,A}
    len::Int
end

function Base.iterate(iter::BufferIterator, state=1)
    state > iter.len && return nothing

    return (get(iter.buffer, iter.buffer.tail - state), state + 1)
end

Base.eltype(::Type{BufferIterator{O,A}}) where {O,A} = Experience{O,A}

Base.length(iter::BufferIterator) = iter.len

take_last(buffer::FixedSizeBuffer, n::Int) = BufferIterator(buffer, n)

# Preciso de um mecanismo para amostrar um minibatch do buffer

function sample!(buffer::FixedSizeBuffer{O,A}, batchsize::Int) where {O,A}
    batch = Vector{Experience{O,A}}(undef, batchsize)

        for i in 1:batchsize
        batch[i] = get(buffer, buffer.last_sample + i)
    end

    buffer.last_sample = mod1(buffer.last_sample + batchsize, buffer.buffersize)

    return batch
end


text:
#=
loglik.jl

    Log-likelihood calculation routines for a linear Gaussian State Space model

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2020/07/29
=#

"""
	loglik(filter)

Compute the log-likelihood for a linear Gaussian State Space model with Kalman
filter output `filter`.

#### Arguments
  - `filter::KalmanFilter`	: Kalman filter output

#### Returns
  - `ll::Real`	: log-likelihood
"""
function loglik(filter::MultivariateFilter)
    (n,T_len)= size(filter.v)

	# Initialize temp. containers
	tmp_nn= Matrix{Float64}(undef, (n,n))
	tmp_n= Vector{Float64}(undef, n)

    # Initialize log-likelihood
    ll= -log(2*π)*T_len*n

    # Log-likelihood
    @inbounds @fastmath for t in 1:T_len
		# Store views
        v_t= view(filter.v,:,t)
		F_t= view(filter.F,:,:,t)

		# Cholesky factorization of Fₜ
		copyto!(tmp_nn, F_t)
		C= cholesky!(Hermitian(tmp_nn))

		# log|Fₜ|
       	ll-= logdet(C)

		# Fₜ⁻¹ｘvₜ
        ldiv!(tmp_n, C, v_t)
		# vₜ'ｘFₜ⁻¹ｘvₜ
		ll-= dot(v_t, tmp_n)
    end

    return .5*ll
end

function loglik(filter::WoodburyFilter)
    (n,T_len)= size(filter.v)

	# Initialize temp. container
	tmp= Matrix{Float64}(undef, (n,n))

    # Initialize log-likelihood
    ll= -log(2*π)*T_len*n

    # Log-likelihood
    @inbounds @fastmath for t in 1:T_len
		# Store views
        v_t= view(filter.v,:,t)
		Fi_t= view(filter.Fi,:,:,t)

		# vₜ'ｘFₜ⁻¹ｘvₜ
		ll-= dot(v_t, Fi_t, v_t)

		# log|Fₜ⁻¹|
		copyto!(tmp, Fi_t)
       	ll+= logdet(cholesky!(Hermitian(tmp)))
    end

    return .5*ll
end

function loglik(filter::UnivariateFilter)
    (n,T_len)= size(filter.v)

    # Initialize log-likelihood
    ll= zero(Float64)

    # Log-likelihood
    @inbounds @fastmath for t in 1:T_len
		for i in 1:n
			F_it= filter.F[i,t]
			if F_it ≠ zero(F_it)
				ll-= log(2*π) + log(F_it) + filter.v[i,t]^2 * inv(F_it)
			end
		end
    end

	return .5*ll
end


text:
module InformedDifferentialEvolution

using FunctionalData, FunctionalDataUtils, Random

export DE

extract(a, field) = map(x->x[field],a)
nanmean(a) = mean(a[!isnan(a)])

function rounder!(a,stepsize::Real)
    for i = 1:length(a)
        a[i] = round(a[i]/stepsize)*stepsize
    end
end

function rounder!(a,stepsize::Matrix)
    for n = 1:len(a), m = 1:sizem(a)
        a[m,n] = round(a[m,n]/stepsize[m])*stepsize[m]
    end
end

DE(costf::Function, mi::Vector, ma::Vector; args...) = DE(costf, col(mi), col(ma); args...)
function DE(costf::Function, mi::Array{T,2}, ma::Array{T,2};
    npop = 100,
    maxiter::Int = 1_000_000,
    timeout::Number = Inf,
	maxstableiter::Int = 100,
    predictors::Array = Any[:default],
    sampler = rand,
    initpop = if sampler == rand
        col(mi) .+ rand(length(mi), npop) .* (col(ma) - col(mi))
    elseif sampler == randn
        (col(mi)+col(ma))/2 .+ clamp(randn(length(mi), npop)./6,-0.5,0.5) .* (col(ma) - col(mi))
    else
        error("InformedDifferentialEvolution: de: invalid value for paramter 'sampler'")
    end,
    recordhistory = false,
    tryallpredictors = false,
    continueabove = -Inf,
    crossoverprob = 0.5,
    diffweight = 0.85,
	roundto = 1e-6,
    data = nothing,
    verbosity::Array = Symbol[],  # :iter, :best, :newbest, :pop, :newbestcost, :initbestcost
    replaceworst = 0.1,
    classicmode = true,
    passpreviouscost = false,
    # :iter show iter number, :newbest show new bests, :pop population
    io = stdout,
    evaluator = nothing  # called as evaluator(pop, costs, bestind, best), can return arbitrary things for history
   ) where T

    pop = copy(initpop)
    rounder!(pop, roundto)
    clamp!(pop,mi,ma)
    costs = zeros(npop)
    newcosts = zeros(npop)
    frompredictor = falses(len(predictors), npop)
    ind1 = randperm(npop)
    ind2 = randperm(npop)
    iter = 1
    ncostevals = 0
	nstableiter = 0
	if isa(roundto, Array) roundto = col(roundto) end

    ###### predictors

    function defaultpredictor!(r::Matrix{T}, pop::Matrix{T}, costs) where T
        if classicmode
            randind() = rand(1:npop)
            for ind0 = 1:npop
                D = rand(1:sizem(pop))

                ind1 = randind()
                while ind1 == ind0 ind1 = randind() end
                ind2 = randind()
                while ind2 == ind0 || ind2 == ind1 ind2 = randind() end
                ind3 = randind()
                while ind3 == ind0 || ind3 == ind1 || ind3 == ind2 ind3 = randind() end

                for m = 1:sizem(pop)
                    if rand() < crossoverprob || m == D
                        r[m,ind0] = pop[m,ind1] + diffweight*(pop[m,ind2]-pop[m,ind3])
                    else
                        r[m,ind0] = pop[m,ind0]
                    end
                end
            end
            r
        else
            for n = 1:len(pop), m = 1:sizem(pop)
                ind1 = rand(1:npop)
                ind2 = rand(1:npop)
                r[m,n] = pop[m,n] + diffweight*(pop[m,ind1]-pop[m,ind2])
            end
            if replaceworst > 0.0
                # @show costs
                ind = sortperm(costs)
                n = round(length(costs)*replaceworst)
                r[:,ind[end-n+1:end]] = r[:, ind[1:n]]
            end
            r
        end
    end

    if isempty(predictors)
         error("no predictors specified, use [:default] for no predictors")
    end
    predictors = convert(Array{Any}, predictors)
    predictors[findall(predictors.==:default)] .= defaultpredictor!
    @assert all(x->isa(x, Function), predictors)  "InformedDifferentialEvolution: all predictors need to be functions or the symbol :default"
    predictedpop = zero(initpop)

    if !isnil(data)
        costf_(a) = costf(a, data)
    else
        costf_ = costf
    end

    for i = 1:npop
        costs[i] = costf_(pop[:,i])
        ncostevals += 1
    end
    bestcost, bestind = findmin(costs)
    best = pop[:,bestind]

    showiter = in(:iter, verbosity)
    showbest = in(:best, verbosity)
    shownewbest = in(:newbest, verbosity)
    shownewbestcost = in(:newbestcost, verbosity)
    showinitbestcost = in(:initbestcost, verbosity)
    showpop = in(:pop, verbosity)
    log(a...) = println(io, a...)

    if recordhistory
        history = Array(Any, maxiter+1)
        history[1] = Dict(:pop => copy(pop), :costs => copy(costs), :bestcost => copy(best), :bestind => bestind, :frompredictor => copy(frompredictor), :ncostevals => copy(ncostevals), :bestind => indmin(costs))
    else
        history = Any[]
    end

    # pv = FD.view(predictedpop)
    t = time()
    while iter <= maxiter && nstableiter < maxstableiter && bestcost > continueabove && time()-t < timeout
        showinitbestcost && iter == 1 && log("Best cost in init pop: $bestcost")
        showiter && log("Iteration: $iter")
        showpop && log("Population:\n$pop")

		nstableiter += 1

        gotbetter = false
        for m = 1:len(predictors)
            predictors[m](predictedpop, pop, costs)
            rounder!(predictedpop, roundto)
            clamp!(predictedpop, mi, ma)

            for n = 1:npop
                # FD.view!(predictedpop, n, pv)
                pv = at(predictedpop,n)
                if passpreviouscost
                    newcost = costf_(pv, costs[n])
                else
                    newcost = costf_(pv)
                end
                ncostevals += 1
                # @show costs[n] newcost pv
                if newcost < costs[n]
                    frompredictor[m,n] = true
                    pop[:,n] = pv
                    costs[n] = newcost
                end
            end
        end
        # println("^ gotbetter")
        oldbestcost = bestcost
        bestcost, bestind = findmin(costs)
        best = col(pop[:,bestind])

		if bestcost < oldbestcost
            shownewbestcost && log("New best at iter $iter with cost $bestcost")
            shownewbest && log("New best at iter $iter with cost $bestcost:\n", best)
			nstableiter = 0
		end

        if recordhistory
            history[iter+1] = Dict(:pop => copy(pop), :costs => copy(costs), :bestcost => copy(best), :frompredictor => copy(frompredictor), :ncostevals => copy(ncostevals), :bestind => bestind)
            if evaluator != nothing
                history[iter+1][:evaluator] = evaluator(pop, costs, bestind, best)
            end
        end
        iter += 1
    end
    if recordhistory
        history = history[1:iter]
        k = collect(keys(first(history)))
        history = Dict(zip(k, map(x->extract(history, x), k)))
    end
    showbest && log("Best:\n$best")
    best, Dict(:bestcost => bestcost, :ncostevals => ncostevals, :history => history, :pop => pop)
end

end # module InformedDifferentialEvolution

module InformedDEHelpers

using InformedDifferentialEvolution
export analyze, singleanalysis

function interpval(values, ticks, at)
    ind = findall(at.>=ticks)
    if isempty(ind)
        ind = 1
    else
        ind = last(ind)
    end
    #@show ind
    if at==ticks[ind]
        r = values[ind]
    elseif at>ticks[end]
        r = NaN
    else
        #@show values ticks at ind
        mixfactor = (at-ticks[ind])/(ticks[ind+1]-ticks[ind])
        r = mixfactor*values[ind+1] + (1-mixfactor)*values[ind]
    end
    #@show r
    r
end
@assert interpval(1:10, 1:10, 1.5) == 1.5
@assert interpval(-1:-1:-10, 1:10, 8.5) == -8.5
@assert map(x->interpval((1:3).^2,1:3,x),1:3) == [1,4,9]

function singleanalysis(f, mi, ma; nruns = 10, npop = 100, maxiter = 100, stepsize = 100, predictors = Any[], tryallpredictors = false, kargs...)
    nsamples = int((maxiter + 1)*max(1,length(predictors))*npop / stepsize)
    #@show nsamples
    costs = fill(NaN,nsamples)
    counters = zeros(nsamples)
    allbests = Any[]
    allstats = Any[]
    for i = 1:nruns
        best, stats = DE(f, mi, ma,
        tryallpredictors = tryallpredictors, recordhistory = true, maxiter = maxiter, npop = npop, predictors = predictors; kargs...)
        #@show stats[:bestcost] stats[:ncostevals]
        #@show stats[:history][:ncostevals]
        maxval = stats[:history][:ncostevals][end]
        for j = 1:maxval/stepsize
            #@show length(stats[:history][:bestcost])
            #@show length(stats[:history][:ncostevals])
            #@show j*stepsize
            isnan(costs[j]) ? costs[j] = 0 : nothing
            costs[j] += interpval(stats[:history][:bestcost], stats[:history][:ncostevals], j*stepsize)
            counters[j] += 1
        end
        push!(allbests, best)
        push!(allstats, stats)
    end
    ind = counters.>0
    costs[ind] ./= counters[ind]

    lines = semilogy(stepsize*(1:length(costs)),costs)
    xlabel("Number of cost function evaluations")
    ylabel("Cost of best idividuum")
    #lines = semilogy(stats[:history][:ncostevals], map(nanmean,stats[:history][:bestcost]))
    allbests, allstats, lines[1]
end

function analyze(f, mi, ma, predictors)
    r = map(x->singleanalysis(f, mi, ma, predictors = x[2]; x[3]...), predictors)
    legend(map(x->x[3],r), map(first, predictors))
end

function demo()
	f(x) = sum(abs.(x))
	predictor(pop,costs) = vcat(mean(pop,1), mean(pop,1))
	mi = [-1,-1]; ma = [1,1]

	analyze(f, mi, ma, [
		("default", [:default],Any[]),
		("pred first", [predictor, :default], Any[]),
		("pred last", [:default, predictor], Any[]),
		("pred only", [predictor], Any[]),
		("tryall", [predictor, :default], [(:tryallpredictors,true)])
	])
end

end # module


text:
# TODO: provide outer constructor checking dimensions?
"""
    RigidSolid(mass::Number, cg::Array{T, 1}, inertia::Array{T, 2})

Rigid solid represented by its mass center of gravity position reference system
and inertia tensor.
"""
struct RigidSolid
    mass::Number
    cg::Array{T, 1} where T<:Number
    inertia::Array{T, 2} where T<:Number
end

# Getters
get_mass(rs::RigidSolid) = rs.mass
get_cg(rs::RigidSolid) = rs.cg
get_inertia(rs::RigidSolid) = rs.inertia

"""
    PointMass(mass::Number, position::Array{T, 1})

Point mass represented by its mass and position. It is represented by a
`RigidSolid` with null inertia tensor.
"""
function PointMass(mass::Number, position::Array{T, 1}) where T<:Number
    RigidSolid(mass, position, zeros(3, 3))
end

"""
    get_inertia(rs::RigidSolid, point::Array{T, 1}) where T<:Number

Calculate inertia tensor of a `RigidSolid` object in a point P using Steiner's
theorem.
"""
function get_inertia(rs::RigidSolid, point::Array{T, 1}) where T<:Number
    steiner_inertia(rs.cg, rs.inertia, rs.mass, point)
end

"""
    +(comp1::RigidSolid, comp2::RigidSolid)

Compound two `RigidSolid` objects. Masses are added, a new cg is calculated and
the inertia tensors are added after being expressed with respect to the new
center of gravity using Steiner's theorem.
"""
function +(comp1::RigidSolid, comp2::RigidSolid)
    mass = comp1.mass + comp2.mass
    cg = (comp1.cg*comp1.mass + comp2.cg*comp2.mass) / (comp1.mass + comp2.mass)
    inertia = get_inertia(comp1, cg) + get_inertia(comp2, cg)
    RigidSolid(mass, cg, inertia)
end


text:
function kpp_psim(xi)

if xi>=0
  y=1+5*xi;
elseif (xi<0) & (xi>=-0.2)
  y=(1-16*xi)^(-1/4);
else
  y=(1.26-8.38*xi)^(-1/3);
end
return y
end


text:
#julia
using Makie, GeometryTypes, Colors
scene = Scene(resolution = (500, 500))
large_sphere = HyperSphere(Point3f0(0), 1f0)
positions = decompose(Point3f0, large_sphere)
linepos = view(positions, rand(1:length(positions), 1000))
lines(linepos, linewidth = 0.1, color = :black)
scatter(positions, strokewidth = 0.02, strokecolor = :white, color = RGBA(0.9, 0.2, 0.4, 0.6))
r = linspace(-1.5, 1.5, 5)
axis(r, r, r)
scene

#julia
using Makie, GeometryTypes, Colors
scene = Scene(resolution = (500, 500))
large_sphere = HyperSphere(Point3f0(0), 1f0)
positions = decompose(Point3f0, large_sphere)
meshscatter(positions, color = RGBA(0.9, 0.2, 0.4, 1))
scene

#julia
using Makie
scene = Scene(resolution = (500, 500))
function xy_data(x, y)
    r = sqrt(x*x + y*y)
    r == 0.0 ? 1f0 : (sin(r)/r)
end
r = linspace(-2, 2, 40)
surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
z = surf_func(20)
surf = surface(r, r, z)

wf = wireframe(r, r, surf[:z] .+ 1.0,
    linewidth = 2f0, color = lift_node(x-> x[5], surf[:colormap])
)
xy = linspace(-2.1, 2.1, 4)
axis(xy, xy, linspace(0, 2, 4))
center!(scene)

io = VideoStream(scene)
for i in linspace(0, 60, 100)
    surf[:z] = surf_func(i)
    recordframe!(io)
end
io

#julia
using Makie, GeometryTypes, GLVisualize, GLWindow
scene = Scene(resolution = (500, 500))
function xy_data(x, y)
    r = sqrt(x*x + y*y)
    r == 0.0 ? 1f0 : (sin(r)/r)
end
r = linspace(-2, 2, 40)
surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
N = 40
r = linspace(-2, 2, 40)
surface(
    r, r, surf_func(10),
    color = GLVisualize.loadasset("doge.png")
)
center!(scene)
scene

#julia
using Makie, GeometryTypes, GLVisualize
scene = Scene(resolution = (500, 500))
x = GLVisualize.loadasset("cat.obj")
Makie.mesh(x.vertices, x.faces, color = :black)
pos = map(x.vertices, x.normals) do p, n
    p => p .+ (normalize(n) .* 0.05f0)
end
linesegment(pos)
scene

#julia
using Makie
scene = Scene(resolution = (500, 500))
x = [0, 1, 2, 0]
y = [0, 0, 1, 2]
z = [0, 2, 0, 1]
color = [:red, :green, :blue, :yellow]
i = [0, 0, 0, 1]
j = [1, 2, 3, 2]
k = [2, 3, 1, 3]

indices = [1, 2, 3, 1, 3, 4, 1, 4, 2, 2, 3, 4]
m = mesh(x, y, z, indices, color = color)
wireframe(m[:mesh], color = :black, linewidth = 2)
r = linspace(-0.5, 2.5, 4)
axis(r, r, r)
center!(scene)
scene

#julia
using Makie, GLVisualize
scene = Scene(resolution = (500, 500))
mesh(GLVisualize.loadasset("cat.obj"))
r = linspace(-0.1, 1, 4)
center!(scene)
scene

#julia
using Makie, GeometryTypes, FileIO, GLVisualize
using GLVisualize: loadasset, assetpath
scene = Scene(resolution = (500, 500))
cat = load(assetpath("cat.obj"), GLNormalUVMesh)
Makie.mesh(cat, color = loadasset("diffusemap.tga"))
center!(scene)

#julia
using Makie, GeometryTypes
scene = Scene(resolution = (500, 500))
Makie.mesh(Sphere(Point3f0(0), 1f0))
center!(scene)
scene

#julia
using Makie, GeometryTypes, GLVisualize
scene = Scene(resolution = (500, 500))
wireframe(GLVisualize.loadasset("cat.obj"))
center!(scene)
scene

#julia
using Makie, GeometryTypes
scene = Scene(resolution = (500, 500))
wireframe(Sphere(Point3f0(0), 1f0))
center!(scene)
scene

#julia
using Makie
scene = Scene(resolution = (500, 500))
heatmap(rand(32, 32))
center!(scene)

#julias
using Makie, FileIO, GeometryTypes, Colors
scene = Scene(resolution = (500, 500), color = :black)
earth = load(download("https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/bluemarble-2048.png"))
image(earth)
center!(scene)

#julia
using Makie, FileIO, GeometryTypes, Colors
scene = Scene(resolution = (500, 500), color = :black)
m = GLNormalUVMesh(Sphere(Point3f0(0), 1f0), 60)
earth = load(download("https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/bluemarble-2048.png"))
Makie.mesh(m, color = earth)
stars = 100_000
scatter((rand(Point3f0, stars) .- 0.5) .* 10,
    glowwidth = 0.005, glow_color = :white, color = RGBA(0.8, 0.9, 0.95, 0.4),
    markersize = rand(linspace(0.0001, 0.01, 100), stars)
)
scene

#julia
using Makie
scene = Scene()
volume(rand(32, 32, 32), algorithm = :iso)
center!(scene)


text:
using Distributions
using Statistics
using Random

# This is a simple lookup table for the normal pdf
#
# Recall that the normal pdf has the form
#
# 1/(sigma sqrt(2pi)) exp(-1/2 ((x-mu) / sigma)^2)
struct NormalLookup
  # PDF values for 0:increment:max
  #
  # WARNING: since the zvalues start at zero, we have to be careful with zero
  # and one indexing
  pdfvals::Vector{Float64}

  # (unit: standard deviations)
  increment::Float64

  function NormalLookup(;increment, max)
    zvals = 0.0:increment:max

    N = Normal(0.0, 1.0)
    pdfvals = map(x->pdf(N, x), zvals)

    new(pdfvals, increment)
  end
end

const min_likelihood = 1e-6
function evaluate_no_norm(l::NormalLookup; error::Float64, stddev::Float64)
  if stddev < 0.0
    error("Negative standard deviation: ", stddev)
  end

  # rounding returns the effective zero-index so we add one
  index = round(Int64, abs(error) / (stddev * l.increment)) + 1

  # Since this is a normal pdf, take the limit for out-of-bounds values
  if index > length(l.pdfvals)
    min_likelihood
  else
    l.pdfvals[index]
  end
end
function evaluate(l::NormalLookup; error::Float64, stddev::Float64)
  evaluate_no_norm(l; error=error, stddev=stddev) / stddev
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
import Base: ∘

using AstroTime: Epoch

using LinearAlgebra: I
using ReferenceFrameRotations: compose_rotation
using StaticArrays: SMatrix, SDiagonal

export
    Rotation,
    ∘,
    compose,
    origin,
    target

struct Rotation{F1, F2, T, T′}
    m::SMatrix{3, 3, T}
    m′::SMatrix{3, 3, T′}
    function Rotation{F1, F2}(m, m′=zeros(3, 3)) where {F1, F2}
        T = eltype(m)
        T′ = eltype(m′)
        new{F1, F2, T, T′}(m, m′)
    end
end

function Rotation(origin::F, target::F, ::Epoch) where {F}
    Rotation{origin::AbstractFrame,
             target::AbstractFrame}(Array{Float64}(I, 3, 3), zeros(3, 3))
end

origin(::Rotation{F}) where {F} = F
target(::Rotation{_F, F}) where {_F, F} = F

function Base.inv(rot::Rotation{F1, F2}) where {F1, F2}
    Rotation{F2, F1}(rot.m', rot.m′')
end

(rot::Rotation)(r, v) = rot.m * r, rot.m′ * r + rot.m * v
(rot::Rotation)(rv) = rot.m * rv[1], rot.m′ * rv[1] + rot.m * rv[2]

function compose(rot1::Rotation{F1, F}, rot2::Rotation{F, F2}) where {F, F1, F2}
    Rotation{F1, F2}(rot2.m * rot1.m, rot2.m′ * rot1.m + rot2.m * rot1.m′)
end

∘(rot1::Rotation, rot2::Rotation) = compose(rot1, rot2)

function Rotation(frame1, frame2, ep::Epoch)
    frames = path_frames(frame1, frame2)
    rot = Rotation(from_sym(frames[1]), from_sym(frames[2]), ep)
    length(frames) == 2 && return rot

    for (f1, f2) in zip(frames[2:end-1], frames[3:end])
        rot = rot ∘ Rotation(from_sym(f1), from_sym(f2), ep)
    end
    rot
end


text:
module WriteToFile

include("./../constraints/geometric/constants.jl")
using Gurobi, JuMP
using DelimitedFiles
using Dates
using JSON
using JSON3
using .Constants

export write_to_file, write_to_json, append_to_file

function write_to_file(m)
    open(resultfile, "a") do io
        write(io, "========== NEW TEST ========== \n")

        write(io, "Date: ", string(Dates.today()), "\n")
        write(io, "JSON FILE: ", "test.json", "\n")
        write(io, "Time: ", "20", "\n")
        write(io, '\n')

        write(io, "Activated: \n")

        write(io, "Rotations Currence:        Y\n")
        #write(io, "Rotations Currence:        N\n")

        #write(io, "Rotations All:             Y\n")
        write(io, "Rotations All:             N\n")

        write(io, "Partial Base Stability:    Y\n")
        #write(io, "Partial Base Stability:    N\n")

        write(io, "Guillotine:                Y\n")
        #write(io, "Guillotine:                N\n")

        write(io, "Multiple Boxes:            Y\n")
        #write(io, "Multiple Boxes:            N\n")

        write(io, "Middle Incentive:          Y\n")
        #write(io, "Middle Incentive:          N\n")

        write(io, "Load Balancing:            Y\n")
        #write(io, "Load Balancing:            N\n")

        #write(io, "Partial Base Load Bearing: Y\n")
        write(io, "Partial Base Load Bearing: N\n")

        write(io, "Simple Load Bearing:       Y\n")
        #write(io, "Simple Load Bearing:       N\n")

        #write(io, "Number Load Bearing:       Y\n")
        write(io, "Simple Load Bearing:       N\n")

        #write(io, "Robot Pattern:             Y\n")
        write(io, "Robot Pattern:             N\n")
        write(io, '\n')

        write(io, "Positions: \n")
        positions =
        [
            [
            value(m[:flbx_front_left_bot_x][i]),
            value(m[:flby_front_left_bot_y][i]),
            value(m[:flbz_front_left_bot_z][i])
            ] for i in 1:length(boxes)
        ]
        writedlm(io, transpose(positions))
        write(io, '\n')

        write(io, "Static stabilities: \n")
        static_stabilities = [value(m[:piecewise_partial_static_stability_for_one_box][i]) for i in 1:length(boxes)]
        writedlm(io, transpose(static_stabilities))
        write(io, '\n')

        write(io, "Area covered: \n")
        area_covered = [(value(m[:piecewise_partial_static_stability_for_one_box][i])/area_of_box[i])*100 for i in 1:length(boxes)]
        writedlm(io, transpose(area_covered))
        write(io, '\n')

        write(io, "Guillotine Cut: \n")
        guillotine =
        [
            [
            value(m[:cut_flbx_left_i][i]),
            value(m[:cut_flbx_j_plus_left_i][i]),
            value(m[:cut_flbx_i_plus_right][i]),
            value(m[:cut_flbx_ij_plus_right][i]),
            value(m[:cut_flby_right][i]),
            value(m[:cut_flby_j_plus_right][i]),
            value(m[:cut_flby_i_plus_left][i]),
            value(m[:cut_flby_ij_plus_left][i]),
            ] for i in 1:length(boxes)
        ]
        writedlm(io, transpose(guillotine))
        write(io, '\n')

        write(io, "========== END OF TEST ========== \n")
        write(io, '\n')
        close(io)
    end
end

function write_to_json(m)
    positions =
        [
            [
            value(m[:flbx_front_left_bot_x][i]),
            value(m[:flby_front_left_bot_y][i]),
            value(m[:flbz_front_left_bot_z][i])
            ] for i in 1:length(boxes)
        ]
    static_stabilities = [value(m[:piecewise_partial_static_stability_for_one_box][i]) for i in 1:length(boxes)]
    area_covered = [(value(m[:piecewise_partial_static_stability_for_one_box][i])/area_of_box[i])*100 for i in 1:length(boxes)]
    guillotine =
        [
            [
            value(m[:cut_flbx_left_i][i]),
            value(m[:cut_flbx_j_plus_left_i][i]),
            value(m[:cut_flbx_i_plus_right][i]),
            value(m[:cut_flbx_ij_plus_right][i]),
            value(m[:cut_flby_right][i]),
            value(m[:cut_flby_j_plus_right][i]),
            value(m[:cut_flby_i_plus_left][i]),
            value(m[:cut_flby_ij_plus_left][i]),
            ] for i in 1:length(boxes)
        ]

    inside_json = []
    for i in 1:length(boxes)
        push!(inside_json, Dict(
                "x" => positions[i][1],
                "y" => positions[i][2],
                "z" => positions[i][3],
                "partial_base_coverage" => static_stabilities[i],
                "area_covered" => area_covered[i],
                "guillotine" =>
                    [
                    guillotine[i][1],
                    guillotine[i][2],
                    guillotine[i][3],
                    guillotine[i][4],
                    guillotine[i][5],
                    guillotine[i][6],
                    guillotine[i][7],
                    guillotine[i][7],
                    ]
                )
        )
    end
    dict = Dict("items" => inside_json)

    open(resultjson, "w") do io
        JSON3.pretty(io, dict)
    end

end

function append_to_file(m, mode::Int64)
    open(resultfile, "a") do io
        write(io, string(filename))
        write(io, '\n')
        write(io, string(mode))
        write(io, '\n')
        write(io, "TerminationStatus: ")
        write(io, string(termination_status(m)))
        write(io, '\n')
        write(io, "Time: ")
        write(io, string(solve_time(m)))
        write(io, string(termination_status(m)))
        if (has_values(m))
            write(io, '\n')
            write(io, "ObjectiveValue: ")
            write(io, string(objective_value(m)))
            write(io, '\n')
            write(io, "ObjectiveBound: ")
            write(io, string(objective_bound(m)))
            write(io, '\n')
            write(io, "Gap: ")
            write(io, string(relative_gap(m)))
            write(io, '\n')
            write(io, "ResultCount: ")
            write(io, string(result_count(m)))
            write(io, '\n')
            write(io, "NodeCount: ")
            write(io, string(node_count(m)))
        end
        write(io, '\n')
        write(io, '\n')
        close(io)
    end
end

end


text: *(t::ReTerm{Float64},v::Mamba.Stochastic{Vector{Float64}}) = t * v.value


text:
using CmdStanExtract
using Test

cnames_dummy = ["x", "y.1", "y.2", "z.1.1", "z.2.1", "z.3.1", "z.1.2", "z.2.2", "z.3.2", "k.1.1.1.1.1"]

key_to_idx = Dict(name => idx for (idx, name) in enumerate(cnames_dummy))

draws = 100
vars = length(cnames_dummy)
chains = 2
chns_dummy = randn(draws, vars, chains)

@testset "CmdStanExtract.jl" begin
    # Write your tests here.
    ex_dict = extract(chns_dummy, cnames_dummy)

    @test size(ex_dict["x"]) == (draws, chains)
    @test size(ex_dict["y"]) == (2, draws, chains)
    @test size(ex_dict["z"]) == (3, 2, draws, chains)
    @test size(ex_dict["k"]) == (1, 1, 1, 1, 1, draws, chains)

    @test ex_dict["x"][2,1] == chns_dummy[2, key_to_idx["x"], 1]
    @test ex_dict["y"][2,3,2] == chns_dummy[3, key_to_idx["y.2"], 2]
    @test ex_dict["z"][3, 1, 10, 1] == chns_dummy[10, key_to_idx["z.3.1"], 1]
    @test ex_dict["k"][1,1,1,1,1,draws,2] == chns_dummy[draws, key_to_idx["k.1.1.1.1.1"], 2]
end


text:
using AntisymmetricMatrices
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

include("antisymmetric.jl")

# write your own tests here


text:
function factored_norm(A::FacElem{NfOrdIdl, NfOrdIdlSet})
  b = Dict{fmpq, fmpz}()
  for (p, k) = A.fac
    n = norm(p)
    add_to_key!(b, fmpq(n), k)
    #if haskey(b, n)
    #  b[n] += k
    #else
    #  b[n] = k
    #end
  end
  bb = FacElem(FlintQQ, b)
  simplify!(bb)
  return bb
end

function norm(A::FacElem{NfOrdIdl, NfOrdIdlSet})
  return evaluate(factored_norm(A))
end

function factored_norm(A::NfOrdFracIdl)
  n = norm(A)
  return FacElem(Dict(fmpq(numerator(n)) => 1, fmpq(denominator(n)) => -1))
end

function factored_norm(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  b = Dict{fmpq, fmpz}()
  for (p, k) = A.fac
    if iszero(k)
      continue
    end
    n = norm(p)
    v = numerator(n)
    add_to_key!(b, fmpq(v), k)
    #if haskey(b, v)
    #  b[v] += k
    #else
    #  b[v] = k
    #end
    v1 = denominator(n)
    add_to_key!(b, fmpq(v1), -k)
    #if haskey(b, v)
    #  b[v] -= k
    #else
    #  b[v] = -k
    #end
  end
  bb = FacElem(FlintQQ, b)
  simplify!(bb)
  return bb
end

function norm(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  return evaluate(factored_norm(A))
end

@doc Markdown.doc"""
    valuation(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, p::NfOrdIdl)
    valuation(A::FacElem{NfOrdIdl, NfOrdIdlSet}, p::NfOrdIdl)
The valuation of $A$ at $P$.
"""
function valuation(A::FacElem{NfOrdIdl, NfOrdIdlSet}, p::NfOrdIdl)
  return sum(valuation(I, p)*v for (I, v) = A.fac if !iszero(v))
end

function valuation(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, p::NfOrdIdl)
  return sum(valuation(I, p)*v for (I, v) = A.fac)
end

@doc Markdown.doc"""
     ideal(O::NfOrd, a::FacElem{nf_elem, AnticNumberField)
The factored fractional ideal $a*O$.
"""
function ideal(O::NfOrd, a::FacElem{nf_elem, AnticNumberField})
  de = Dict{NfOrdFracIdl, fmpz}()
  for (e, k) = a.fac
    if !iszero(k)
      I = ideal(O, e)
      add_to_key!(de, I, k)
    end
  end
  return FacElem(FractionalIdealSet(O), de)
end

function ==(A::NfOrdIdl, B::FacElem{NfOrdIdl, NfOrdIdlSet})
  C = inv(B)*A
  return isone(C)
end
==(B::FacElem{NfOrdIdl, NfOrdIdlSet}, A::NfOrdIdl) = A == B

function ==(A::NfOrdFracIdl, B::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  C = A*inv(B)
  return isone(C)
end
==(B::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, A::NfOrdIdl) = A == B

function isone(A::NfOrdFracIdl)
  B = simplify(A)
  return B.den == 1 && isone(B.num)
end

function ==(A::FacElem{NfOrdFracIdl,NfOrdFracIdlSet}, B::FacElem{NfOrdFracIdl,NfOrdFracIdlSet})
  @assert check_parent(A, B) "Elements must have same parent"
  return isone(A*inv(B))
end
function ==(A::FacElem{NfOrdIdl,NfOrdIdlSet}, B::FacElem{NfOrdIdl,NfOrdIdlSet})
  @assert check_parent(A, B) "Elements must have same parent"
  return isone(A*inv(B))
end
function ==(A::FacElem{NfOrdIdl,NfOrdIdlSet}, B::FacElem{NfOrdFracIdl,NfOrdFracIdlSet})
  @assert order(base_ring(A)) === order(base_ring(B)) "Elements must be defined over the same order"
  return isone(A*inv(B))
end

==(A::FacElem{NfOrdFracIdl,NfOrdFracIdlSet}, B::FacElem{NfOrdIdl,NfOrdIdlSet}) = B==A

==(A::NfOrdFracIdl, B::FacElem{NfOrdIdl,NfOrdIdlSet}) = isone(A*inv(B))

function *(A::FacElem{NfOrdIdl,NfOrdIdlSet}, B::FacElem{NfOrdFracIdl,NfOrdFracIdlSet})
  @assert order(base_ring(A)) === order(base_ring(B)) "Elements must be defined over the same order"
  C = deepcopy(B)
  for (i,k) = A.fac
    C *= FacElem(Dict(i//1 => k))
  end
  return C
end
*(A::FacElem{NfOrdFracIdl,NfOrdFracIdlSet}, B::FacElem{NfOrdIdl,NfOrdIdlSet}) = B*A

function isone(A::FacElem{NfOrdIdl, NfOrdIdlSet})
  if all(x -> iszero(x), values(A.fac))
    return true
  end
  simplify!(A)
  return length(A.fac) == 1 && (isone(first(keys(A.fac))) || iszero(first(values(A.fac))))
end

function isone(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  A = simplify(A)
  return length(A.fac) == 1 && (isone(first(keys(A.fac))) || iszero(first(values(A.fac))))
end

function factor(Q::FacElem{NfOrdIdl, NfOrdIdlSet})
  if !all(isprime, keys(Q.fac))
    S = factor_coprime(Q)
    fac = Dict{NfOrdIdl, Int}()
    for (p, e)=S
      lp = factor(p)
      for (q, v) in lp
        fac[q] = Int(v*e)
      end
    end
  else
    fac = Dict(p=>Int(e) for (p,e) = Q.fac)
  end
  return fac
end

function FacElem(Q::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, O::NfOrdIdlSet)
  D = Dict{NfOrdIdl, fmpz}()
  for (I, v) = Q.fac
    if iszero(v)
      continue
    end
    if isone(I.den)
      add_to_key!(D, I.num, v)
    else
      n,d = integral_split(I)
      add_to_key!(D, n, v)
      add_to_key!(D, d, -v)
    end
  end
  return FacElem(O, D)
end

@doc Markdown.doc"""
    factor_coprime(Q::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}) -> Dict{NfOrdIdl, Int}
A coprime factorisation of $Q$: each ideal in $Q$ is split using \code{integral_split} and then
a coprime basis is computed.
This does {\bf not} use any factorisation.
"""
function factor_coprime(Q::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  D = FacElem(Q, NfOrdIdlSet(order(base_ring(Q))))
  S = factor_coprime(D)
  return S
end

@doc Markdown.doc"""
     factor(Q::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}) -> Dict{NfOrdIdl, Int}
The factorisation of $Q$, by refining a coprime factorisation.
"""
function factor(Q::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  S = factor_coprime(Q)
  fac = Dict{NfOrdIdl, Int}()
  for (p, e) = S
    lp = factor(p)
    for (q, v) in lp
      fac[q] = Int(v*e)
    end
  end
  return fac
end

#TODO: expand the coprime stuff to automatically also get the exponents
@doc Markdown.doc"""
    simplify(x::FacElem{NfOrdIdl, NfOrdIdlSet}) -> FacElem
    simplify(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}) -> FacElem

Uses ```coprime_base``` to obtain a simplified version of $x$, ie.
in the simplified version all base ideals will be pariwise coprime
but not neccessarily prime!.
"""
function simplify(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return z
end

function factor_over_coprime_base(x::FacElem{NfOrdIdl, NfOrdIdlSet}, coprime_base::Vector{NfOrdIdl})
  ev = Dict{NfOrdIdl, fmpz}()
  if isempty(coprime_base)
    return ev
  end
  OK = order(coprime_base[1])
  for p in coprime_base
    if isone(p)
      continue
    end
    P = minimum(p)
    @vprint :CompactPresentation 3 "Computing valuation at an ideal lying over $P"
    assure_2_normal(p)
    v = fmpz(0)
    for (b, e) in x
      if iszero(e)
        continue
      end
      if divisible(norm(b, copy = false), P)
        v += valuation(b, p)*e
      end
    end
    @vprint :CompactPresentation 3 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
    if !iszero(v)
      ev[p] = v
    end
  end
  return ev
end

function simplify!(x::FacElem{NfOrdIdl, NfOrdIdlSet}; refine::Bool = false)
  if length(x.fac) <= 1
    return nothing
  elseif all(x -> iszero(x), values(x.fac))
    x.fac = Dict{NfOrdIdl, fmpz}()
    return nothing
  end
  base_x = NfOrdIdl[y for (y, v) in x if !iszero(v)]
  cp = coprime_base(base_x, refine = refine)
  ev = factor_over_coprime_base(x, cp)
  x.fac = ev
  return nothing
end

function simplify(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return z
end

function simplify!(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  de = factor_coprime(x)
  if length(de)==0
    de = Dict(ideal(order(base_ring(parent(x))), 1) => fmpz(1))
  end
  x.fac = Dict((i//1, k) for (i,k) = de)
end

@doc Markdown.doc"""
    factor_coprime(x::FacElem{NfOrdIdl, NfOrdIdlSet}) -> Dict{NfOrdIdl, Int}
Computed a partial factorisation of $x$, ie. writes $x$ as a product
of pariwise coprime integral ideals.
"""
function factor_coprime(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return Dict{NfOrdIdl, Int}(p=>Int(v) for (p,v) = z.fac)
end

function factor_coprime!(x::FacElem{NfOrdIdl, NfOrdIdlSet}; refine::Bool = false)
  simplify!(x, refine = refine)
  return Dict{NfOrdIdl, Int}(p => Int(v) for (p,v) = x.fac)
end


text:
abstract AbstractNameManager
typealias NameType Union{Base.Symbol, AbstractString}
typealias NameCounter Dict{Base.Symbol, Int}

import Base: get!

# Default implementation for generating a name for a symbol.
# When a name is specified by the user, it will be used. Otherwise, a name
# is automatically generated based on the hint string.
function _default_get_name!(counter :: NameCounter, name :: NameType, hint :: NameType)
  if isa(name, Base.Symbol) || !isempty(name)
    return Symbol(name)
  end

  hint = Symbol(hint)
  if !haskey(counter, hint)
    counter[hint] = 0
  end
  name = Symbol("$hint$(counter[hint])")
  counter[hint] += 1
  return name
end

type BasicNameManager <: AbstractNameManager
  counter :: NameCounter
end
BasicNameManager() = BasicNameManager(NameCounter())

function get!(manager :: BasicNameManager, name :: NameType, hint :: NameType)
  _default_get_name!(manager.counter, name, hint)
end

type PrefixNameManager <: AbstractNameManager
  prefix  :: Base.Symbol
  counter :: NameCounter
end
PrefixNameManager(prefix :: NameType) = PrefixNameManager(Symbol(prefix), NameCounter())

function get!(manager :: PrefixNameManager, name :: NameType, hint :: NameType)
  name = _default_get_name!(manager.counter, name, hint)
  return Symbol("$(manager.prefix)$name")
end

DEFAULT_NAME_MANAGER = BasicNameManager()


text:
include("../../../PlotUtils/PlotUtils.jl")
include(joinpath(@__DIR__, "../../imports.jl"))
using PyPlot
const plt = PlotUtils.plt

function tomu(delta)
  mu0 = -cumsum(delta[0])
  mu1 = cumsum(delta[1])
  return [mu0; mu1]
end

function etas_tomat(etas)
  eta0s = [eta[0] for eta in etas]
  eta1s = [eta[1] for eta in etas]
  return (cat(eta0s..., dims=4),
          cat(eta1s..., dims=4))
end

# Results directory
results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-7-13/"

# Name of output file
OUTPUT_FILE = "output.bson"

# Path to all output files
output_paths = [joinpath(root, OUTPUT_FILE)
                for (root, _, files) in walkdir(results_dir)
                if OUTPUT_FILE in files]

# Load outputs
path_to_output = output_paths[2]
output = BSON.load(path_to_output)
samples = output[:samples][1]
m = output[:d].data.m
yimps = [o[:theta__y_imputed]
         for o in output[:samples][2]]
std(yimps[1][1][Bool.(m[1])])
std(yimps[10][1][Bool.(m[1])])

extract(sym) = [s[sym] for s in samples]

# mu*
mus = [tomu(d) for d in extract(:theta__delta)]
mu_mat = Matrix(hcat(mus...)')

# sig2
sig2s = extract(:theta__sig2)
sig2_mat = Matrix(hcat(sig2s...)')

# W
Ws = extract(:theta__W)
Wmat = cat(Ws..., dims=3)

# Z
Zs = extract(:theta__Z)

# R
Rs = hcat([vec(sum(W .> 0, dims=2)) for W in Ws]...)

# eta
eta0_mat, eta1_mat = etas_tomat(extract(:theta__eta))

# Plot (sig2, mu) vs eta
ell = 1
for i in 1:2, j in 1:20
  println((i, j))

  plt.subplot(2, 1, 1)
  plt.plot(eta0_mat[i, j, ell, :], mu_mat[:, 1], lw=0.5)
  plt.ylabel(L"$\mu^\star_{0, 1}$")
  plt.xlabel(latexstring("\\eta^0_{1, $(j), 1}"))

  plt.subplot(2, 1, 2)
  plt.plot(eta0_mat[i, j, ell, :], sig2_mat[:, i], lw=0.5)
  plt.ylabel(latexstring("\\sigma^2_{$(i)}"))
  plt.xlabel(latexstring("\\eta^0_{1, $(j), 1}"))

  plt.savefig("img/mu-vseta0_$(i)_$(j)_1.pdf", bbox_inches="tight")
  plt.close()
end

# plt.plot(output[:dden][10][1, 13])
# plt.plot(output[:dden][end][1, 13])
# plt.savefig("img/dden11.pdf")
# plt.close()


text:
# Multivariate optimization
function check_kwargs(kwargs, fallback_method)
    kws = Dict{Symbol, Any}()
    method = nothing
    for kwarg in kwargs
        if kwarg[1] != :method
            kws[kwarg[1]] = kwarg[2]
        else
            method = kwarg[2]
        end
    end

    if method == nothing
        method = fallback_method
    end
    kws, method
end

default_options(method::AbstractOptimizer) = Dict{Symbol, Any}()

function add_default_opts!(opts::Dict{Symbol, Any}, method::AbstractOptimizer)
    for newopt in default_options(method)
        if !haskey(opts, newopt[1])
            opts[newopt[1]] = newopt[2]
        end
    end
end

fallback_method(f) = NelderMead()
fallback_method(f, g!) = LBFGS()
fallback_method(f, g!, h!) = Newton()

fallback_method(d::OnceDifferentiable) = LBFGS()
fallback_method(d::TwiceDifferentiable) = Newton()

# promote the objective (tuple of callables or an AbstractObjective) according to method requirement
promote_objtype(method, initial_x, autodiff::Symbol, inplace::Bool, args...) = error("No default objective type for $method and $args.")
# actual promotions, notice that (args...) captures FirstOrderOptimizer and NonDifferentiable, etc
promote_objtype(method::ZerothOrderOptimizer, x, autodiff::Symbol, inplace::Bool, args...) = NonDifferentiable(args..., x, real(zero(eltype(x))))
promote_objtype(method::FirstOrderOptimizer,  x, autodiff::Symbol, inplace::Bool, f) = OnceDifferentiable(f, x, real(zero(eltype(x))); autodiff = autodiff)
promote_objtype(method::FirstOrderOptimizer,  x, autodiff::Symbol, inplace::Bool, args...) = OnceDifferentiable(args..., x, real(zero(eltype(x))); inplace = inplace)
promote_objtype(method::FirstOrderOptimizer,  x, autodiff::Symbol, inplace::Bool, f, g, h) = OnceDifferentiable(f, g, x, real(zero(eltype(x))); inplace = inplace)
promote_objtype(method::SecondOrderOptimizer, x, autodiff::Symbol, inplace::Bool, f) = TwiceDifferentiable(f, x, real(zero(eltype(x))); autodiff = autodiff)
promote_objtype(method::SecondOrderOptimizer, x, autodiff::Symbol, inplace::Bool, f::NotInplaceObjective) = TwiceDifferentiable(f, x, real(zero(eltype(x))))
promote_objtype(method::SecondOrderOptimizer, x, autodiff::Symbol, inplace::Bool, f::InplaceObjective) = TwiceDifferentiable(f, x, real(zero(eltype(x))))
promote_objtype(method::SecondOrderOptimizer, x, autodiff::Symbol, inplace::Bool, f, g) = TwiceDifferentiable(f, g, x, real(zero(eltype(x))); inplace = inplace, autodiff = autodiff)
promote_objtype(method::SecondOrderOptimizer, x, autodiff::Symbol, inplace::Bool, f, g, h) = TwiceDifferentiable(f, g, h, x, real(zero(eltype(x))); inplace = inplace)
# no-op
promote_objtype(method::ZerothOrderOptimizer, x, autodiff::Symbol, inplace::Bool, nd::NonDifferentiable)  = nd
promote_objtype(method::ZerothOrderOptimizer, x, autodiff::Symbol, inplace::Bool, od::OnceDifferentiable) = od
promote_objtype(method::FirstOrderOptimizer,  x, autodiff::Symbol, inplace::Bool, od::OnceDifferentiable) = od
promote_objtype(method::ZerothOrderOptimizer, x, autodiff::Symbol, inplace::Bool, td::TwiceDifferentiable) = td
promote_objtype(method::FirstOrderOptimizer,  x, autodiff::Symbol, inplace::Bool, td::TwiceDifferentiable) = td
promote_objtype(method::SecondOrderOptimizer, x, autodiff::Symbol, inplace::Bool, td::TwiceDifferentiable) = td

# if no method or options are present
function optimize(f,         initial_x::AbstractArray; inplace = true, autodiff = :finite, kwargs...)
    method = fallback_method(f)
    checked_kwargs, method = check_kwargs(kwargs, method)
    d = promote_objtype(method, initial_x, autodiff, inplace, f)
    add_default_opts!(checked_kwargs, method)

    options = Options(; checked_kwargs...)
    optimize(d, initial_x, method, options)
end
function optimize(f, g, initial_x::AbstractArray; inplace = true, autodiff = :finite, kwargs...)

    method = fallback_method(f, g)
    checked_kwargs, method = check_kwargs(kwargs, method)
    d = promote_objtype(method, initial_x, autodiff, inplace, f, g)
    add_default_opts!(checked_kwargs, method)

    options = Options(; checked_kwargs...)
    optimize(d, initial_x, method, options)
end
function optimize(f, g, h, initial_x::AbstractArray; inplace = true, autodiff = :finite, kwargs...)

    method = fallback_method(f, g, h)
    checked_kwargs, method = check_kwargs(kwargs, method)
    d = promote_objtype(method, initial_x, autodiff, inplace, f, g, h)
    add_default_opts!(checked_kwargs, method)

    options = Options(; checked_kwargs...)
    optimize(d, initial_x, method, options)
end

# no method supplied with objective
function optimize(d::T, initial_x::AbstractArray, options::Options) where T<:AbstractObjective
    optimize(d, initial_x, fallback_method(T), options)
end
# no method supplied with inplace and autodiff keywords becauase objective is not supplied
function optimize(f, initial_x::AbstractArray, options::Options; inplace = true, autodiff = :finite)
    method = fallback_method(f)
    d = promote_objtype(method, initial_x, autodiff, inplace, f)
    optimize(d, initial_x, method, options)
end
function optimize(f, g, initial_x::AbstractArray, options::Options; inplace = true, autodiff = :finite)

    method = fallback_method(f, g!)
    d = promote_objtype(method, initial_x, autodiff, inplace, f, g)
    optimize(d, initial_x, method, options)
end
function optimize(f, g, h, initial_x::AbstractArray, options::Options; inplace = true, autodiff = :finite)

    method = fallback_method(f, g!, h!)
    d = promote_objtype(method, initial_x, autodiff, inplace, f, g, h)

    optimize(d, initial_x, method, options)
end

# potentially everything is supplied (besides caches)
function optimize(f, initial_x::AbstractArray, method::AbstractOptimizer,
                     options::Options = Options(;default_options(method)...); inplace = true, autodiff = :finite)

    d = promote_objtype(method, initial_x, autodiff, inplace, f)
    optimize(d, initial_x, method, options)
end
function optimize(f, g, initial_x::AbstractArray, method::AbstractOptimizer,
         options::Options = Options(;default_options(method)...); inplace = true, autodiff = :finite)

    d = promote_objtype(method, initial_x, autodiff, inplace, f, g)

    optimize(d, initial_x, method, options)
end
function optimize(f, g, h, initial_x::AbstractArray{T}, method::AbstractOptimizer,
         options::Options = Options(;default_options(method)...); inplace = true, autodiff = :finite) where T

    d = promote_objtype(method, initial_x, autodiff, inplace, f, g, h)

    optimize(d, initial_x, method, options)
end

function optimize(d::D, initial_x::AbstractArray, method::SecondOrderOptimizer,
                  options::Options = Options(;default_options(method)...); autodiff = :finite, inplace = true) where {D <: Union{NonDifferentiable, OnceDifferentiable}}
    d = promote_objtype(method, initial_x, autodiff, inplace, d)
    optimize(d, initial_x, method, options)
end


text:
# initializing and optimizing embeddings

function initialize_embedding(rng::AbstractRNG, graph::AbstractMatrix{T}, n_components, ::Val{:spectral}) where {T}
    local embed
    try
        embed = spectral_layout(graph, n_components)
        # expand
        expansion = 10 / maximum(embed)
        embed .= (embed .* expansion) .+ (1//10000) .* randn.(rng, T)
    catch e
        @info "$e\nError encountered in spectral_layout; defaulting to random layout"
        embed = initialize_embedding(rng, graph, n_components, Val(:random))
    end
    return embed
end

function initialize_embedding(rng::AbstractRNG, graph::AbstractMatrix{T}, n_components, ::Val{:random}) where {T}
    20 .* rand(rng, T, (n_components, size(graph, 1))) .- 10
end

"""
    initialize_embedding(rng::AbstractRNG, graph::AbstractMatrix{<:Real}, ref_embedding::AbstractMatrix{T<:AbstractFloat}) -> embedding

Initialize an embedding of points corresponding to the columns of the `graph`, by taking weighted means of
the columns of `ref_embedding`, where weights are values from the rows of the `graph`.

The resulting embedding will have shape `(size(ref_embedding, 1), size(graph, 2))`, where `size(ref_embedding, 1)`
is the number of components (dimensions) of the `reference embedding`, and `size(graph, 2)` is the number of
samples in the resulting embedding. Its elements will have type T.
"""
function initialize_embedding(::AbstractRNG, graph::AbstractMatrix{<:Real}, ref_embedding::AbstractMatrix{T})::Matrix{T} where {T<:AbstractFloat}
    (ref_embedding * graph) ./ (sum(graph, dims=1) .+ eps(T))
end

"""
    spectral_layout(graph, embed_dim) -> embedding

Initialize the graph layout with spectral embedding.
"""
function spectral_layout(graph::SparseMatrixCSC{T},
                         embed_dim::Integer) where {T<:Real}
    graph_f64 = convert.(Float64, graph)
    D_ = Diagonal(dropdims(sum(graph_f64; dims=2); dims=2))
    D = inv(sqrt(D_))
    # normalized laplacian
    L = Symmetric(I - D*graph*D)

    k = embed_dim+1
    num_lanczos_vectors = max(2k+1, round(Int, sqrt(size(L, 1))))
    # get the 2nd - embed_dim+1th smallest eigenvectors
    eigenvals, eigenvecs = eigs(L; nev=k,
                                   ncv=num_lanczos_vectors,
                                   which=:SM,
                                   tol=1e-4,
                                   v0=ones(Float64, size(L, 1)),
                                   maxiter=size(L, 1)*5)
    layout = permutedims(eigenvecs[:, 2:k])::Array{Float64, 2}
    return convert.(T, layout)
end

"""
    optimize_embedding(rng, graph, query_embedding, ref_embedding, n_epochs, initial_alpha, min_dist, spread, gamma, neg_sample_rate, _a=nothing, _b=nothing; move_ref=false) -> embedding

Optimize an embedding by minimizing the fuzzy set cross entropy between the high and low dimensional simplicial sets using stochastic gradient descent.
Optimize "query" samples with respect to "reference" samples.

# Arguments
- `graph`: a sparse matrix of shape (n_samples, n_samples)
- `query_embedding`: a matrix (n_components, n_samples) the embedded data points to be optimized ("query" samples)
- `ref_embedding`: a matrix (n_components, n_samples) the embedded data points to be optimized against ("reference" samples)
- `n_epochs`: the number of training epochs for optimization
- `initial_alpha`: the initial learning rate
- `gamma`: the repulsive strength of negative samples
- `neg_sample_rate`: the number of negative samples per positive sample
- `_a`: this controls the embedding. If the actual argument is `nothing`, this is determined automatically by `min_dist` and `spread`.
- `_b`: this controls the embedding. If the actual argument is `nothing`, this is determined automatically by `min_dist` and `spread`.

# Keyword Arguments
- `move_ref::Bool = false`: if true, also improve the embeddings in `ref_embedding`, else fix them and only improve embeddings in `query_embedding`.
"""
function optimize_embedding(rng::AbstractRNG,
                            graph::SparseMatrixCSC{T},
                            query_embedding::AbstractMatrix{S},
                            ref_embedding::AbstractMatrix{S},
                            n_epochs::Integer,
                            initial_alpha::Real,
                            min_dist::Real,
                            spread::Real,
                            gamma::Real,
                            neg_sample_rate::Integer,
                            _a::Union{Real, Nothing} = nothing,
                            _b::Union{Real, Nothing} = nothing;
                            move_ref::Bool=false) where {T <: Real, S <: Real}
    a, b = fit_ab(convert(S,min_dist), convert(S, spread), _a, _b)
    self_reference = query_embedding === ref_embedding

    a_2_b = -2 * a * b
    gamma_2_b = 2 * gamma * b

    alpha = initial_alpha
    for e in 1:n_epochs
        @inbounds for i in 1:size(graph, 2)
            for ind in nzrange(graph, i)
                j = rowvals(graph)[ind]
                p = nonzeros(graph)[ind]

                if rand(rng) <= p
                    xi = view(query_embedding, :, i)
                    xj = view(ref_embedding, :, j)

                    sdist = evaluate(SqEuclidean(), xi, xj)
                    sdist = max(sdist, eps(S))
                    delta = (a_2_b * sdist^b)/(sdist*(1 + a*sdist^b))

                    @simd for d in eachindex(xi)
                        grad = clamp(delta * (xi[d] - xj[d]), -4, 4)
                        xi[d] += alpha * grad
                        if move_ref
                            xj[d] -= alpha * grad
                        end
                    end

                    for _ in 1:neg_sample_rate
                        k = rand(rng, 1:size(ref_embedding, 2))
                        if i == k && self_reference
                            continue
                        end
                        xk = view(ref_embedding, :, k)

                        sdist = evaluate(SqEuclidean(), xi, xk)
                        sdist = max(sdist, eps(S))
                        delta = gamma_2_b / ((S(1//1000) + sdist)*(1 + a*sdist^b))

                        @simd for d in eachindex(xi)
                            grad = clamp(delta * (xi[d] - xk[d]), -4, 4)
                            xi[d] += alpha * grad
                        end
                    end

                end
            end
        end
        alpha = initial_alpha*(1 - e//n_epochs)
    end

    return query_embedding
end