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 Dierckx  # cubic b-spline for airfoil data
# import Roots  # solve residual equation
import Parameters: @unpack
import PyCall
import Interpolations
import ForwardDiff

export Rotor, OperatingPoint, Outputs
export af_from_file, af_from_data
export simple_op, windturbine_op
export solve, thrusttorque, nondim

include("openmdao.jl")

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

"""
    Rotor(r, chord, theta, af, Rhub, Rtip, B, turbine, precone=0.0)

Define rotor geometry

**Arguments**
- `r::Array{Float64, 1}`: radial locations along turbine (Rhub < r < Rtip)
- `chord::Array{Float64, 1}`: corresponding local chord lengths
- `theta::Array{Float64, 1}`: corresponding twist angles (radians)
- `af::Array{function, 1}`: a function of the form: cl, cd = af(alpha, Re, Mach)
- `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
- `precone::Float64`: precone angle
"""
struct Rotor{TF, TAF, TI, TB, TA<:AbstractArray{TF, 1}, TAAF<:AbstractArray{TAF, 1}}

    r::TA
    chord::TA
    theta::TA
    af::TAAF
    Rhub::TF
    Rtip::TF
    B::TI
    turbine::TB
    precone::TF

end

# convenience constructor for no precone
Rotor(r, chord, theta, af, Rhub, Rtip, B, turbine) = Rotor(r, chord, theta, af, Rhub, Rtip, B, turbine, 0.0)

# constructor for all keyword arguments, all scalars, ignoring extras
function Rotor(; r::TF, chord::TF, theta::TF, af::TAF, Rhub::TF, Rtip::TF, B::TI, turbine::TB, precone::TF, kwargs...) where {TF, TAF, TI, TB}
    Rotor(TF[r], TF[chord], TF[theta], TAF[af], Rhub, Rtip, B, turbine, precone)
end

"""
    OperatingPoint(Vx, Vy, pitch, 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::Array{Float64, 1}`: velocity in x-direction along blade
- `Vy::Array{Float64, 1}`: velocity in y-direction along blade
- `pitch::Float64`: pitch angle (rad).  defined same direction as twist.
- `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)
"""
struct OperatingPoint{TF, TA<:AbstractArray{TF, 1}}
    Vx::TA
    Vy::TA
    pitch::TF
    rho::TF
    mu::TF
    asound::TF
end

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

# constructor for all keyword arguments, ignoring extras
function OperatingPoint(; Vx::TF, Vy::TF, pitch::TF, rho::TF, mu::TF, asound::TF, kwargs...)  where {TF}
    OperatingPoint(TF[Vx], TF[Vy], pitch, rho, mu, asound)
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::Array{Float64, 1}`: normal force per unit length
- `Tp::Array{Float64, 1}`: tangential force per unit length
- `a::Array{Float64, 1}`: axial induction factor
- `ap::Array{Float64, 1}`: tangential induction factor
- `u::Array{Float64, 1}`: axial induced velocity
- `v::Array{Float64, 1}`: tangential induced velocity
- `phi::Array{Float64, 1}`: inflow angle
- `alpha::Array{Float64, 1}`: angle of attack
- `W::Array{Float64, 1}`: inflow velocity
- `cl::Array{Float64, 1}`: lift coefficient
- `cd::Array{Float64, 1}`: drag coefficient
- `cn::Array{Float64, 1}`: normal force coefficient
- `ct::Array{Float64, 1}`: tangential force coefficient
- `F::Array{Float64, 1}`: hub/tip loss correction
- `G::Array{Float64, 1}`: effective hub/tip loss correction for induced velocities: u = Vx * a * G, v = Vy * ap * G
"""
struct Outputs{TF}

    Np::Array{TF, 1}
    Tp::Array{TF, 1}
    a::Array{TF, 1}
    ap::Array{TF, 1}
    u::Array{TF, 1}
    v::Array{TF, 1}
    phi::Array{TF, 1}
    alpha::Array{TF, 1}
    W::Array{TF, 1}
    cl::Array{TF, 1}
    cd::Array{TF, 1}
    cn::Array{TF, 1}
    ct::Array{TF, 1}
    F::Array{TF, 1}
    G::Array{TF, 1}

end

# convenience constructor to initialize
function Outputs(etype, nr)
    return Outputs(
        zeros(etype, nr), zeros(etype, nr), zeros(etype, nr), zeros(etype, nr),
        zeros(etype, nr), zeros(etype, nr), zeros(etype, nr), zeros(etype, nr),
        zeros(etype, nr), zeros(etype, nr), zeros(etype, nr), zeros(etype, nr),
        zeros(etype, nr), zeros(etype, nr), zeros(etype, nr)
    )
end

struct ResidualPartials{TF}
    phi::TF
    r::TF
    chord::TF
    theta::TF
    Vx::TF
    Vy::TF
    rho::TF
    mu::TF
    asound::TF
    Rhub::TF
    Rtip::TF
    precone::TF
end

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

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

"""
    af_from_file(filename)

Read an airfoil file.
Currently only reads one Reynolds 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 a function of the form `cl, cd = func(alpha, Re, M)` although Re and M are currently ignored.
"""
function af_from_file(filename; use_interpolations_jl=false)

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

    open(filename) do f

        # skip header
        readline(f)

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

    end

    return af_from_data(alpha*pi/180.0, cl, cd; use_interpolations_jl=use_interpolations_jl)
end

"""
    af_from_data(alpha, cl, cd)

Create an AirfoilData object directly from alpha, cl, and cd arrays.
alpha should be in radians.

`af_from_file` calls this function indirectly.  Uses a cubic B-spline
(if the order of the data permits it).  A small amount of smoothing of
lift and drag coefficients is also applied to aid performance
for gradient-based optimization.

Returns a function of the form `cl, cd = func(alpha, Re, M)` although Re and M are currently ignored.
"""
function af_from_data(alpha, cl, cd, spl_k=3; use_interpolations_jl=false)

    # # TODO: update once smoothing is implemented: https://github.com/JuliaMath/Interpolations.jl/issues/254
    # afcl = CubicSplineInterpolation(alpha, cl)
    # afcd = CubicSplineInterpolation(alpha, cd)
    # af = AirfoilData(afcl, afcd)

    k = min(length(alpha)-1, spl_k)  # can't use cubic spline is number of entries in alpha is small

    # 1D interpolations for now.  ignoring Re dependence (which is very minor)
    afcl_1d = Dierckx.Spline1D(alpha, cl; k=k, s=0.1)
    afcd_1d = Dierckx.Spline1D(alpha, cd; k=k, s=0.001)

    if use_interpolations_jl

        # Check if the original alpha is uniformly distributed.
        dalpha = alpha[2:end] - alpha[1:end-1]
        tol = 1e-8
        if all(@. abs(dalpha - dalpha[1]) < tol)
            # Uniform spacing, so no need to reinterpolate the Dierckx
            # interpolation onto a uniform grid.

            # println("detected uniform angle of attack sampling.")
            n_uniform = length(alpha)
            alpha_uniform = LinRange(minimum(alpha), maximum(alpha), n_uniform)
            afcl_1d = Interpolations.CubicSplineInterpolation(alpha_uniform, cl, extrapolation_bc=Interpolations.Periodic())
            afcd_1d = Interpolations.CubicSplineInterpolation(alpha_uniform, cd, extrapolation_bc=Interpolations.Periodic())

        else
            # Evaluate the Dierckx interpolation object on a uniform grid.
            n_uniform = 2*length(alpha)
            alpha_uniform = LinRange(minimum(alpha), maximum(alpha), n_uniform)
            cl_uniform = afcl_1d(alpha_uniform)
            cd_uniform = afcd_1d(alpha_uniform)

            # Use the uniform data to create a new interpolations objects using
            # Interpolations.jl.
            afcl_1d = Interpolations.CubicSplineInterpolation(alpha_uniform, cl_uniform, extrapolation_bc=Interpolations.Periodic())
            afcd_1d = Interpolations.CubicSplineInterpolation(alpha_uniform, cd_uniform, extrapolation_bc=Interpolations.Periodic())
        end

    end

    afeval(alpha, Re, M) = afcl_1d(alpha), afcd_1d(alpha)  # ignore Re, M

    return afeval
end

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

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

"""
(private) residual function
modifies outputs if setoutputs=true
"""
function residual!(phi, rotor, inflow, outputs, idx, setoutputs=false)

    # unpack inputs
    r = rotor.r[idx]
    chord = rotor.chord[idx]
    theta = rotor.theta[idx]
    af = rotor.af[idx]
    Rhub = rotor.Rhub
    Rtip = rotor.Rtip
    B = rotor.B
    Vx = inflow.Vx[idx]
    Vy = inflow.Vy[idx]
    pitch = inflow.pitch
    rho = inflow.rho

    # check if turbine or propeller and change input sign if necessary
    setsign = rotor.turbine ? 1 : -1
    theta *= setsign
    pitch *= setsign
    Vx *= setsign

    # 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 / inflow.mu

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

    # airfoil cl/cd
    cl, cd = af(alpha, Re, Mach)

    # 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
    # F::typeof(Ftip) = 1.0

    # 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

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

    #     u = 0.0
    #     v = k0p*abs(Vx)
    #     a = 0.0
    #     ap = 0.0

    # 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

    # end

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

    # fill outputs
    if setoutputs

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

        # ---- swap output signs as needed -----
        Tp *= setsign
        v *= setsign

        # 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.
        G = (1.0 - sqrt(1.0 - 4*a*(1.0 - a)*F))/(2*a)
        u *= G
        v *= G

        outputs.Np[idx] = Np
        outputs.Tp[idx] = Tp
        outputs.a[idx] = a
        outputs.ap[idx] = ap
        outputs.u[idx] = u
        outputs.v[idx] = v
        outputs.phi[idx] = phi
        outputs.alpha[idx] = alpha
        outputs.W[idx] = sqrt(W2)
        outputs.cl[idx] = cl
        outputs.cd[idx] = cd
        outputs.cn[idx] = cn
        outputs.ct[idx] = ct
        outputs.F[idx] = F
        outputs.G[idx] = G
    end

    return R

end

# function residual(B, af, turbine, inputs::AbstractArray)
#     phi = inputs[1]
#     r = inputs[2]
#     chord = inputs[3]
#     theta = inputs[4]
#     Vx, Vy, rho, mu, asound = inputs[5], inputs[6], inputs[7], inputs[8], inputs[9]
#     Rhub, Rtip, precone = inputs[10], inputs[11], inputs[12]

#     section = Section(r, chord, theta, af)
#     inflow = Inflow(Vx, Vy, rho, mu, asound)
#     rotor = Rotor(Rhub, Rtip, B, turbine, precone)

#     R, outputs = residual(phi, section, inflow, rotor)

#     return R
# end

# function residual_outputs(B, af, turbine, inputs::AbstractArray)
#     phi = inputs[1]
#     r = inputs[2]
#     chord = inputs[3]
#     theta = inputs[4]
#     Vx, Vy, rho, mu, asound = inputs[5], inputs[6], inputs[7], inputs[8], inputs[9]
#     Rhub, Rtip, precone = inputs[10], inputs[11], inputs[12]

#     section = Section(r, chord, theta, af)
#     inflow = Inflow(Vx, Vy, rho, mu, asound)
#     rotor = Rotor(Rhub, Rtip, B, turbine, precone)

#     R, outputs = residual(phi, section, inflow, rotor)

#     return [getfield(outputs, i) for i in fieldnames(Outputs)]
# end

# function residual_partials(phi, section, inflow, rotor)
#     # unpack inputs
#     @unpack r, chord, theta, af = section
#     @unpack Vx, Vy, rho, mu, asound = inflow
#     @unpack Rhub, Rtip, B, turbine, precone = rotor

#     # Get a version of the residual function that's compatible with ForwardDiff.
#     function R(inputs)
#         return residual(B, af, turbine, inputs)
#     end

#     # Do it.
#     x = [phi, r, chord, theta, Vx, Vy, rho, mu, asound, Rhub, Rtip, precone]
#     return ResidualPartials(ForwardDiff.gradient(R, x)...)

# end

# function output_partials(phi, section, inflow, rotor)
#     # unpack inputs
#     @unpack r, chord, theta, af = section
#     @unpack Vx, Vy, rho, mu, asound = inflow
#     @unpack Rhub, Rtip, B, turbine, precone = rotor

#     # Get a version of the residual function that's compatible with ForwardDiff.
#     function R(inputs)
#         return residual_outputs(B, af, turbine, inputs)
#     end

#     # Do it.
#     x = [phi, r, chord, theta, Vx, Vy, rho, mu, asound, Rhub, Rtip, precone]
#     return ForwardDiff.jacobian(R, x)
# end

function pack_inputs(phi, rotor, inflow, idx)
    r = rotor.r[idx]
    chord = rotor.chord[idx]
    theta = rotor.theta[idx]
    Rhub = rotor.Rhub
    Rtip = rotor.Rtip
    precone = rotor.precone
    Vx = inflow.Vx[idx]
    Vy = inflow.Vy[idx]
    pitch = inflow.pitch
    rho = inflow.rho
    mu = inflow.mu
    asound = inflow.asound
    return [phi, r, chord, theta, Rhub, Rtip, precone, Vx, Vy, pitch, rho, mu, asound]
end

function unpack_inputs(x)
    i = 1
    phi = x[i]
    i += 1
    r = x[i]
    i += 1
    chord = x[i]
    i += 1
    theta = x[i]
    i += 1
    Rhub = x[i]
    i += 1
    Rtip = x[i]
    i += 1
    precone = x[i]
    i += 1
    Vx = x[i]
    i += 1
    Vy = x[i]
    i += 1
    pitch = x[i]
    i += 1
    rho = x[i]
    i += 1
    mu = x[i]
    i += 1
    asound = x[i]

    kwargs = Dict(
        :phi => phi,
        :r => r,
        :chord => chord,
        :theta => theta,
        :Rhub => Rhub,
        :Rtip => Rtip,
        :precone => precone,
        :Vx => Vx,
        :Vy => Vy,
        :pitch => pitch,
        :rho => rho,
        :mu => mu,
        :asound => asound)
    return kwargs
end

function pack_outputs(outputs, idx)
    Np = outputs.Np[idx]
    Tp = outputs.Tp[idx]
    a = outputs.a[idx]
    ap = outputs.ap[idx]
    u = outputs.u[idx]
    v = outputs.v[idx]
    alpha = outputs.alpha[idx]
    W = outputs.W[idx]
    cl = outputs.cl[idx]
    cd = outputs.cd[idx]
    cn = outputs.cn[idx]
    ct = outputs.ct[idx]
    F = outputs.F[idx]
    G = outputs.G[idx]
    return [Np, Tp, a, ap, u, v, alpha, W, cl, cd, cn, ct, F, G]
end

function unpack_jacobian(x::AbstractArray{T, N}) where {T, N}
    out = Dict{Tuple{Symbol, Symbol}, T}()
    i = 1
    Np = unpack_inputs(x[i, :])
    for (key, val) in Np
        out[(:Np, key)] = val
    end
    i += 1
    Tp = unpack_inputs(x[i, :])
    for (key, val) in Tp
        out[(:Tp, key)] = val
    end
    i += 1
    a = unpack_inputs(x[i, :])
    for (key, val) in a
        out[(:a, key)] = val
    end
    i += 1
    ap = unpack_inputs(x[i, :])
    for (key, val) in ap
        out[(:ap, key)] = val
    end
    i += 1
    u = unpack_inputs(x[i, :])
    for (key, val) in u
        out[(:u, key)] = val
    end
    i += 1
    v = unpack_inputs(x[i, :])
    for (key, val) in v
        out[(:v, key)] = val
    end
    i += 1
    alpha = unpack_inputs(x[i, :])
    for (key, val) in alpha
        out[(:alpha, key)] = val
    end
    i += 1
    W = unpack_inputs(x[i, :])
    for (key, val) in W
        out[(:W, key)] = val
    end
    i += 1
    cl = unpack_inputs(x[i, :])
    for (key, val) in cl
        out[(:cl, key)] = val
    end
    i += 1
    cd = unpack_inputs(x[i, :])
    for (key, val) in cd
        out[(:cd, key)] = val
    end
    i += 1
    cn = unpack_inputs(x[i, :])
    for (key, val) in cn
        out[(:cn, key)] = val
    end
    i += 1
    ct = unpack_inputs(x[i, :])
    for (key, val) in ct
        out[(:ct, key)] = val
    end
    i += 1
    F = unpack_inputs(x[i, :])
    for (key, val) in F
        out[(:F, key)] = val
    end
    i += 1
    G = unpack_inputs(x[i, :])
    for (key, val) in G
        out[(:G, key)] = val
    end
    return out
end

function residual_partials(phi, rotor, inflow, outputs, idx)

    B = rotor.B
    af = rotor.af[idx]
    turbine = rotor.turbine

    function R(_x)
        kwargs = unpack_inputs(_x)
        _phi = kwargs[:phi]
        _rotor = Rotor(; B=B, af=af, turbine=turbine, kwargs...)
        _inflow = OperatingPoint(; kwargs...)
        _outputs = Outputs(eltype(_x), 1)
        return residual!(_phi, _rotor, _inflow, _outputs, 1, false)
    end

    x = pack_inputs(phi, rotor, inflow, idx)
    return unpack_inputs(ForwardDiff.gradient(R, x))
end

function output_partials(phi, rotor, inflow, outputs, idx)
    B = rotor.B
    af = rotor.af[idx]
    turbine = rotor.turbine

    function out(_x)
        kwargs = unpack_inputs(_x)
        _phi = kwargs[:phi]
        _rotor = Rotor(; B=B, af=af, turbine=turbine, kwargs...)
        _inflow = OperatingPoint(; kwargs...)
        _outputs = Outputs(eltype(_x), 1)
        _R = residual!(_phi, _rotor, _inflow, _outputs, 1, true)
        return pack_outputs(_outputs, 1)
    end

    x = pack_inputs(phi, rotor, inflow, idx)
    return unpack_jacobian(ForwardDiff.jacobian(out, x))
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::Rotor, op::OperatingPoint)

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

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

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

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

    # unpack
    setsign = rotor.turbine ? 1 : -1
    Vx = op.Vx * setsign  # TODO: ideally should do sign swap only in one place
    Vy = op.Vy
    # theta = rotor.theta * setsign
    # pitch = op.pitch * setsign

    # ---- 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

    nr = length(rotor.r)

    # initialize outputs
    outputs = Outputs(eltype(rotor.r[1]), nr)

    for i = 1:nr

        startfrom90 = false

        if Vx[i] > 0 && Vy[i] > 0
            order = (q1, q2, q3, q4)
        elseif Vx[i] < 0 && Vy[i] > 0
            order = (q2, q1, q4, q3)
        elseif Vx[i] > 0 && Vy[i] < 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, op, outputs, i)

        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

            # 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 = Roots.fzero(R, phiL, phiU)

                so = PyCall.pyimport("scipy.optimize")
                phistar = so.brentq(R, phiL, phiU)
                residual!(phistar, rotor, op, outputs, i, true)  # call once more to set outputs
                break
            end
        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

# function solve(rotor::Rotor, op::AbstractArray{OperatingPoint, 1})
#     n = length(op)
#     output_array = Array{Outputs, n}

#     for i = 1:n
#         output_array[i] = solve(rotor, op[i])
#     end
# end

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

"""
    simple_op(Vinf, Omega, r, rho, pitch=0.0, 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{Float64, 1}`: radial location where inflow is computed (m)
- `rho::Float`: air density (kg/m^3)
- `pitch::Float`: pitch (rad)
- `mu::Float`: air viscosity (Pa * s)
- `asounnd::Float`: air speed of sound (m/s)
- `precone::Float64`: precone angle (rad)
"""
function simple_op(Vinf, Omega, r, rho, pitch=0.0, mu=1.0, asound=1.0, precone=0.0)

    Vx = Vinf * cos(precone) * ones(length(r))
    Vy = Omega .* r * cos(precone)

    return OperatingPoint(Vx, Vy, pitch, 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)
- `pitch::Float64`: pitch angle (rad)
- `r::Array{Float64, 1}`: radial locations 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::Float`: air viscosity (Pa * s)
- `asounnd::Float`: air speed of sound (m/s)
"""
function windturbine_op(Vhub, Omega, pitch, 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, pitch, rho, mu, asound)

end

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

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

"""
    thrusttorque(rotor::Rotor, outputs::Outputs)

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

**Arguments**
- `rotor::Rotor`: rotor object
- `outputs::Outputs`: output data along blade

**Returns**
- `T::Array{Float64, 1}`: thrust (along x-dir see Documentation).
- `Q::Array{Float64, 1}`: torque (along x-dir see Documentation).
"""
function thrusttorque(rotor::Rotor, outputs::Outputs)

    # add hub/tip for complete integration.  loads go to zero at hub/tip.
    rfull = [rotor.Rhub; rotor.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 * trapz(rfull, thrust)
    Q = rotor.B * trapz(rfull, torque)

    return T, Q
end

function thrusttorque(B, inputs::AbstractArray)

    num_radial = div(length(inputs)-3, 3)

    r = inputs[1:num_radial]
    Np = inputs[num_radial+1:2*num_radial]
    Tp = inputs[2*num_radial+1:3*num_radial]

    Rhub = inputs[3*num_radial+1]
    Rtip = inputs[3*num_radial+2]
    precone = inputs[3*num_radial+3]

    # Dummy values for stuff we don't need for the integrated loads. I guess I
    # could just reuse r, but whatever.
    array_dummy = similar(r)
    af_dummy = 0.

    rotor = Rotor(Rhub, Rtip, B, false, precone)
    sections = Section.(r, array_dummy, array_dummy, af_dummy)
    outputs = Outputs.(Np, Tp,
                       array_dummy,
                       array_dummy,
                       array_dummy,
                       array_dummy,
                       array_dummy,
                       array_dummy,
                       array_dummy,
                       array_dummy,
                       array_dummy)
    T, Q = thrusttorque(rotor, sections, outputs)

    return [T, Q]

end

function thrusttorque_partials(rotor, sections, Np, Tp)

    B = rotor.B
    function R(inputs)
        return thrusttorque(B, inputs)
    end

    num_radial = length(sections)
    x = zeros(typeof(sections[1].r), 3*num_radial+3)
    @. x[1:num_radial] = getfield(sections, :r)
    @. x[num_radial+1:2*num_radial] = Np
    @. x[2*num_radial+1:3*num_radial] = Tp
    x[3*num_radial+1] = rotor.Rhub
    x[3*num_radial+2] = rotor.Rtip
    x[3*num_radial+3] = rotor.precone

    return ForwardDiff.jacobian(R, x)
end

"""
    thrusttorque(rotor::Rotor, outputs::AbstractArray{Outputs{TF}, 1}) where TF <: Number

Integrate the thrust/torque across the blade given an array of output data.
Generally used for azimuthal averaging of thrust/torque.
"""
function thrusttorque(rotor::Rotor, outputs::AbstractArray{Outputs{TF}, 1}) where TF <: Number

    T = 0.0
    Q = 0.0
    n = length(outputs)

    for i = 1:n
        Tsub, Qsub = thrusttorque(rotor, outputs[i])
        T += Tsub / n
        Q += Qsub / n
    end

    return T, Q
end

"""
(private) trapezoidal integration
"""
function trapz(x, y)  # integrate y w.r.t. x

    integral = 0.0
    for i = 1:length(x)-1
        integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1])
    end
    return integral
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:
"""
    friedman_test_statistic(ranks::Vector,n_datasets::Int,n_models::Int)

Friedman test - do all the models perform the same?
"""
friedman_test_statistic(R::Vector,n::Int,k::Int) = 12*n/(k*(k+1))*(sum(R.^2) - k*(k+1)^2/4)

"""
    friedman_critval(α::Real, k::Int)

Critical value of the Friedman test for k models.
"""
friedman_critval(α::Real, k::Int) = crit_chisq(α/2, k-1)
crit_chisq(α::Real, df::Int) = quantile(Chisq(df), 1-α)

"""
    nemenyi_cd(n_models::Int, n_datasets::Int, α::Real)

Nemenyi critical difference value. If average ranks of two models are larger than this,
their performance difference is statistically siginificant.
"""
nemenyi_cd(k::Int, n::Int, α::Real) = sqrt(k*(k+1)/(6*n))*crit_srd(α, k, Inf)/sqrt(2)
crit_srd(α::Real, k::Real, df::Real) =
    (isnan(k) | isnan(df)) ? NaN : quantile(StudentizedRange(df, k), 1-α)

"""
    string2file(f, s)
Save string s to file f.
"""
function string2file(f, s)
    open(f, "w") do _f
        write(_f, s)
    end
end

"""
    eol(s)
Replaces the "& " at the end of s with a tabular end of line.
"""
function eol(s)
    return string(s[1:end-2], " \\\\ \n")
end

"""
    wspad(s, n)

Pads s with n white spaces.
"""
function wspad(s, n)
    return string(s, repeat(" ", n))
end

"""
    ranks2tikzcd(ranks, algnames, nemenyi_cd, caption = ""; scale=1.0,
        vpad=0.3, label = "", pos = "h", levelwidth=0.06, vdist=0.1)

Draws a critical difference diagram in tikz.
"""
function ranks2tikzcd(ranks, algnames, c, caption = ""; scale=1.0,
    vpad=0.3, label = "", pos = "h", levelwidth=0.06, vdist=0.1)
    n = length(ranks)
    @assert n == length(algnames)
    ranks = reshape(ranks, n)

    # header
    s = "\\begin{tikzpicture}[scale=$(scale)] \n"
    s = wspad(s, 2)
    #axis
    s = string(s, drawaxis(ranks))
    # nodes
    s = string(s, drawnodes(ranks, algnames; vpad=vpad))
    # levels
    s = string(s, drawlevels(ranks, c; levelwidth=levelwidth, vdist=vdist))
    # end of figure
    s = wspad(s, 1)
    s = string(s, "\\end{tikzpicture} \n")

    return s
end

function drawaxis(ranks)
    mx = ceil(maximum(ranks))
    mn = floor(minimum(ranks))

    s = "\\draw ($(mn),0) -- ($(mx),0); \n"
    s = wspad(s, 2)
    s = string(s, "\\foreach \\x in {$(Int(mn)),...,$(Int(mx))} \\draw (\\x,0.10) -- (\\x,-0.10) node[anchor=north]{\$\\x\$}; \n")
    return s
end

function drawnodes(ranks, algnames; vpad=0.3)
    isort = sortperm(ranks)
    ranks = ranks[isort]
    algnames = algnames[isort]

    mx = ceil(ranks[end])
    mn = floor(ranks[1])
    nor = length(ranks) # number of ranks
    nos = Int(ceil(nor/2)) # number of ranks on one side (start with left, there may be one more if nor is odd)

    s = ""
    for (i, r) in enumerate(ranks)
        s = wspad(s,2)
        if i<=nos
            s = string(s, "\\draw ($(r),0) -- ($r,$(i*vpad-0.1)) -- ($(mn-0.1), $(i*vpad-0.1)) node[anchor=east] {$(algnames[i])}; \n")
        else
            s = string(s, "\\draw ($(r),0) -- ($r,$((nor-i)*vpad+0.2)) -- ($(mx+0.1), $((nor-i)*vpad+0.2)) node[anchor=west] {$(algnames[i])}; \n")
        end
    end

    return s
end

function drawlevels(ranks, cv; levelwidth=0.06, vdist=0.1)
    ranks = sort(ranks)

    mx = ceil(ranks[end])
    mn = floor(ranks[1])
    nor = length(ranks) # number of ranks

    # do this differently
    ranks = sort(ranks)

    mx = ceil(ranks[end])
    mn = floor(ranks[1])
    nor = length(ranks) # number of ranks

    # do a matrix that says which ranks are connected
    con_mat = zeros(nor, nor)
    for i in 1:nor
        for j in 1:nor
            con_mat[i,j] = abs(ranks[j] - ranks[i]) < cv
        end
    end

    # now
    levels = [] # each element is a tuple containing the (start, end, height) values
    lastlasti = 0
    for i in 1:nor
        if length(levels) == 0
            global levelh = vdist
        elseif lastlasti < i
            global levelh = vdist
        else
            lastlevelh = levels[end][end]
            if lastlevelh == vdist
                global levelh = lastlevelh + vdist
            else
                global levelh = lastlevelh + vdist
            end
        end
        startl = ranks[i]
        lasti = findlast(con_mat[i,:] .== 1)
        if lasti != nothing && lasti != lastlasti && lasti != i
            endl = ranks[lasti]
            push!(levels, (startl, endl, levelh))
        end
        lastlasti = lasti
    end

    # draw them
    s = ""
    for level in levels
        s = wspad(s,2)
        s = string(s, "\\draw[line width=$(levelwidth)cm,color=black,draw opacity=1.0] ($(level[1]-0.03),$(level[3])) -- ($(level[2]+0.03),$(level[3])); \n")
    end

    return s
end

ext: jl


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

"""
    AbstractJoiner

A method for joining spatial objects.
"""
abstract type AbstractJoiner end

"""
    join(object₁, object₂, joiner)

Join spatial object `object₁` and `object₂` with `joiner` method.
"""
Base.join(object₁::AbstractSpatialObject,
          object₂::AbstractSpatialObject,
          joiner::AbstractJoiner) =
  @error "not implemented"

#------------------
# IMPLEMENTATIONS
#------------------
include("joining/variable.jl")


text:
module Polynomials

##############################################################################
##
## TODO
##
##############################################################################

# - Should remove module and rely on package-level module as per convention
# - Methods e.g., promote, convert, macro +, -, *, /

##############################################################################
##
## Exports
##
##############################################################################

export Polynomial, degree

##############################################################################
##
## Implementation
##
##############################################################################

import Base.show,
       Base.length,
       Base.convert,
       Base.divrem,
       Base.promote_rule,
       Base.promote,
       Base.abs,
       Base.div,
       Base.rem,
       Base.zero,
       Base.zeros

# TODO: Rethink this subclassing; polynomial is a field element, not a number.
type Polynomial{T} <: Number
  coefficients::Array{T}

  function Polynomial()
    return new([0])
  end

  function Polynomial(c::Array)
    i = length(c)
    while i > 0
      if c[i] != 0
        break
      end
      i -= 1
    end

    a = Array(T, i)
    for j in 1:i
      a[j] = c[j]
    end

    new(a)
  end
end

convert{T<:Number}(::Type{Polynomial{T}}, x::T) = Polynomial(x)
convert{T<:Number, S<:Number}(::Type{Polynomial{T}}, x::S) = Polynomial(convert(T, x))
convert{T<:Number, S<:Number}(::Type{Polynomial{T}}, x::Polynomial{S}) = Polynomial([convert(T,y) for y in x.coefficients])

promote_rule{T<:Number}(::Type{Polynomial{T}}, ::Type{T}) = Polynomial
promote_rule{T<:Number, S<:Number}(::Type{Polynomial{T}}, ::Type{S}) = Polynomial{promote_type(T,S)}
promote_rule{T<:Number, S<:Number}(::Type{Polynomial{T}}, ::Type{Polynomial{S}}) = Polynomial{promote_type(T,S)}

function Polynomial{T}(c::Array{T})
  return Polynomial{T}([c])
end

function zero{T}(p::Polynomial{T})
  return Polynomial{T}([])
end

function zero{T}(a::Type{Polynomial{T}})
  return Polynomial{T}([])
end

function isZero(p::Polynomial)
  return p.coefficients == []
end

function abs(p::Polynomial)
  return length(p.coefficients)
end

function length(p::Polynomial)
  return length(p.coefficients)
end

function degree(p::Polynomial)
  return length(p.coefficients) - 1
end

function leadingCoefficient(p::Polynomial)
  return p.coefficients[length(p.coefficients)]
end

function show(io::IO, p::Polynomial)
  if isZero(p)
    print(io, "0")
  else
    print(io, string(p.coefficients[1]))
    for i in 2:length(p)
      if p.coefficients[i] != 0
        print(io, string(" + ", p.coefficients[i], "x^", i - 1))
      end
    end
  end
end

function ==(a::Polynomial, b::Polynomial)
  length(a) == length(b) && (isZero(a ) || !(false in [a.coefficients[i] == b.coefficients[i] for i in 1:length(a)]))
end

function -(p::Polynomial)
  return Polynomial([-i for i in p.coefficients])
end

function -(a::Polynomial, b::Polynomial)
  return a+(-b)
end

function +{T}(a::Polynomial{T}, b::Polynomial{T})
  newcoefficients = zeros(T, max(length(a), length(b)))
  i = 1
  while i <= min(length(a), length(b))
    newcoefficients[i] += (a.coefficients[i] + b.coefficients[i])
    i += 1
  end
  while i <= length(a)
    newcoefficients[i] += a.coefficients[i]
    i += 1
  end
  while i <= length(b)
    newcoefficients[i] += b.coefficients[i]
    i += 1
  end
  return Polynomial{T}(newcoefficients)
end

function *{T}(a::Polynomial{T}, b::Polynomial{T})
  if isZero(a) || isZero(b)
    return Polynomial{T}(Array(T, 0))
  else
    newcoefficients = zeros(T, length(a) + length(b) - 1)
    for i in 0:length(a)-1
      for j in 0:length(b)-1
        newcoefficients[i + j + 1] += a.coefficients[i + 1] * b.coefficients[j + 1]
      end
    end
    return Polynomial{T}(newcoefficients)
  end
end

function divrem{T}(a::Polynomial{T}, b::Polynomial{T})
  quotient = Polynomial{T}(Array(T, 0))
  remainder = a
  divisorDeg = degree(b)
  divisorLC = leadingCoefficient(b)

  while degree(remainder) >= divisorDeg
    monomialExponent = degree(remainder) - divisorDeg
    monomialZeros = zeros(T, monomialExponent)
    monomialDivisor = Polynomial{T}([monomialZeros, [convert(T, leadingCoefficient(remainder)/divisorLC)]])

    quotient += monomialDivisor
    remainder -= (monomialDivisor*b)
  end

  return (quotient, remainder)
end

function rem{T}(a::Polynomial{T}, b::Polynomial{T})
  q,r = divrem(a,b)
  return r
end

function /{T}(a::Polynomial{T}, b::Polynomial{T})
  q,r = divrem(a,b)
  return q
end

function ^{T}(a::Polynomial{T}, n::Integer)
  ret = Polynomial{T}([1])
  for i in 1:n
    ret = ret*a
  end
  return ret
end

end # module Polynomial


text: # parameters that govern counterfactual


text:
type PossibleExtremaReached
    x0::Real
end

type Dithering
    x0::Real
end

"""

SOLVE button from the HP-34C
follows roughly algorithm described http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf
though some modifications were made.

Goal is to return a value `x` with either:
* `f(x) == 0.0` or
* `f(prevfloat(x)) * f(nextfloat(x)) < 0`.
if a bracket is found that can be done.

Otherwise, values with `f(x) <= ftol` are given status of answers
and others are marked as a possible extrema.

If the algorithm gets stuck, a Dithering error is thrown.

While this is slower than fzero with order 2, 5, 8, or 16 it is
more robust to the initial condition.

"""
function SOLVE{T<:AbstractFloat}(f, x0::T;
               ftol::Real=10*eps(one(x0)),
               xtol::Real=zero(x0),
               xtolrel::Real=zero(x0),
               maxeval::Int=30,
               verbose::Bool=false)

    fx0 = f(x0)
    abs(fx0) <= ftol && return(x0)

    try
        ## get initial guess. Use Steffensen guess if reasonable size
        stepsize = max(1/100, min(abs(fx0), abs(x0/100)))
        if fx0 < 0
            g(x) = -f(x)
            secant_method_no_bracket(g, x0 - stepsize, x0, ftol=ftol, xtol=xtol, xtolrel=xtolrel, maxeval=maxeval, verbose=verbose)
        else
            secant_method_no_bracket(f, x0 - stepsize, x0, ftol=ftol, xtol=xtol, xtolrel=xtolrel, maxeval=maxeval, verbose=verbose)
        end
    catch ex
        if isa(ex, StateConverged)
            return convert(T,ex.x0)
        elseif isa(ex, PossibleExtremaReached)
            guess = convert(T, ex.x0)
            if sign(f(nextfloat(guess))) * sign(f(prevfloat(guess))) <= 0
                return guess
            end
            # heuristic for case where it thinks it is at a minimum
            if abs(f(guess)) < ftol
                return guess
            else
                rethrow(ex)
            end
        else
            rethrow(ex)
        end
    end
    error("huh? Shouldn't get here")
end

"""

Solve f(x) = 0 if we have a bracket [a,b] and starting point a < c < b

Use Algorithm 4.2 of Alefeld, Potra, Shi unless it dithers, in which case the root is found by bisection.

"""
function have_bracket(f, a, b, c=(0.5)*(a+b);
                      xtol::Real=zero(c), xtolrel::Real=zero(c), maxeval=10, verbose=false)
    ## try a42a unless it fails, then go to failsafe
    verbose && println("Have a bracket, switching to bracketing method")
    fa, fb = f(a), f(b)
    fa == 0.0 && return(a)
    fb == 0.0 && return(b)
    isinf(fa) && throw(ConvergenceFailed("f(a) must be finite"))
    isinf(fb) && throw(ConvergenceFailed("f(b) must be finite"))

    if sign(fa) * sign(fb) > 0
        throw(ConvergenceFailed("Interval [$a, $b] is not a bracket"))
    end

    a,b = sort([a,b])
    if !(a <= c <= b)
        throw(ConvergenceFailed("Value $c is not in interval ($a, $b)"))
    end

    try
        x0 = a42a(f, a, b, c, xtol=xtol, maxeval=maxeval, verbose=verbose)
    catch e
        if isa(e, StateConverged)
            return e.x0
        else
            verbose && println("Dithering, switching to bisection method")
            return find_zero(f, a, b, verbose=verbose)
        end
    end
end

## assume f(b) > 0
## steffensen of secant line method, but we don't assume a bracket
## we "bend" large trajectories
function secant_method_no_bracket{T<:AbstractFloat}(f, a::T, b::T;
                                  ftol=10*eps(one(float(a))), xtol::Real=10*eps(one(float(a))),
                                  xtolrel::Real=10*eps(one(float(a))),
                                  maxeval::Int=100, verbose::Bool=false)

    alpha, beta = a, b ## beta keeps smallest values
    falpha, fbeta = f(alpha), f(beta)

    ## checks -- have a bracket
    if sign(falpha) * sign(fbeta) < 0
        a = a42(f, alpha, beta, xtol=xtol, maxeval=maxeval, verbose=verbose)
        throw(StateConverged(a))
    end

    ## keep fbeta smallest
    if falpha < fbeta
        alpha, beta, falpha, fbeta = beta, alpha, fbeta, falpha
    end

    gamma = Inf

    ctr, quad_ctr = 0, 0

    while fbeta > 0
        if ctr > maxeval
            if abs(f(beta)) < ftol
                throw(StateConverged(beta))
            else
                ## try a higher order method? May get us there
                o = thukral8_itr(f, beta, xtol=float(xtol), xtolrel=float(xtolrel), ftol = float(ftol))
                for (x,fx) in o
                    nothing
                end
                xn = o.xn[end]
                if o.state == :converged && !isnan(xn) && !isinf(xn)
                    throw(StateConverged(xn))
                else
                    throw(Dithering(beta))
                end
            end
        else
            (ctr = ctr + 1)
        end
        ## use steffensen or secant line depending on fbeta's size
        if abs(fbeta) < abs(beta - alpha)
            stepsize = fbeta * fbeta / (f(beta + fbeta) - fbeta)  # steffensen
        else
            stepsize =  fbeta * (beta - alpha)  / (fbeta - falpha) # secant line
        end
        gamma = beta - stepsize

        ## bend in some cases. XXX What is best bending scheme? XXX
        if abs(gamma-beta) >= 5e2 * abs(beta - alpha)
            stepsize = cbrt( (beta - alpha) / (fbeta - falpha) * fbeta)
#            gamma = beta + (1/4) * (beta - alpha) / (fbeta - falpha) * fbeta
            gamma = beta - stepsize
        end

        fgamma = f(gamma)

        # try
        #     fgamma = f(gamma) ## what to do with domain error?
        # catch e
        #     if isa(e, DomainError)
        #         error("domain error -- fix me")
        #     end
        # end

        if abs(fgamma) <= ftol
            throw(StateConverged(gamma))
        end
        if fgamma < 0
            x0 = have_bracket(f, gamma, beta, maxeval=20, verbose=verbose)
            throw(StateConverged(x0))
        end

        ## increasing value, try quadratic
        if fgamma >= fbeta
            quad_ctr = quad_ctr + 1
            if quad_ctr > 5*10
                abs(fbeta) <  ftol ? throw(StateConverged(beta)) : throw(Dithering(beta))
            end
            ## secant line stopped decreasing
            if alpha == gamma || gamma == beta
                throw(PossibleExtremaReached(gamma)) ## poorly titled error
            end

            ## quadratic interpolation to get gamma
            denom = (beta - alpha) * (fbeta - fgamma)  - (beta - fgamma) * (fbeta - falpha)
            if denom == 0.0
                throw(PossibleExtremaReached(alpha))
            end
            gamma = beta -  ((beta - alpha)^2 * (fbeta - fgamma) - (beta - gamma)^2 * (fbeta - falpha))/denom/2

            fgamma = f(gamma)

            if fgamma == 0
                throw(StateConverged(gamma))
            end
            if fgamma < 0
                x0 = have_bracket(f, gamma, beta, maxeval=20, verbose=verbose)
                throw(StateConverged(x0))
            end
        else
            quad_ctr = 0        # reset
        end

        if fgamma > fbeta
            alpha, beta, falpha, fbeta = gamma, beta, fgamma, fbeta
        elseif fgamma < fbeta
            alpha, beta, falpha, fbeta = beta, gamma, fbeta, fgamma
        else
            throw(PossibleExtremaReached(beta))
        end
        verbose && println(@sprintf("x_%s:%18.15f\t f(x_%s):%18.15f", ctr, beta, ctr, f(beta)))
    end
    if abs(f(beta)) < ftol
        throw(StateConverged(beta))
    else
        throw(ConvergenceFailed("More than $maxeval iterations before convergence"))
    end
end


text:
cell_to_vtkcell(::Type{Line}) = VTKCellTypes.VTK_LINE
cell_to_vtkcell(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE

cell_to_vtkcell(::Type{Quadrilateral}) = VTKCellTypes.VTK_QUAD
cell_to_vtkcell(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD
cell_to_vtkcell(::Type{Triangle}) = VTKCellTypes.VTK_TRIANGLE
cell_to_vtkcell(::Type{QuadraticTriangle}) = VTKCellTypes.VTK_QUADRATIC_TRIANGLE
cell_to_vtkcell(::Type{Cell{2,8,4}}) = VTKCellTypes.VTK_QUADRATIC_QUAD

cell_to_vtkcell(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON
cell_to_vtkcell(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA
cell_to_vtkcell(::Type{QuadraticTetrahedron}) = VTKCellTypes.VTK_QUADRATIC_TETRA

"""
    vtk_grid(filename::AbstractString, grid::Grid)

Create a unstructured VTK grid from a `Grid`. Return a `DatasetFile`
which data can be appended to, see `vtk_point_data` and `vtk_cell_data`.
"""
function WriteVTK.vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; compress::Bool=true) where {dim,C,T}
    cls = MeshCell[]
    for cell in grid.cells
        celltype = Ferrite.cell_to_vtkcell(typeof(cell))
        push!(cls, MeshCell(celltype, collect(cell.nodes)))
    end
    coords = reshape(reinterpret(T, getnodes(grid)), (dim, getnnodes(grid)))
    return vtk_grid(filename, coords, cls; compress=compress)
end

"""
    vtk_point_data(vtk, data::Vector{<:Tensor}, name)

Write the tensor field data to the vtk file. Only writes the tensor values available in `data`.
In the vtu-file, ordering of the tensor components is coulumn-wise (just like Julia).
[1 2
 3 4] => 1, 3, 2, 4
"""
function WriteVTK.vtk_point_data(
    vtk::WriteVTK.DatasetFile,
    data::Union{
        Vector{Tensor{order,dim,T,M}},
        Vector{SymmetricTensor{order,dim,T,M}}
        },
    name::AbstractString
    ) where {order,dim,T,M}

    npoints = length(data)
    out = zeros(T, M, npoints)
    out[1:M, :] .= reshape(reinterpret(T, data), (M, npoints))
    return vtk_point_data(vtk, out, name)
end

"""
    vtk_point_data(vtk, data::Vector{<:Vec}, name)

Write the vector field data to the vtk file.
"""
function WriteVTK.vtk_point_data(vtk::WriteVTK.DatasetFile, data::Vector{Vec{dim,T}}, name::AbstractString) where {dim,T}
    npoints = length(data)
    out = zeros(T, (dim == 2 ? 3 : dim), npoints)
    out[1:dim, :] .= reshape(reinterpret(T, data), (dim, npoints))
    return vtk_point_data(vtk, out, name)
end

function vtk_nodeset(vtk::WriteVTK.DatasetFile, grid::Grid{dim}, nodeset::String) where {dim}
    z = zeros(getnnodes(grid))
    z[collect(getnodeset(grid, nodeset))] .= 1.0
    vtk_point_data(vtk, z, nodeset)
end

"""
    vtk_cellset(vtk, grid::Grid)

Export all cell sets in the grid. Each cell set is exported with
`vtk_cell_data` with value 1 if the cell is in the set, and 0 otherwise.
"""
function vtk_cellset(vtk::WriteVTK.DatasetFile, grid::AbstractGrid, cellsets=keys(grid.cellsets))
    z = zeros(getncells(grid))
    for cellset in cellsets
        z .= 0.0
        z[collect(getcellset(grid, cellset))] .= 1.0
        vtk_cell_data(vtk, z, cellset)
    end
    return vtk
end

"""
    vtk_cellset(vtk, grid::Grid, cellset::String)

Export the cell set specified by `cellset` as cell data with value 1 if
the cell is in the set and 0 otherwise.
"""
vtk_cellset(vtk::WriteVTK.DatasetFile, grid::AbstractGrid, cellset::String) =
    vtk_cellset(vtk, grid, [cellset])

import Ferrite.field_offset
function WriteVTK.vtk_point_data(vtkfile, dh::MixedDofHandler, u::Vector, suffix="")

    fieldnames = Ferrite.getfieldnames(dh)  # all primary fields

    for name in fieldnames
        @debug println("exporting field $(name)")
        field_dim = getfielddim(dh, name)
        space_dim = field_dim == 2 ? 3 : field_dim
        data = fill(NaN, space_dim, getnnodes(dh.grid))  # set default value

        for fh in dh.fieldhandlers
            # check if this fh contains this field, otherwise continue to the next
            field_pos = findfirst(i->i == name, getfieldnames(fh))
            if field_pos == 0 && continue end

            cellnumbers = sort(collect(fh.cellset))  # TODO necessary to have them ordered?
            offset = field_offset(fh, name)

            for cellnum in cellnumbers
                cell = dh.grid.cells[cellnum]
                n = ndofs_per_cell(dh, cellnum)
                eldofs = zeros(Int, n)
                _celldofs = celldofs!(eldofs, dh, cellnum)
                counter = 1

                for node in cell.nodes
                    for d in 1:field_dim
                        data[d, node] = u[_celldofs[counter + offset]]
                        @debug println("  exporting $(u[_celldofs[counter + offset]]) for dof#$(_celldofs[counter + offset])")
                        counter += 1
                    end
                    if field_dim == 2
                        # paraview requires 3D-data so pad with zero
                        data[3, node] = 0
                    end
                end
            end
        end
        vtk_point_data(vtkfile, data, string(name, suffix))
    end

    return vtkfile
end


text:
"""
    AbstractRF

Abstract type for representing radiofrequency (RF) pulses.
"""
abstract type AbstractRF end

"""
    InstantaneousRF(α, θ = 0) <: AbstractRF

Represents an idealized instantaneous RF pulse with flip angle `α` and phase
`θ`.
"""
struct InstantaneousRF{T<:Real} <: AbstractRF
    α::T
    θ::T
end

InstantaneousRF(α, θ = zero(α)) = InstantaneousRF(promote(α, θ)...)

Base.show(io::IO, rf::InstantaneousRF) = print(io, "InstantaneousRF(", rf.α, ", ", rf.θ, ")")
Base.show(io::IO, ::MIME"text/plain", rf::InstantaneousRF{T}) where {T} =
    print(io, "Instantaneous RF pulse with eltype $T:\n α = ", rf.α, " rad\n θ = ", rf.θ, " rad")

"""
    RF(waveform, Δt, [Δθ], [grad]) <: AbstractRF

Represents an RF pulse with the given (possibly complex-valued) waveform (G) and
time step `Δt` (ms). `Δθ` is additional phase added to the waveform (defaults to
`0`), and `grad` is the B0 gradient that is turned on during the RF pulse
(defaults to `Gradient(0, 0, 0)`, i.e., turned off).

# Properties
- `α::Vector{<:Real}`: Instantaneous flip angles (rad) at each time point;
  computed from the magnitude of `waveform`
- `θ::Vector{<:Real}`: Instantaneous phase (rad) at each time point; computed
  from the phase of `waveform`
- `Δt::Real`: Time step (ms)
- `Δθ_initial::Real`: Phase added to `θ` before any phase-cycling increment has
  been applied
- `Δθ::Ref{<:Real}`: Phase to be added to `θ`; can be updated to simulate
  phase-cycling/RF spoiling
- `grad`: Gradient applied during the RF pulse
  - `::Gradient`: Constant gradient
  - `::Vector{<:Gradient}`: Time-varying gradient
"""
struct RF{T<:Real,G<:Union{<:Gradient,<:AbstractVector{<:Gradient}}} <: AbstractRF
    α::Vector{T}
    θ::Vector{T}
    Δt::Float64
    Δθ_initial::T
    Δθ::Ref{T} # Type Ref to enable RF-spoiling, which requires updating Δθ
    grad::G

    function RF(
        α::AbstractVector{<:Real},
        θ::AbstractVector{<:Real},
        Δt::Real,
        Δθ::Real,
        grad::Union{<:Gradient,<:AbstractVector{<:Gradient}}
    )

        length(α) == length(θ) || error("α and θ must have the same number of elements")
        grad isa AbstractVector && (length(grad) == length(α) ||
            error("grad is a vector but has a different number of elements than α"))
        T = promote_type(eltype(α), eltype(θ), typeof(Δθ))
        new{T,typeof(grad)}(α, θ, Δt, Δθ, Δθ, grad)

    end
end

# waveform in Gauss, Δt in ms
RF(waveform, Δt, Δθ, grad) = RF(GAMMA .* abs.(waveform) .* (Δt / 1000), angle.(waveform), Δt, Δθ, grad)
RF(waveform, Δt, Δθ::Real) = RF(waveform, Δt, Δθ, Gradient(0, 0, 0))
RF(waveform, Δt, grad) = RF(waveform, Δt, 0, grad)
RF(waveform, Δt) = RF(waveform, Δt, 0, Gradient(0, 0, 0))

Base.show(io::IO, rf::RF) = print(io, "RF(", rf.α, ", ", rf.θ, ", ", rf.Δt, ", ", rf.Δθ_initial, ", ", rf.grad, ")")

function Base.show(io::IO, ::MIME"text/plain", rf::RF{T,G}) where {T,G}

    print(io, "RF{$T,$G}:")
    print(io, "\n α = ", rf.α, " rad")
    print(io, "\n θ = ", rf.θ, " rad")
    print(io, "\n Δt = ", rf.Δt, " ms")
    print(io, "\n Δθ (initial) = ", rf.Δθ_initial, " rad")
    print(io, "\n Δθ (current) = ", rf.Δθ[], " rad")
    print(io, "\n grad = ", rf.grad, " G/cm")

end

Base.length(rf::RF) = length(rf.α)

"""
    duration(rf)

Return the duration (ms) of the RF pulse.
"""
duration(::InstantaneousRF) = 0
duration(rf::RF) = length(rf) * rf.Δt

struct ExcitationWorkspace{T1,T2,T3,T4,T5}
    Af::T1
    Bf::T2
    Ae::T3
    tmpA1::T4
    tmpA2::T4
    tmpB1::T2
    tmpB2::T2
    bm_workspace::T5
end

function ExcitationWorkspace(
    spin::AbstractSpin,
    bm_workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

    ExcitationWorkspace(typeof(spin), bm_workspace)

end

function ExcitationWorkspace(
    spin::Union{Type{Spin{T}},Type{SpinMC{T,N}}},
    bm_workspace = spin <: Spin ? nothing : BlochMcConnellWorkspace(spin)
) where {T,N}

    Ae = ExcitationMatrix{T}()
    if spin <: Spin
        Af = FreePrecessionMatrix{T}()
        Bf = Magnetization{T}()
        tmpA1 = BlochMatrix{T}()
        tmpA2 = BlochMatrix{T}()
        tmpB1 = Magnetization{T}()
        tmpB2 = Magnetization{T}()
    else
        Af = BlochMcConnellMatrix{T}(N)
        Bf = MagnetizationMC{T}(N)
        tmpA1 = BlochMcConnellMatrix{T}(N)
        tmpA2 = BlochMcConnellMatrix{T}(N)
        tmpB1 = MagnetizationMC{T}(N)
        tmpB2 = MagnetizationMC{T}(N)
    end
    ExcitationWorkspace(Af, Bf, Ae, tmpA1, tmpA2, tmpB1, tmpB2, bm_workspace)

end

"""
    rotatetheta!(A, θ, α)

Simulate left-handed rotation by angle `α` about an axis in the x-y plane that
makes angle `θ` with the negative y-axis, overwriting `A`.
"""
function rotatetheta!(A, θ, α)

    (sinθ, cosθ) = sincos(θ)
    (sinα, cosα) = sincos(α)

    A.a11 = sinθ^2 + cosα * cosθ^2
    A.a21 = sinθ * cosθ - cosα * sinθ * cosθ
    A.a31 = -sinα * cosθ
    A.a12 = sinθ * cosθ - cosα * sinθ * cosθ
    A.a22 = cosθ^2 + cosα * sinθ^2
    A.a32 = sinα * sinθ
    A.a13 = sinα * cosθ
    A.a23 = -sinα * sinθ
    A.a33 = cosα

    return nothing

end

"""
    excite(spin, rf::InstantaneousRF, [nothing])
    excite(spin, rf::RF, [workspace])

Simulate excitation for the given spin. Returns `(A, B)` such that `A * M + B`
applies excitation to the magnetization `M`. If `isnothing(B)` (as is the case
for `InstantaneousRF`s), then `A * M` applies excitation to `M`.

For `RF` objects, `workspace isa ExcitationWorkspace`. For `SpinMC` objects, use
`workspace = ExcitationWorkspace(spin, nothing)` to use an approximate matrix
exponential to solve the Bloch-McConnell equation.

For an in-place version, see [`excite!`](@ref).

# Examples
```jldoctest
julia> s = Spin(1, 1000, 100, 3.75); s.M
Magnetization vector with eltype Float64:
 Mx = 0.0
 My = 0.0
 Mz = 1.0

julia> (A,) = excite(s, InstantaneousRF(π/2, π/4)); A * s.M
Magnetization vector with eltype Float64:
 Mx = 0.7071067811865476
 My = -0.7071067811865475
 Mz = 6.123233995736766e-17
```
"""
function excite(spin::AbstractSpin, rf::InstantaneousRF, ::Nothing = nothing)

    A = ExcitationMatrix{eltype(spin)}()
    excite!(A, spin, rf)
    return (A, nothing)

end

"""
    excite!(A, [nothing], spin, rf::InstantaneousRF, [nothing])
    excite!(A, B, spin, rf::RF, [workspace])

Simulate excitation, overwriting `A` and `B` (in-place version of
[`excite`](@ref)).
"""
function excite!(A::ExcitationMatrix, spin::AbstractSpin, rf::InstantaneousRF)

    rotatetheta!(A.A, rf.θ, rf.α)

end

function excite!(A::ExcitationMatrix, ::Nothing, spin::AbstractSpin, rf::InstantaneousRF, ::Nothing = nothing)

    excite!(A, spin, rf)

end

function excite(spin::AbstractSpin, rf::RF, workspace = ExcitationWorkspace(spin))

    T = eltype(spin)
    if spin isa Spin
        A = BlochMatrix{T}()
        B = Magnetization{T}()
    else
        N = spin.N
        A = BlochMcConnellMatrix{T}(N)
        B = MagnetizationMC{T}(N)
    end
    excite!(A, B, spin, rf, workspace)
    return (A, B)

end

function excite!(
    A::Union{<:BlochMatrix,<:BlochMcConnellMatrix},
    B::Union{<:Magnetization,<:MagnetizationMC},
    spin::AbstractSpin,
    rf::RF{T,G},
    workspace::ExcitationWorkspace = ExcitationWorkspace(spin)
) where {T,G}

    Δtby2 = rf.Δt / 2
    if G <: AbstractVector
        freeprecess!(workspace.Af, workspace.Bf, spin, Δtby2, rf.grad[1], workspace.bm_workspace)
    else
        freeprecess!(workspace.Af, workspace.Bf, spin, Δtby2, rf.grad, workspace.bm_workspace)
    end
    excite!(workspace.Ae, spin, InstantaneousRF(rf.α[1], rf.θ[1] + rf.Δθ[]))

    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Af, workspace.Bf, workspace.Ae, nothing)
    combine!(A, B, workspace.tmpA1, workspace.tmpB1, workspace.Af, workspace.Bf)

    for t = 2:length(rf)

        if G <: AbstractVector
            freeprecess!(workspace.Af, workspace.Bf, spin, Δtby2, rf.grad[t], workspace.bm_workspace)
        end
        excite!(workspace.Ae, spin, InstantaneousRF(rf.α[t], rf.θ[t] + rf.Δθ[]))

        combine!(workspace.tmpA1, workspace.tmpB1, A, B, workspace.Af, workspace.Bf)
        combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Ae, nothing)
        combine!(A, B, workspace.tmpA2, workspace.tmpB2, workspace.Af, workspace.Bf)

    end

end

"""
    excite!(spin, ...)

Apply excitation to the given spin, overwriting the spin's magnetization vector.
"""
function excite!(spin::AbstractSpin, args...)

    (A, B) = excite(spin, args...)
    BtoM = copy(spin.M)
    applydynamics!(spin, BtoM, A, B)

end


text:
export solve

"""
```julia
solve(nls::AbstractNLS,
      θ_init::AbstractVector,
      conf::Abstract_Solver_Conf)::Abstract_Solver_Result
```

Generic interface to solve an [`AbstractNLS`](@ref) problem.

The used algorithm is defined through [`Abstract_Solver_Conf`](@ref)
specializations.

The method returns a [`Abstract_Solver_Result`](@ref) specialization.
"""
solve(nls::AbstractNLS,
      θ_init::AbstractVector,
      conf::Abstract_Solver_Conf) = @assert(false,"To implement")


text:
using BitsFields
using Test

field1span  =  6
field1shift =  0
field2span  = 10
field2shift = field1span

field1 = BitField(UInt16, field1span, field1shift, :field1)
field2 = BitField(UInt16, field2span, field2shift, :field2)

bitfields = BitFields(field1, field2)

workingbits = zero(bitfields)

field1value = 0x15
field2value = 0x02f6

set!(bitfields[1], field1value, workingbits)
set!(bitfields[2], field2value, workingbits)

@test get(bitfields[2], workingbits) == UInt16(0x02f6)

@test get(bitfields, workingbits)    == [ UInt16(0x15), UInt16(0x02f6) ]

#                    span      shift      name         bits used
bitfield1 = BitField(  8,            0, :field1);   #      8
bitfield2 = BitField(  4,            8, :field2);   #     12
bitfield3 = BitField( 12,          4+8, :field3);   #     24
bitfield4 = BitField( 16,       12+4+8, :field4);   #     40
bitfield5 = BitField( 20,    16+12+4+8, :field5);   #     60
bitfield6 = BitField(  4, 20+16+12+4+8, :field6);   #     64

bitfields = BitFields(bitfield1, bitfield2, bitfield3,
                      bitfield4, bitfield5, bitfield6);

valfield1 = 0xae;
valfield2 = 0x06;
valfield3 = 0x05f3;
valfield4 = 0x3113;
valfield5 = 0x7654;
valfield6 = 0x03;

target = zero(bitfields)

set!(bitfields[1], valfield1, target)
set!(bitfields[2], valfield2, target)
set!(bitfields[3], valfield3, target)
set!(bitfields[4], valfield4, target)
set!(bitfields[5], valfield5, target)
set!(bitfields[6], valfield6, target)

@test get(bitfield1, target) == valfield1
@test get(bitfield2, target) == valfield2
@test get(bitfield3, target) == valfield3
@test get(bitfield4, target) == valfield4
@test get(bitfield5, target) == valfield5
@test get(bitfield6, target) == valfield6

fields = get(bitfields, target)

result = [ 0x00000000000000ae,
           0x0000000000000006,
           0x00000000000005f3,
           0x0000000000003113,
           0x0000000000007654,
           0x0000000000000003 ]

@test fields == result

@test fields[1] == valfield1
@test fields[2] == valfield2
@test fields[3] == valfield3
@test fields[4] == valfield4
@test fields[5] == valfield5
@test fields[6] == valfield6

namedfields = NamedTuple(bitfields)

@test namedfields.field1 == bitfields[1]
@test namedfields.field2 == bitfields[2]
@test namedfields.field3 == bitfields[3]
@test namedfields.field4 == bitfields[4]
@test namedfields.field5 == bitfields[5]
@test namedfields.field6 == bitfields[6]


text:
#Simulates svd for bigfloat matrices by converting to Float64/Complex128.
Base.LinAlg.svdfact!(A::Matrix{BigFloat}; thin=true) = svdfact!(convert.(Float64, A); thin=true)
Base.LinAlg.svdfact!(A::Array{Complex{BigFloat},2}; thin=true) = svdfact!(convert.(Complex128, A); thin=true)


text:
NIDmm_Misc_Dict = Dict()

"""
dmm = connect!(obj::T) where {T<:NIDmm}
"""
function connect!(obj::T) where {T<:NIDmm}
	obj_ptr = ViPSession(0)
	obj_Int = [UInt8.(collect(obj.address)); UInt8(0)] # terminate with NULL char
	status = niDMM_init(obj_Int, Int16(0), Int16(0), obj_ptr)
	obj.handle = obj_ptr[]
	obj.initialized = true
	return obj
end

# close
"""
status = disconnect!(obj::T) where {T<:NIDmm}
"""
function disconnect!(obj::T) where {T<:NIDmm}
	status = niDMM_close(obj.handle)
	obj.initialized = false
	return status
end

# IDN
function get_idn(obj::T) where {T<:NIDmm}
	chlist = [UInt8.(collect("")); UInt8(0)] # terminate with NULL char
	bufSize = Int32(256)
	res    = zeros(UInt8,bufSize)
	pres   = Ref(res)
	status = niDMM_GetAttributeViString(obj.handle,
		chlist, NIDMM_ATTR_INSTRUMENT_MODEL, bufSize, res)
	return (status, unsafe_string(pointer(res)))
end


text:
#=
fibonacci:
- Julia version: 1.4.1
- Author: Baur Urazalinov
- Date: 05-05-2020
=#

function fib(n)
    if n == 0
        return 0
    end
    if n == 1
        return 1
    end
    return fib(n-1)+fib(n-2)
end


text:
## tests of PhyloNetworklm

@testset "phyloNetworklm on small network" begin

tree_str= "(A:2.5,((B:1,#H1:0.5::0.1):1,(C:1,(D:0.5)#H1:0.5::0.9):1):0.5);"
net = readTopology(tree_str)
preorder!(net)
# printEdges(net)
# plot(net, useEdgeLength = true,  showEdgeLength=true, showGamma=true)

# Rk: Is there a way to check that the branch length are coherent with
# one another (Especialy for hybrids) ?
# Not yet (CA, 2016-12-01).
# Would be great to add functions to calculate distance node - root.
# several such distances depending on path: 2 parent choices at each hybrid

# Ancestral state reconstruction with ready-made matrices
params = ParamsBM(10, 1)
srand(2468) # sets the seed for reproducibility, to debug potential error
sim = simulate(net, params)
Y = sim[:Tips]
X = ones(4, 1)
phynetlm = phyloNetworklm(X, Y, net)
@show phynetlm
# Naive version (GLS)
ntaxa = length(Y)
Vy = phynetlm.Vy
Vyinv = inv(Vy)
XtVyinv = X' * Vyinv
logdetVy = logdet(Vy)
betahat = inv(XtVyinv * X) * XtVyinv * Y
fittedValues =  X * betahat
resids = Y - fittedValues
sigma2hat = 1/ntaxa * (resids' * Vyinv * resids)
# log likelihood
loglik = - 1 / 2 * (ntaxa + ntaxa * log(2 * pi) + ntaxa * log(sigma2hat) + logdetVy)
# null version
nullX = ones(ntaxa, 1)
nullXtVyinv = nullX' * Vyinv
nullresids = Y - nullX * inv(nullXtVyinv * nullX) * nullXtVyinv * Y
nullsigma2hat = 1/ntaxa * (nullresids' * Vyinv * nullresids)
nullloglik = - 1 / 2 * (ntaxa + ntaxa * log(2 * pi) + ntaxa * log(nullsigma2hat) + logdetVy)

@test coef(phynetlm) ≈ betahat
@test nobs(phynetlm) ≈ ntaxa
@test residuals(phynetlm) ≈ resids
@test model_response(phynetlm) ≈ Y
@test predict(phynetlm) ≈ fittedValues
@test dof_residual(phynetlm) ≈ ntaxa-length(betahat)
@test sigma2_estim(phynetlm) ≈ sigma2hat
@test loglikelihood(phynetlm) ≈ loglik
@test vcov(phynetlm) ≈ sigma2hat*ntaxa/(ntaxa-length(betahat))*inv(XtVyinv * X)
@test stderr(phynetlm) ≈ sqrt.(diag(sigma2hat*ntaxa/(ntaxa-length(betahat))*inv(XtVyinv * X)))
@test dof(phynetlm) ≈ length(betahat)+1
@test deviance(phynetlm) ≈ sigma2hat * ntaxa
@test nulldeviance(phynetlm) ≈ nullsigma2hat * ntaxa
@test nullloglikelihood(phynetlm) ≈ nullloglik
@test loglikelihood(phynetlm) ≈ nullloglikelihood(phynetlm)
@test deviance(phynetlm) ≈ nulldeviance(phynetlm)
@test r2(phynetlm) ≈ 1-sigma2hat / nullsigma2hat atol=1e-15
@test adjr2(phynetlm) ≈ 1 - (1 - (1-sigma2hat/nullsigma2hat))*(ntaxa-1)/(ntaxa-length(betahat)) atol=1e-15
@test aic(phynetlm) ≈ -2*loglik+2*(length(betahat)+1)
@test aicc(phynetlm) ≈ -2*loglik+2*(length(betahat)+1)+2(length(betahat)+1)*((length(betahat)+1)+1)/(ntaxa-(length(betahat)+1)-1)
@test bic(phynetlm) ≈ -2*loglik+(length(betahat)+1)*log(ntaxa)

# with data frames
dfr = DataFrame(trait = Y, tipNames = sim.M.tipNames)
fitbis = phyloNetworklm(@formula(trait ~ 1), dfr, net)
#@show fitbis

@test coef(phynetlm) ≈ coef(fitbis)
@test vcov(phynetlm) ≈ vcov(fitbis)
@test nobs(phynetlm) ≈ nobs(fitbis)
@test residuals(phynetlm)[fitbis.model.ind] ≈ residuals(fitbis)
@test model_response(phynetlm)[fitbis.model.ind] ≈ model_response(fitbis)
@test predict(phynetlm)[fitbis.model.ind] ≈ predict(fitbis)
@test dof_residual(phynetlm) ≈ dof_residual(fitbis)
@test sigma2_estim(phynetlm) ≈ sigma2_estim(fitbis)
@test stderr(phynetlm) ≈ stderr(fitbis)
@test confint(phynetlm) ≈ confint(fitbis)
@test loglikelihood(phynetlm) ≈ loglikelihood(fitbis)
@test dof(phynetlm) ≈ dof(fitbis)
@test deviance(phynetlm) ≈ deviance(fitbis)
@test nulldeviance(phynetlm) ≈ nulldeviance(fitbis)
@test nullloglikelihood(phynetlm) ≈ nullloglikelihood(fitbis)
@test r2(phynetlm) ≈ r2(fitbis) atol=1e-15
@test adjr2(phynetlm) ≈ adjr2(fitbis) atol=1e-15
@test aic(phynetlm) ≈ aic(fitbis)
@test aicc(phynetlm) ≈ aicc(fitbis)
@test bic(phynetlm) ≈ bic(fitbis)
tmp = (@test_warn "You fitted the data against a custom matrix" mu_estim(phynetlm))
@test tmp ≈ mu_estim(fitbis)

## fixed values parameters
fitlam = phyloNetworklm(@formula(trait ~ 1), dfr, net, model = "lambda", fixedValue=1.0)
@show fitlam

@test lambda_estim(fitlam) ≈ 1.0
@test coef(fitlam) ≈ coef(fitbis)
@test vcov(fitlam) ≈ vcov(fitbis)
@test nobs(fitlam) ≈ nobs(fitbis)
@test residuals(fitlam)[fitbis.model.ind] ≈ residuals(fitbis)
@test model_response(fitlam)[fitbis.model.ind] ≈ model_response(fitbis)
@test predict(fitlam)[fitbis.model.ind] ≈ predict(fitbis)
@test dof_residual(fitlam) ≈ dof_residual(fitbis)
@test sigma2_estim(fitlam) ≈ sigma2_estim(fitbis)
@test stderr(fitlam) ≈ stderr(fitbis)
@test confint(fitlam) ≈ confint(fitbis)
@test loglikelihood(fitlam) ≈ loglikelihood(fitbis)
@test dof(fitlam) ≈ dof(fitbis) + 1
@test deviance(fitlam) ≈ deviance(fitbis)
@test nulldeviance(fitlam) ≈ nulldeviance(fitbis)
@test nullloglikelihood(fitlam) ≈ nullloglikelihood(fitbis)
@test r2(fitlam) ≈ r2(fitbis) atol=1e-15
@test adjr2(fitlam) ≈ adjr2(fitbis) - 0.5 atol=1e-15
@test aic(fitlam) ≈ aic(fitbis) + 2
#@test aicc(fitlam) ≈ aicc(fitbis)
@test bic(fitlam) ≈ bic(fitbis) + log(nobs(fitbis))
@test mu_estim(fitlam) ≈ mu_estim(fitbis)

fitSH = phyloNetworklm(@formula(trait ~ 1), dfr, net, model = "scalingHybrid", fixedValue = 1.0)
@test loglikelihood(fitlam) ≈ loglikelihood(fitSH)
@test aic(fitlam) ≈ aic(fitSH)

## Pagel's Lambda
fitlam = (@test_warn "Maximum lambda value" phyloNetworklm(@formula(trait ~ 1), dfr, net, model = "lambda"))
@test lambda_estim(fitlam) ≈ 1.24875

## Scaling Hybrid
fitSH = phyloNetworklm(@formula(trait ~ 1), dfr, net, model = "scalingHybrid")
@test lambda_estim(fitSH) ≈ 4.057891910001937 atol=1e-8

end

###############################################################################
### With shifts
###############################################################################
@testset "Shifts and Transgressive Evolution" begin

net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");
preorder!(net)

## Simulate
params = ParamsBM(10, 0.1, shiftHybrid([3.0, -3.0],  net))
srand(2468) # sets the seed for reproducibility, to debug potential error
sim = simulate(net, params)
Y = sim[:Tips]

## Construct regression matrix
dfr_shift = regressorShift(net.edge[[8,17]], net)
dfr_shift[:sum] = vec(sum(Array(dfr_shift[:,find(names(dfr_shift) .!= :tipNames)]), 2))
dfr_hybrid = regressorHybrid(net)

@test dfr_shift[:shift_8] ≈ dfr_hybrid[:shift_8]
@test dfr_shift[:shift_17] ≈ dfr_hybrid[:shift_17]
@test dfr_shift[:sum] ≈ dfr_hybrid[:sum]

## Data
dfr = DataFrame(trait = Y, tipNames = sim.M.tipNames)
dfr = join(dfr, dfr_hybrid, on=:tipNames)

## Simple BM
fitShift = phyloNetworklm(@formula(trait ~ shift_8 + shift_17), dfr, net)
@show fitShift

## Test against fixed values lambda models
fitlam = phyloNetworklm(@formula(trait ~ shift_8 + shift_17), dfr, net, model = "lambda", fixedValue = 1.0)

@test lambda_estim(fitlam) ≈ 1.0
@test coef(fitlam) ≈ coef(fitShift)
@test vcov(fitlam) ≈ vcov(fitShift)
@test nobs(fitlam) ≈ nobs(fitShift)
@test residuals(fitlam) ≈ residuals(fitShift)
@test model_response(fitlam) ≈ model_response(fitShift)
@test predict(fitlam) ≈ predict(fitShift)
@test dof_residual(fitlam) ≈ dof_residual(fitShift)
@test sigma2_estim(fitlam) ≈ sigma2_estim(fitShift)
@test stderr(fitlam) ≈ stderr(fitShift)
@test confint(fitlam) ≈ confint(fitShift)
@test loglikelihood(fitlam) ≈ loglikelihood(fitShift)
@test dof(fitlam) ≈ dof(fitShift) + 1
@test deviance(fitlam)  ≈ deviance(fitShift)
@test nulldeviance(fitlam)  ≈ nulldeviance(fitShift)
@test nullloglikelihood(fitlam)  ≈ nullloglikelihood(fitShift)
@test r2(fitlam) ≈ r2(fitShift) atol=1e-15
#@test adjr2(fitlam) ≈ adjr2(fitShift) - 0.5 atol=1e-15
@test aic(fitlam) ≈ aic(fitShift) + 2
#@test aicc(fitlam)  ≈ aicc(fitShift)
@test bic(fitlam) ≈ bic(fitShift) + log(nobs(fitShift))
@test mu_estim(fitlam)  ≈ mu_estim(fitShift)

fitSH = phyloNetworklm(@formula(trait ~ shift_8 + shift_17), dfr, net, model = "scalingHybrid", fixedValue = 1.0)
@test loglikelihood(fitlam) ≈ loglikelihood(fitSH)
@test aic(fitlam) ≈ aic(fitSH)

## ftest against own naive implementation
modnull = phyloNetworklm(@formula(trait ~ 1), dfr, net)
modhom = phyloNetworklm(@formula(trait ~ sum), dfr, net)
modhet = phyloNetworklm(@formula(trait ~ sum + shift_8), dfr, net)

table1 = ftest(modhet, modhom, modnull)
table2 = PhyloNetworks.anova(modnull, modhom, modhet)

@test table1.fstat[1] ≈ table2[:F][2]
@test table1.fstat[2] ≈ table2[:F][1]
@test table1.pval[1].v ≈ table2[Symbol("Pr(>F)")][2]
@test table1.pval[2].v ≈ table2[Symbol("Pr(>F)")][1]
# ## Replace next 4 lines with previous ones when GLM.ftest available
# @test table1[:F][2] ≈ table2[:F][2]
# @test table1[:F][1] ≈ table2[:F][1]
# @test table1[Symbol("Pr(>F)")][1] ≈ table2[Symbol("Pr(>F)")][1]
# @test table1[Symbol("Pr(>F)")][2] ≈ table2[Symbol("Pr(>F)")][2]

# Check that it is the same as doing shift_8 + shift_17
modhetbis = phyloNetworklm(@formula(trait ~ shift_8 + shift_17), dfr, net)

table2bis = PhyloNetworks.anova(modnull, modhom, modhetbis)

@test table2[:F] ≈ table2bis[:F]
@test table2[Symbol("Pr(>F)")] ≈ table2bis[Symbol("Pr(>F)")]
@test table2[:dof_res] ≈ table2bis[:dof_res]
@test table2[:RSS] ≈ table2bis[:RSS]
@test table2[:dof] ≈ table2bis[:dof]
@test table2[:SS] ≈ table2bis[:SS]

end

###############################################################################
#### Other Network
###############################################################################
@testset "phyloNetworklm and ancestralStateReconstruction" begin

# originally: "(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);"
# followed by changes in net.edge[?].length values to make the network ultrametric
net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");
# plot(net, useEdgeLength = true,  showEdgeNumber=true)

#### Simulate correlated data in data frames ####
b0 = 1
b1 = 10
srand(5678)
sim = simulate(net, ParamsBM(1, 1))
A = sim[:Tips]
B = b0 + b1 * A + simulate(net,  ParamsBM(0, 0.1))[:Tips]

# With Matrices
X = hcat(ones(12), A)
fit_mat = phyloNetworklm(X, B, net)
#@show fit_mat

# Naive version (GLS)
ntaxa = length(B)
Vy = fit_mat.Vy
Vyinv = inv(Vy)
XtVyinv = X' * Vyinv
logdetVy = logdet(Vy)
betahat = inv(XtVyinv * X) * XtVyinv * B
fittedValues =  X * betahat
resids = B - fittedValues
sigma2hat = 1/ntaxa * (resids' * Vyinv * resids)
# log likelihood
loglik = - 1 / 2 * (ntaxa + ntaxa * log(2 * pi) + ntaxa * log(sigma2hat) + logdetVy)
# null version
nullX = ones(ntaxa, 1)
nullXtVyinv = nullX' * Vyinv
nullresids = B - nullX * inv(nullXtVyinv * nullX) * nullXtVyinv * B
nullsigma2hat = 1/ntaxa * (nullresids' * Vyinv * nullresids)
nullloglik = - 1 / 2 * (ntaxa + ntaxa * log(2 * pi) + ntaxa * log(nullsigma2hat) + logdetVy)
@test coef(fit_mat) ≈ betahat
@test nobs(fit_mat) ≈ ntaxa
@test residuals(fit_mat) ≈ resids
@test model_response(fit_mat) ≈ B
@test predict(fit_mat) ≈ fittedValues
@test dof_residual(fit_mat) ≈ ntaxa-length(betahat)
@test sigma2_estim(fit_mat) ≈ sigma2hat
@test loglikelihood(fit_mat) ≈ loglik
@test vcov(fit_mat) ≈ sigma2hat*ntaxa/(ntaxa-length(betahat)).*inv(XtVyinv * X)
@test stderr(fit_mat) ≈ sqrt.(diag(sigma2hat*ntaxa/(ntaxa-length(betahat)).*inv(XtVyinv * X)))
@test dof(fit_mat) ≈ length(betahat)+1
@test deviance(fit_mat) ≈ sigma2hat * ntaxa
@test nulldeviance(fit_mat) ≈ nullsigma2hat * ntaxa
@test nullloglikelihood(fit_mat) ≈ nullloglik
@test r2(fit_mat) ≈ 1-sigma2hat / nullsigma2hat atol=1e-15
@test adjr2(fit_mat) ≈ 1 - (1 - (1-sigma2hat/nullsigma2hat))*(ntaxa-1)/(ntaxa-length(betahat)) atol=1e-15
@test aic(fit_mat) ≈ -2*loglik+2*(length(betahat)+1)
@test aicc(fit_mat) ≈ -2*loglik+2*(length(betahat)+1)+2(length(betahat)+1)*((length(betahat)+1)+1)/(ntaxa-(length(betahat)+1)-1)
@test bic(fit_mat) ≈ -2*loglik+(length(betahat)+1)*log(ntaxa)

## perfect user using right format and formula
dfr = DataFrame(trait = B, pred = A, tipNames = sim.M.tipNames)
phynetlm = phyloNetworklm(@formula(trait ~ pred), dfr, net)
#@show phynetlm

@test coef(phynetlm) ≈ coef(fit_mat)
@test vcov(phynetlm) ≈ vcov(fit_mat)
@test nobs(phynetlm) ≈ nobs(fit_mat)
@test residuals(phynetlm) ≈ residuals(fit_mat)
@test model_response(phynetlm) ≈ model_response(fit_mat)
@test predict(phynetlm) ≈ predict(fit_mat)
@test dof_residual(phynetlm) ≈ dof_residual(fit_mat)
@test sigma2_estim(phynetlm) ≈ sigma2_estim(fit_mat)
@test stderr(phynetlm) ≈ stderr(fit_mat)
@test confint(phynetlm) ≈ confint(fit_mat)
@test loglikelihood(phynetlm) ≈ loglikelihood(fit_mat)
@test dof(phynetlm) ≈ dof(fit_mat)
@test deviance(phynetlm) ≈ deviance(fit_mat)
@test nulldeviance(phynetlm) ≈ nulldeviance(fit_mat)
@test nullloglikelihood(phynetlm) ≈ nullloglikelihood(fit_mat)
@test r2(phynetlm) ≈ r2(fit_mat)
@test adjr2(phynetlm) ≈ adjr2(fit_mat)
@test aic(phynetlm) ≈ aic(fit_mat)
@test aicc(phynetlm) ≈ aicc(fit_mat)
@test bic(phynetlm) ≈ bic(fit_mat)

# unordered data
srand(1234)
dfr = dfr[sample(1:12, 12, replace=false), :]
fitbis = phyloNetworklm(@formula(trait ~ pred), dfr, net)

@test coef(phynetlm) ≈ coef(fitbis)
@test vcov(phynetlm) ≈ vcov(fitbis)
@test nobs(phynetlm) ≈ nobs(fitbis)
@test residuals(phynetlm)[fitbis.model.ind] ≈ residuals(fitbis)
@test model_response(phynetlm)[fitbis.model.ind] ≈ model_response(fitbis)
@test predict(phynetlm)[fitbis.model.ind] ≈ predict(fitbis)
@test dof_residual(phynetlm) ≈ dof_residual(fitbis)
@test sigma2_estim(phynetlm) ≈ sigma2_estim(fitbis)
@test stderr(phynetlm) ≈ stderr(fitbis)
@test confint(phynetlm) ≈ confint(fitbis)
@test loglikelihood(phynetlm) ≈ loglikelihood(fitbis)
@test dof(phynetlm) ≈ dof(fitbis)
@test deviance(phynetlm) ≈ deviance(fitbis)
@test nulldeviance(phynetlm) ≈ nulldeviance(fitbis)
@test nullloglikelihood(phynetlm) ≈ nullloglikelihood(fitbis)
@test r2(phynetlm) ≈ r2(fitbis)
@test adjr2(phynetlm) ≈ adjr2(fitbis)
@test aic(phynetlm) ≈ aic(fitbis)
@test aicc(phynetlm) ≈ aicc(fitbis)
@test bic(phynetlm) ≈ bic(fitbis)
@test mu_estim(phynetlm) ≈ mu_estim(fitbis)

# unnamed ordered data
dfr = DataFrame(trait = B, pred = A)
fitter = (@test_warn "As requested (no_names=true)" phyloNetworklm(@formula(trait ~ pred), dfr, net, no_names=true))

@test coef(phynetlm) ≈ coef(fitter)
@test vcov(phynetlm) ≈ vcov(fitter)
@test nobs(phynetlm) ≈ nobs(fitter)
@test residuals(phynetlm) ≈ residuals(fitter)
@test model_response(phynetlm) ≈ model_response(fitter)
@test predict(phynetlm) ≈ predict(fitter)
@test dof_residual(phynetlm) ≈ dof_residual(fitter)
@test sigma2_estim(phynetlm) ≈ sigma2_estim(fitter)
@test stderr(phynetlm) ≈ stderr(fitter)
@test confint(phynetlm) ≈ confint(fitter)
@test loglikelihood(phynetlm) ≈ loglikelihood(fitter)
@test dof(phynetlm) ≈ dof(fitter)
@test deviance(phynetlm) ≈ deviance(fitter)
@test nulldeviance(phynetlm) ≈ nulldeviance(fitter)
@test nullloglikelihood(phynetlm) ≈ nullloglikelihood(fitter)
@test r2(phynetlm) ≈ r2(fitter)
@test adjr2(phynetlm) ≈ adjr2(fitter)
@test aic(phynetlm) ≈ aic(fitter)
@test aicc(phynetlm) ≈ aicc(fitter)
@test bic(phynetlm) ≈ bic(fitter)

# unnamed un-ordered data
dfr = dfr[sample(1:12, 12, replace=false), :]
@test_throws ErrorException phyloNetworklm(@formula(trait ~ pred), dfr, net) # Wrong pred

### Add NAs
dfr = DataFrame(trait = B, pred = A, tipNames = tipLabels(sim))
dfr[:pred] = allowmissing(dfr[:pred])
dfr[[2, 8, 11], :pred] = missing
fitna = phyloNetworklm(@formula(trait ~ pred), dfr, net)
#@show fitna

dfr = dfr[sample(1:12, 12, replace=false), :]
fitnabis = phyloNetworklm(@formula(trait ~ pred), dfr, net)

@test coef(fitna) ≈ coef(fitnabis)
@test vcov(fitna) ≈ vcov(fitnabis)
@test nobs(fitna) ≈ nobs(fitnabis)
@test sort(residuals(fitna)) ≈ sort(residuals(fitnabis))
@test sort(model_response(fitna)) ≈ sort(model_response(fitnabis))
@test sort(predict(fitna)) ≈ sort(predict(fitnabis))
@test dof_residual(fitna) ≈ dof_residual(fitnabis)
@test sigma2_estim(fitna) ≈ sigma2_estim(fitnabis)
@test stderr(fitna) ≈ stderr(fitnabis)
@test confint(fitna) ≈ confint(fitnabis)
@test loglikelihood(fitna) ≈ loglikelihood(fitnabis)
@test dof(fitna) ≈ dof(fitnabis)
@test deviance(fitna) ≈ deviance(fitnabis)
@test nulldeviance(fitna) ≈ nulldeviance(fitnabis)
@test nullloglikelihood(fitna) ≈ nullloglikelihood(fitnabis)
@test r2(fitna) ≈ r2(fitnabis)
@test adjr2(fitna) ≈ adjr2(fitnabis)
@test aic(fitna) ≈ aic(fitnabis)
@test aicc(fitna) ≈ aicc(fitnabis)
@test bic(fitna) ≈ bic(fitnabis)

## Tests against fixed values parameters
fitlam = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "lambda", fixedValue = 1.0)
#@show fitlam

@test lambda_estim(fitlam) ≈ 1.0
@test coef(fitlam) ≈ coef(fitnabis)
@test vcov(fitlam) ≈ vcov(fitnabis)
@test nobs(fitlam) ≈ nobs(fitnabis)
@test residuals(fitlam) ≈ residuals(fitnabis)
@test model_response(fitlam) ≈ model_response(fitnabis)
@test predict(fitlam) ≈ predict(fitnabis)
@test dof_residual(fitlam) ≈ dof_residual(fitnabis)
@test sigma2_estim(fitlam) ≈ sigma2_estim(fitnabis)
@test stderr(fitlam) ≈ stderr(fitnabis)
@test confint(fitlam) ≈ confint(fitnabis)
@test loglikelihood(fitlam) ≈ loglikelihood(fitnabis)
@test dof(fitlam) ≈ dof(fitnabis) + 1
@test deviance(fitlam) ≈ deviance(fitnabis)
@test nulldeviance(fitlam) ≈ nulldeviance(fitnabis)
@test nullloglikelihood(fitlam) ≈ nullloglikelihood(fitnabis)
@test r2(fitlam) ≈ r2(fitnabis) atol=1e-15
@test adjr2(fitlam)-1 ≈ (adjr2(fitnabis)-1)*(nobs(fitnabis)-dof(fitnabis)+1)/(nobs(fitnabis)-dof(fitlam)+1) atol=1e-15
@test aic(fitlam) ≈ aic(fitnabis) + 2
#@test aicc(fitlam) ≈ aicc(fitnabis)
@test bic(fitlam) ≈ bic(fitnabis) + log(nobs(fitnabis))
@test mu_estim(fitlam) ≈ mu_estim(fitnabis)

fitSH = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "scalingHybrid", fixedValue = 1.0)
@test loglikelihood(fitlam) ≈ loglikelihood(fitSH)
@test aic(fitlam) ≈ aic(fitSH)

## Pagel's Lambda
fitlam = (@test_warn "Maximum lambda value" phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "lambda"))
#@show fitlam
@test lambda_estim(fitlam) ≈ 1.1135518305 atol=1e-10

## scaling Hybrid
fitSH = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "scalingHybrid")
@show fitSH
@test lambda_estim(fitSH) ≈ -52.81305448333567 atol=1e-8

### Ancestral State Reconstruction
params = ParamsBM(3, 1)
sim = simulate(net, params)
Y = sim[:Tips]
# From known parameters
ancestral_traits = ancestralStateReconstruction(net, Y, params)
# BLUP
dfr = DataFrame(trait = Y, tipNames = tipLabels(sim))
phynetlm = phyloNetworklm(@formula(trait~1), dfr, net)
blup = (@test_warn "These prediction intervals show uncertainty in ancestral values" ancestralStateReconstruction(phynetlm));
# plot(net, blup)
@show blup

# BLUP same, using the function dirrectly
blup_bis = (@test_warn "These prediction intervals show uncertainty in ancestral values" ancestralStateReconstruction(dfr, net));

@test expectations(blup)[:condExpectation] ≈ expectations(blup_bis)[:condExpectation]
@test expectations(blup)[:nodeNumber] ≈ expectations(blup_bis)[:nodeNumber]
@test blup.traits_tips ≈ blup_bis.traits_tips
@test blup.TipNumbers ≈ blup_bis.TipNumbers
@test predint(blup) ≈ predint(blup_bis)
@test predintPlot(blup)[:PredInt] == predintPlot(blup_bis)[:PredInt]
@test predintPlot(blup, withExp=true)[:PredInt] == predintPlot(blup_bis, withExp=true)[:PredInt]
@test expectationsPlot(blup)[:PredInt] == expectationsPlot(blup_bis)[:PredInt]

dfr = DataFrame(trait = Y, tipNames = tipLabels(sim), reg = Y)
@test_throws ErrorException ancestralStateReconstruction(dfr, net) # cannot handle a predictor

# Unordered
dfr2 = dfr[sample(1:12, 12, replace=false), :]
phynetlm = phyloNetworklm(@formula(trait~1), dfr2, net)
blup2 = (@test_warn "These prediction intervals show uncertainty in ancestral values" ancestralStateReconstruction(phynetlm))

@test expectations(blup)[:condExpectation][1:length(blup.NodeNumbers)] ≈ expectations(blup2)[:condExpectation][1:length(blup.NodeNumbers)]
@test blup.traits_tips[phynetlm.model.ind] ≈ blup2.traits_tips
@test blup.TipNumbers[phynetlm.model.ind] ≈ blup2.TipNumbers
@test predint(blup)[1:length(blup.NodeNumbers), :] ≈ predint(blup2)[1:length(blup.NodeNumbers), :]

# With unknown tips
dfr[:trait] = allowmissing(dfr[:trait])
dfr[[2, 4], :trait] = missing
phynetlm = phyloNetworklm(@formula(trait~1), dfr, net)
blup = (@test_warn "These prediction intervals show uncertainty in ancestral values" ancestralStateReconstruction(phynetlm))
# plot(net, blup)

# Unordered
dfr2 = dfr[[1, 2, 5, 3, 4, 6, 7, 8, 9, 10, 11, 12], :]
phynetlm = phyloNetworklm(@formula(trait~1), dfr, net)
blup2 = (@test_warn "These prediction intervals show uncertainty in ancestral values" ancestralStateReconstruction(phynetlm))

@test expectations(blup)[:condExpectation] ≈ expectations(blup2)[:condExpectation]
@test predint(blup) ≈ predint(blup2)
@test predintPlot(blup)[:PredInt] == predintPlot(blup2)[:PredInt]
@test predintPlot(blup, withExp=true)[:PredInt] == predintPlot(blup2, withExp=true)[:PredInt]

# Test mark on missing
ee = expectationsPlot(blup)
predMiss = ee[indexin([n.number for n in net.leaf][[2,4]], ee[:nodeNumber]),:PredInt]
for pp = predMiss
    @test pp[end] == '*'
end

end

#################
## Data with no phylogenetic signal
#################
@testset "lambda when no signal" begin

net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");
# plot(net, useEdgeLength = true,  showEdgeNumber=true)

#### Simulate correlated data in data frames ####
b0 = 1
b1 = 10
srand(5678)
A = randn(size(tipLabels(net), 1))
B = b0 + b1 * A + randn(size(tipLabels(net), 1))
dfr = DataFrame(trait = B, pred = A, tipNames = tipLabels(net))

## Network
phynetlm = (@test_warn "Maximum lambda value" phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "lambda"))

@test lambda_estim(phynetlm) ≈ 0.5894200143 atol=1e-8

## Major Tree
tree = majorTree(net)
phynetlm = (@test_warn "Maximum lambda value" phyloNetworklm(@formula(trait ~ pred), dfr, tree, model = "lambda"))

@test lambda_estim(phynetlm) ≈ 0.5903394415 atol=1e-6

## scaling Hybrid
lmtree = phyloNetworklm(@formula(trait ~ pred), dfr, tree, model = "BM")
lmnet = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "BM")
lmSHzero = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "scalingHybrid", fixedValue = 0.0)
lmSHone = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "scalingHybrid", fixedValue = 1.0)

@test loglikelihood(lmtree) ≈ loglikelihood(lmSHzero)
@test loglikelihood(lmnet) ≈ loglikelihood(lmSHone)

lmSH = phyloNetworklm(@formula(trait ~ pred), dfr, net, model = "scalingHybrid")

@test lambda_estim(lmSH) ≈ 23.46668204551696 atol=1e-5

end

############################
## Against no regressor
###########################
@testset "phloNetworklm with no regressor" begin

net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");

params = ParamsBM(10, 1)
srand(2468) # sets the seed for reproducibility, to debug potential error
sim = simulate(net, params)
Y = sim[:Tips]
phynetlm = phyloNetworklm(zeros(length(Y),0), Y, net)
#@show phynetlm
# Naive version (GLS)
ntaxa = length(Y)
Vy = phynetlm.Vy
Vyinv = inv(Vy)
logdetVy = logdet(Vy)
fittedValues =  zeros(length(Y))
resids = Y - fittedValues
sigma2hat = 1/ntaxa * (resids' * Vyinv * resids)
# log likelihood
loglik = - 1 / 2 * (ntaxa + ntaxa * log(2 * pi) + ntaxa * log(sigma2hat) + logdetVy)
# null version
nullX = ones(ntaxa, 1)
nullXtVyinv = nullX' * Vyinv
nullresids = Y - nullX * inv(nullXtVyinv * nullX) * nullXtVyinv * Y
nullsigma2hat = 1/ntaxa * (nullresids' * Vyinv * nullresids)
nullloglik = - 1 / 2 * (ntaxa + ntaxa * log(2 * pi) + ntaxa * log(nullsigma2hat) + logdetVy)

@test nobs(phynetlm) ≈ ntaxa
@test residuals(phynetlm) ≈ resids
@test model_response(phynetlm) ≈ Y
@test predict(phynetlm) ≈ fittedValues
@test dof_residual(phynetlm) ≈ ntaxa
@test sigma2_estim(phynetlm) ≈ sigma2hat
@test loglikelihood(phynetlm) ≈ loglik
@test deviance(phynetlm) ≈ sigma2hat * ntaxa
@test nulldeviance(phynetlm) ≈ nullsigma2hat * ntaxa
@test nullloglikelihood(phynetlm) ≈ nullloglik
@test r2(phynetlm) ≈ 1-sigma2hat / nullsigma2hat atol=1e-14
@test adjr2(phynetlm) ≈ 1 - (1 - (1-sigma2hat/nullsigma2hat))*(ntaxa-1)/(ntaxa) atol=1e-14
@test aic(phynetlm) ≈ -2*loglik+2*(1)
@test aicc(phynetlm) ≈ -2*loglik+2*(1)+2(1)*((1)+1)/(ntaxa-(1)-1)
@test bic(phynetlm) ≈ -2*loglik+(1)*log(ntaxa)

# with data frames
dfr = DataFrame(trait = Y, tipNames = sim.M.tipNames)
fitbis = phyloNetworklm(@formula(trait ~ -1), dfr, net)
@show fitbis
#@test coef(phynetlm) ≈ coef(fitbis)
#@test vcov(phynetlm) ≈ vcov(fitbis)
@test nobs(phynetlm) ≈ nobs(fitbis)
@test residuals(phynetlm)[fitbis.model.ind] ≈ residuals(fitbis)
@test model_response(phynetlm)[fitbis.model.ind] ≈ model_response(fitbis)
@test predict(phynetlm)[fitbis.model.ind] ≈ predict(fitbis)
@test dof_residual(phynetlm) ≈ dof_residual(fitbis)
@test sigma2_estim(phynetlm) ≈ sigma2_estim(fitbis)
#@test stderr(phynetlm) ≈ stderr(fitbis)
#@test confint(phynetlm) ≈ confint(fitbis)
@test loglikelihood(phynetlm) ≈ loglikelihood(fitbis)
#@test dof(phynetlm) ≈ dof(fitbis)
@test deviance(phynetlm) ≈ deviance(fitbis)
@test nulldeviance(phynetlm) ≈ nulldeviance(fitbis)
@test nullloglikelihood(phynetlm) ≈ nullloglikelihood(fitbis)
@test r2(phynetlm) ≈ r2(fitbis) atol=1e-15
@test adjr2(phynetlm) ≈ adjr2(fitbis) atol=1e-15
@test aic(phynetlm) ≈ aic(fitbis)
@test aicc(phynetlm) ≈ aicc(fitbis)
@test bic(phynetlm) ≈ bic(fitbis)
#@test mu_estim(phynetlm)  mu_estim(fitbis)

end


text:
using LinearAlgebra, SparseArrays, ForwardDiff
using Convex, ECOS

# Dimensions
n = 2 # state
m = 1 # control
T = 2 # horizon
num_var = n * T + m * (T - 1)
num_con = n * T
num_dec = num_var + num_con
z0 = rand(num_dec)

function traj(z)
    x1 = view(z, 1:n)
    u1 = view(z, n .+ (1:m))
    x2 = view(z, n + m .+ (1:n))
    λ1 = view(z, 2 * n + m .+ (1:n))
    λ2 = view(z, 3 * n + m .+ (1:n))
    return x1, u1, x2, λ1, λ2
end

x1, u1, x2, λ1, λ2 = traj(z0)

# Dynamics
x_init = ones(n)

A = [1.0 1.0;
     0.0 1.0]

B = [0.1;
     0.0]

In = Diagonal(ones(n))

# Objective
Q = [1.0 0.0; 0.0 1.0]
R = Array([1.0e-1])

# Convex.jl
y = Variable(num_var) # decision variables
objective = 0.5 * quadform(y[1:n], Q) + 0.5 * quadform(y[n + m .+ (1:n)], Q) + 0.5 * R[1] * square(y[n .+ (1:m)][1])

problem = minimize(objective) # setup problem

# constraint setup
problem.constraints += A * y[1:n] + B * y[n .+ (1:m)][1] - y[n + m .+ (1:n)] == 0.0
problem.constraints += y[1:n] - x_init == 0.0

solve!(problem, ECOS.Optimizer) # solve

# Convex solution
x1_sol = y.value[1:n]
u1_sol = y.value[n .+ (1:m)]
x2_sol = y.value[n + m .+ (1:n)]
λ1_sol = -1.0 * problem.constraints[1].dual
λ2_sol = -1.0 * problem.constraints[2].dual

z_sol = vec([x1_sol; u1_sol; x2_sol; λ1_sol; λ2_sol])

# kkt system r + Hz = 0
function r(z)
    x1, u1, x2, λ1, λ2 = traj(z)
    [Q * x1 + A' * λ1 + λ2;
     (R * u1[1])[1] + B' * λ1;
     Q * x2 - λ1;
     A * x1 + B * u1[1] - x2;
     x1 - x_init]
end
r(z_sol)

H = [Q zeros(n, m) zeros(n, n) A' In;
     zeros(m, n) R zeros(m, n) B' zeros(m, n);
     zeros(n, n) zeros(n, m) Q -In zeros(n, n);
     A B -In zeros(n, n) zeros(n, n);
     In zeros(n, m) zeros(n, n) zeros(n, n) zeros(n, n)]

r(z0)
ForwardDiff.jacobian(r, z0) - H

# function con(z)
#     x1, u1, x2, λ1, λ2 = traj(z)
#
#     [A * x1 + B * u1[1] - x2
#      x1 - ones(2);]
# end
#
# C = [A B -In;
#      In zeros(n, m + n)]
#
# con(z0)
# ForwardDiff.jacobian(con, z0)[:, 1:num_var] - C

# initialize variables
x1 = rand(n) #copy(x_init) #
u1 = rand(m)
x2 = A * x1 + B * u1[1]
λ1 = rand(n) # copy(λ1_sol)
λ2 = rand(n) # copy(λ2_sol)

z = [x1; u1; x2; λ1; λ2]

# kkt solve
Δz = -1.0 * H \ r(z)
Δx1, Δu1, Δx2, Δλ1, Δλ2 = traj(Δz)

# iLQR solve
function compute_Δu1(Δx1)
    Quu = R[1] + B' * Q * B
    Qux = B' * Q * A
    Qu = (R * u1[1])[1] + B' * Q * x2
    -1.0 * Quu \ (Qux * Δx1 + Qu)
end

Δu1[1] - compute_Δu1(Δx1)


text:
module TestTypeChangingAccumulators

using Adapt
using CUDA
using Folds
using Folds.Testing: Testing, parse_tests
using FoldsCUDA
using Test
using Transducers

missing_if_odd(x) = isodd(x) ? missing : x

function test_sum()
    @test Folds.sum(missing_if_odd, 0:2:2^10, CUDAEx()) == sum(0:2:2^10)
    @test Folds.sum(missing_if_odd, 0:2:2^20, CUDAEx()) == sum(0:2:2^20)
    @test Folds.sum(missing_if_odd, 0:2^10, CUDAEx()) === missing
    @test Folds.sum(missing_if_odd, 0:2^20, CUDAEx()) === missing
end

partition_length_maximum(xs, ex = PreferParallel()) = Transducers.fold(
    max,
    xs |> ReducePartitionBy(identity, Map(_ -> 1)'(+), 0),
    ex;
    init = typemin(Int),
)

function test_partition_length_maximum()
    @testset "2^$e" for e in [5, 10, 15, 20, 25]
        xs = CUDA.rand(Bool, 2^e)
        @test partition_length_maximum(xs) == partition_length_maximum(collect(xs))
    end
end

rawdata = """
reduce(TeeRF(min, max), (2x for x in 1:10 if isodd(x)))
"""

#=
prod(y for x in 1:11 if isodd(x) for y in 1:x:x^2)
=#

testdata = parse_tests(rawdata, @__MODULE__)

upload(x) = adapt(CuArray, x)

test_with_sequential() = Testing.test_with_sequential(testdata, [CUDAEx()]; upload = upload)

end  # module


text:
module TrafficAssign

# package code goes here
using LightGraphs, Optim, BinDeps
using Distributed, Printf, LinearAlgebra, SparseArrays

TNTP_SHA = "ca933f74c03af688dfffa33f4b00826cba672479"

include("load_network.jl")
include("frank_wolfe.jl")

export
        load_ta_network, download_tntp, read_ta_network,
        ta_frank_wolfe,
        TA_Data, solve_TAP

end # module

max_stars_repo_name: salomonw/tnet

max_issues_repo_name: salomonw/tnet

max_forks_repo_name: salomonw/tnet


text:
import Pkg

Pkg.add(path="https://github.com/jmurga/Analytical.jl")
Pkg.add("CSV")
Pkg.add("JLD2")
Pkg.add("DataFrames")
Pkg.add("GZip")
Pkg.add("RCall")
Pkg.add("ProgressMeter")
Pkg.add("Fire")
Pkg.add("ParallelUtilities")