text:
"""
    couple_models(aero::LiftingLine, stru::GEBT, dyn::RigidBody, flap::LiftingLineFlaps)

Create a coupled model using a lifting line aerodynamic model, a
geometrically exact beam theory model, a rigid body dynamics model, and a lifting
line control surface model.  This model introduces additional parameters
corresponding to the freestream air density ``\\rho_\\infty``, followed by the
external loads ``F_{x,i}, F_{y,i}, F_{z,i}, M_{x,i}, M_{y,i}, M_{z,i}`` or
displacements ``u_{x,i}, u_{y,i}, u_{z,i}, \\theta_{x,i}, \\theta_{y,i},
\\theta_{z,i}`` for each node, followed by the constant distributed loads
``f_{x,i}, f_{y,i}, f_{z,i}, m_{x,i}, m_{y,i}, m_{z,i}`` applied on each beam
element (excluding aerodynamic loads), followed by the control surface
deflections ``\\delta_1, \\delta_2, \\dots, \\delta_N``.

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

    return (aero, stru, dyn, flap)
end

# --- traits --- #

function number_of_additional_parameters(aero::LiftingLine, stru::GEBT, dyn::RigidBody,
    flap::LiftingLineFlaps)

    return 1 + 6*length(stru.icol_point) + 6*length(stru.icol_elem) + length(flap.gains)
end

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

    return InPlace()
end

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

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

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

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

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

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

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

# --- methods --- #

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

    # extract number of state variables, inputs, and parameters
    nxa = number_of_states(aero) # number of aerodynamic states
    nxs = number_of_states(stru) # number of structural states
    nxd = number_of_states(dyn) # number of rigid body states
    nxf = number_of_states(flap) # number of control surface states
    nya = number_of_inputs(aero) # number of aerodynamic inputs
    nys = number_of_inputs(stru) # number of structural inputs
    nyd = number_of_inputs(dyn) # number of rigid body inputs
    nyf = number_of_inputs(flap) # number of control surface inputs
    npa = number_of_parameters(aero) # number of aerodynamic parameters
    nps = number_of_parameters(stru) # number of structural parameters
    npd = number_of_parameters(dyn) # number of rigid body parameters
    npf = number_of_parameters(flap) # number of control surface parameters
    npadd = number_of_additional_parameters(aero, stru, dyn, flap) # number of additional parameters

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

    # separate state variables, inputs, and parameters
    dxa = view(dx, ixa) # aerodynamic rates
    xa = view(x, ixa) # aerodynamic states
    ya = view(y, iya) # aerodynamic inputs
    pa = view(p, ipa) # aerodynamic parameters
    dxs = view(dx, ixs) # structural rates
    xs = view(x, ixs) # structural states
    ys = view(y, iys) # structural inputs
    ps = view(p, ips) # structural parameters
    dxd = view(dx, ixd) # rigid body rates
    xd = view(x, ixd) # rigid body states
    yd = view(y, iyd) # rigid body inputs
    pd = view(p, ipd) # rigid body parameters
    dxf = view(dx, ixf) # structural rates
    xf = view(x, ixf) # structural states
    yf = view(y, iyf) # structural inputs
    pf = view(p, ipf) # structural parameters
    padd = view(p, ipadd) # additional parameters for coupled model
    dxas = view.(Ref(dxa), ixas) # aerodynamic rates for each section
    xas = view.(Ref(xa), ixas) # aerodynamic states for each section
    yas = view.(Ref(ya), iyas) # aerodynamic inputs for each section
    pas = view.(Ref(pa), ipas) # aerodynamic parameters for each section
    dxfs = view.(Ref(dxf), ixfs) # control surface rates for each section
    xfs = view.(Ref(xf), ixfs) # control surface states for each section
    yfs = view.(Ref(yf), iyfs) # control surface inputs for each section
    pfs = view.(Ref(pf), ipfs) # control surface parameters for each section

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

    # rigid body states
    xr, yr, zr, ϕr, θr, ψr, ur, vr, wr, pr, qr, rr = xd

    # rigid body state rates
    dxr, dyr, dzr, dϕr, dθr, dψr, dur, dvr, dwr, dpr, dqr, drr = dxd

    # extract global parameters
    ρ = padd[1] # freestream air density
    δ = padd[SVector{NG}(1 + 6*npoint + 6*nelem + 1 : 1 + 6*npoint + 6*nelem + NG)]
    dδ = zero(δ)

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

    # rigid body linear and angular acceleration
    dV = SVector(dur, dvr, dwr)
    dΩ = SVector(dpr, dqr, drr)

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

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

    # initialize rigid body properties
    mass = 0.0
    Ir = @SMatrix zeros(3,3)

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

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

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

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

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

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

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

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

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

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

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

        # local freestream linear and angular accelerations
        dvi = -R*(CtCab'*dV + CtCab'*cross(dΩ, pe) + GXBeam.element_linear_velocity(element, dP_elem, dH_elem))
        dωi = R*CtCab'*dΩ + R*GXBeam.element_angular_velocity(element, dP_elem, dH_elem)

        # element inertial properties
        poff = 3*npoint + 36*(i-1)
        μ = ps[poff + 31] # mass per unit length
        xm2 = ps[poff + 32] # center of mass location
        xm3 = ps[poff + 33] # center of mass location
        i22 = ps[poff + 34] # rotational inertia
        i33 = ps[poff + 35] # rotational inertia
        i23 = ps[poff + 36] # product of inertia
        me = ΔL*μ # element mass
        Ie = ΔL*(@SMatrix [i22+i33 0 0; 0 i22 -i23; 0 -i23 i33]) # local frame inertia
        Ie = (CtCab*R')*Ie*(R*CtCab') # body frame inertia
        Ie = me*(pe'*pe*I - pe*pe') + Ie # body frame inertia about origin

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

        # section state variable rates
        duai = SVector{Nxai}(dxas[i]) # aerodynamic state variables
        dusi = vcat(dvi, dωi) # structural state variables
        dufi = SVector{Nxfi}(dxfs[i]) # control surface state variables
        duci = SMatrix{1,NG}(getindex.(flap.gains, i)) * dδ # controller state variables
        dui = vcat(duai, dusi, dufi, duci)

        # section state variables
        uai = SVector{Nxai}(xas[i]) # aerodynamic state variables
        usi = vcat(vi, ωi) # structural state variables
        ufi = SVector{Nxfi}(xfs[i]) # control surface state variables
        uci = SMatrix{1,NG}(getindex.(flap.gains, i)) * δ # controller state variables
        ui = vcat(uai, usi, ufi, uci)

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

        # section inputs
        yi = get_coupling_inputs(section_models, dui, ui, pi, t)

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

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

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

        # add apparent forces due to body frame linear and angular acceleration
        fi += me*dV - cross(me*dΩ, pe)
        mi += me*dΩ

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

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

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

        # add element mass and inertia to total mass and inertia
        mass += me
        Ir += Ie
    end

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

    # save rigid body model inputs
    y[iyd[1]] = mass
    y[iyd[2]] = Ixx = Ir[1,1]
    y[iyd[3]] = Iyy = Ir[2,2]
    y[iyd[4]] = Izz = Ir[3,3]
    y[iyd[5]] = Ixz = Ir[1,3]
    y[iyd[6]] = Ixy = Ir[1,2]
    y[iyd[7]] = Iyz = Ir[2,3]
    y[iyd[8:10]] = Ftot
    y[iyd[11:13]] = Mtot

    return y
end

# --- performance overloads --- #

# TODO

# --- convenience methods --- #

function set_additional_parameters!(padd, model::LiftingLine, stru::GEBT,
    dyn::RigidBody, flap::LiftingLineFlaps{NF,NG,TF,TG}; rho, point_conditions,
    element_loads, delta) where {NF,NG,TF,TG}

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

    padd[1] = rho

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

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

    padd[1 + 6*np + 6*ne + 1 : 1 + 6*np + 6*ne + NG] = delta

    return padd
end

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

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

    rho = padd[1]

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

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

    delta = view(padd, 1 + 6*np + 6*ne + 1 : 1 + 6*np + 6*ne + NG)

    return (rho = rho, point_conditions = point_conditions, element_loads = element_loads,
        delta = delta)
end


text:
function printDate(year, month, day)

    print("Today's date is $year-$month-$day\n")

end

myFun(first, second) = first^second

ext: jl


text:
module JDFManager

using DataFrameTools
using SHA
using Dates

export cached_df_fn, from_hash, to_hash, to_metaname, from_metaname, set_jdf_dir

jdf_dir = tempdir()

function set_jdf_dir(d)
	global jdf_dir
	jdf_dir = d
end

cache_dir(dir, hash) = joinpath(jdf_dir, dir, hash * ".jdf")

function cached_df_fn(metaname, fn; dir="JDF", metadata="", useCache=true)
	if useCache
		df = from_metaname(metaname, dir=dir)
		if df != nothing
			return df
		end
	end
	df = fn()
	if df == nothing
		return df
	end
	to_metaname(df, metaname, dir=dir, metadata=metadata)
	df
end

from_metaname(metaname; dir="JDF") = from_hash(bytes2hex(sha256(metaname)), dir=dir, metaname=metaname)

function from_hash(hash; dir="JDF", metaname="")
	jdf = cache_dir(dir, hash)
	if isdir(jdf)
		println(stderr, "Cached ", metaname, " - ", joinpath(dir, hash * ".jdf"))
		return DataFrameTools.df_read(jdf)
	end
end

to_metaname(df, metaname; dir="JDF", metadata="") = to_hash(df, bytes2hex(sha256(metaname)), dir=dir, metafilename=metaname)

function to_hash(df, hash; dir="JDF", metafilename="", metadata="")
	if df == nothing
		return
	end
	jdf = cache_dir(dir, hash)
	println(stderr, "Caching ", metafilename, " - ", joinpath(dir, hash * ".jdf"))
	println(stderr, "DateTime -> Date")
	DataFrameTools.df_write(jdf, roundtime(df))
	if metafilename != ""
		open(joinpath(jdf, replace(metafilename, ":"=>".") ), "w+") do fid
			write(fid, metadata)
		end
	end
end

function roundtime(df)
	for column in names(df)
		if typeof(df[!, column]) == Vector{Union{Missing, Dates.DateTime}}
			df[!, column] = map(dt-> ismissing(dt) ? missing : convert(Date, dt), df[!, column])
		end
	end
	df
end

######
end


text:
# function IntersectionOfBoundaries(s₁::Array{Float64, 2}, s₂::Array{Float64, 2},
#                                     convexexp1in2::Array{Float64, 2},
#                                     convexexp2in1::Array{Float64, 2},
#                                     ordered_vertices1::RowVector{Int},
#                                     ordered_vertices2::RowVector{Int},
#                                     num1in2::Int,
#                                     num2in1::Int,
#                                     ncomm::Int,
#                                     tol::Float64)
#
#
# end

function IntersectionOfBoundaries_NoStorage(s₁,
                                    s₂,
                                    convexexp1in2::AbstractArray{Float64, 2},
                                    convexexp2in1::AbstractArray{Float64, 2},
                                    ordered_vertices1::AbstractVector{Int},
                                    ordered_vertices2::AbstractVector{Int},
                                    num1in2::Int,
                                    num2in1::Int,
                                    ncomm::Int,
                                    tol::Float64)

    n = size(s₁, 1)
    IntVert = zeros(Float64, 0)
    ConvexExpIntVert = zeros(Float64, 0)
    Z = zeros(Float64, 2*n+2, 1)

    Indices = 1:n+1

    for label1 = 2^num1in2+1:2^(n+1) - 2
        b1 = binary_while(label1, n)
        num_vert1 = sum(b1)

        if num_vert1 > 1
            for label2 = 2^num2in1+1:2^(n+1)-2
                b2 = binary_while(label2, n)

                num_vert2 = sum(b2)

                if num_vert2 > 1
                    no_common_vert = true
                    #b1[1:ncomm] + b2[1:ncomm]   contains, 0,1 or 2
                    #If there is any vertex shared by the boundaries, a 2 will appear somewhere
                    if ncomm > 0
                        no_common_vert = maximum(b1[1:ncomm] + b2[1:ncomm]) < 2
                    end

                    if num_vert1 + num_vert2 <= n+2 && no_common_vert

                        if num_vert1 >= num_vert2
                            r = num_vert1
                            s = num_vert2
                            TargetVertices = view(s₂, :, :)
                            ReferenceBoundary = view(ordered_vertices1, Indices[findall(x->x!=0, b1)])
                            TargetBoundary = view(ordered_vertices2, Indices[findall(x->x!=0, b2)])
                            beta = view(convexexp2in1, :, :)
                            Gamma = view(convexexp2in1, setdiff(1:n+1,ReferenceBoundary), TargetBoundary)
                            Rank = rank(Gamma)
                            Rank0 = rank([Gamma; ones(1, s)])
                            no_vanishing_column = minimum(maximum(abs.(Gamma), dims=1))
                            Switch = 0
                        else
                            s = num_vert1
                            r = num_vert2
                            TargetVertices = view(s₁,:,:)
                            ReferenceBoundary = view(ordered_vertices2, Indices[findall(x->x!=0, b2)])
                            TargetBoundary = view(ordered_vertices1, Indices[findall(x->x!=0, b1)])
                            beta = view(convexexp1in2, :, :)
                            Gamma = view(convexexp1in2, setdiff(1:n+1,ReferenceBoundary), TargetBoundary)
                            Rank= rank(Gamma)
                            Rank0 = rank([Gamma;ones(1,s)])
                            no_vanishing_column = minimum(maximum(abs.(Gamma), dims=1))
                            Switch = 1
                        end

                        if Rank0-Rank == 1 && Rank == s-1 && no_vanishing_column > 0
                            lambda = QR(Gamma, tol)
                            alpha = [1 .- ones(Int,1,r-1)*view(beta,ReferenceBoundary[2:r],TargetBoundary) ;
                                    view(beta,ReferenceBoundary[2:r],TargetBoundary)]*lambda
                            alpha[abs.(alpha) .<= tol] .= 0
                            if min(minimum(alpha), minimum(lambda)) > 0 && sum(findall(x->x!=0, alpha))[1] > 1
                                # Filtering out non minimal boundaries
                                NewPoint = view(TargetVertices,:,TargetBoundary)*lambda
                                append!(IntVert,NewPoint)
                                Z .= 0
                                if Switch == 0
                                    Z[ReferenceBoundary] .= view(alpha,:)
                                    Z[TargetBoundary .+ (n + 1)] .= view(lambda,:)
                                else
                                    Z[ReferenceBoundary .+ (n + 1)] .= view(alpha,:)
                                    Z[TargetBoundary] .= view(lambda,:)
                                end

                                append!(ConvexExpIntVert,Z)
                            end
                        end
                    end
                end
            end
        end
    end

    return copy(transpose(reshape(IntVert, n, div(length(IntVert), n)))),
        copy(transpose(reshape(ConvexExpIntVert, 2*n+2, div(length(ConvexExpIntVert), 2*n+2))))
end

function binarydecomp(m::Int, n::Int)
    binary = zeros(Int,n+1)
    binary[n+1] = floor(Int,m/2^n)

    @inbounds for i = 1:n
        m = m-binary[n+1-i+1]*2^(n-i+1)
        binary[n+1-i] = floor(Int, m/2^(n-i))
    end

    return binary
end

function binary_while(m::Int, n::Int)
    binary = zeros(Int, n + 1)
    first1 = 0
    if m > 0
        first1 = floor(Int, log2(m))
        binary[first1+1] = 1
        m = m - 2^first1
        while m > 0
            first1 = floor(Int, log2(m))
            binary[first1 + 1] = 1
            m = m - 2^first1
        end
    end

    return binary
end


text:
import AbstractPlotting.MakieLayout.LColorbar

function AbstractPlotting.MakieLayout.LColorbar(parent, plots::Vector{<: AbstractPlot}; kwargs...)
    colormap = plots[1].colormap
    # compute colorrange
    min = minimum(p.colorrange[][1] for p in plots)
    max = maximum(p.colorrange[][2] for p in plots)
    colorrange = (min, max)

    for p in plots
        p.colorrange = colorrange
    end

    LColorbar(parent,
        colormap = plots[1].colormap,
        limits = colorrange;
        kwargs...)
end

max_stars_repo_name: greimel/AlgebraOfGraphics

max_issues_repo_name: greimel/AlgebraOfGraphics

max_forks_repo_name: greimel/AlgebraOfGraphics


text:
using Documenter, StatsBase, Statistics, Random, LinearAlgebra

# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

makedocs(
    sitename = "StatsBase.jl",
    modules = [StatsBase],
    pages = ["index.md",
             "weights.md",
             "scalarstats.md",
             "robust.md",
             "deviation.md",
             "cov.md",
             "counts.md",
             "ranking.md",
             "sampling.md",
             "empirical.md",
             "signalcorr.md",
             "multivariate.md",
             "misc.md",
             "statmodels.md",
             "transformations.md"],
    strict=true,
    checkdocs=:exports
)

deploydocs(
    repo = "github.com/JuliaStats/StatsBase.jl.git"
)


text:
module Sides

using Random

export
    Side,
    left,
    right,
    flipsign_if_right

@enum Side left right
flipsign_if_right(x::Number, side::Side) = ifelse(side == right, -x, x)
Base.rand(rng::AbstractRNG, ::Type{Side}) = ifelse(rand(rng, Bool), left, right)
Base.:(-)(side::Side) = ifelse(side == right, left, right)

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ de098278-8e74-11eb-13e3-49c1b6e06e7d
using Pkg; Pkg.activate("MLJ_env", shared=true)

# ╔═╡ ac463e7a-8b59-11eb-229e-db560e17c5f5
begin
	using Test
	using PlutoUI
end

# ╔═╡ 8c80e072-8b59-11eb-3c21-a18fe43c4536
md"""
## Databases and SQL

ref. from book **"Data Science from Scratch"**, Chap 24

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ e7373726-8b59-11eb-2a2b-b5138e4f5268
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ f5ee64b2-8b59-11eb-2751-0778efd589cd
html"""

"""

# ╔═╡ 81290d1c-8ce2-11eb-3340-337957fd81b7
html"""

"""

# ╔═╡ 8ff1bb20-8ce2-11eb-1de6-fd84daec8930
md"""
#### Create table and insert rows

"""

# ╔═╡ d3ee2138-8ce2-11eb-0b29-659a3be01512
md"""
We will represent our table in memory by a dictionary (not really space efficient).

Also, we will:
  - Support the primary key for our table (some table do not and will not have a primary key).
  - Allow for NULL value for any column of our table (but for the primary key).
"""

# ╔═╡ 6d759180-928b-11eb-1cc1-0593c7f3b0c2
begin
	const Row = Dict{Symbol, Any}
	const GRow = Union{Row, Vector{Pair{Symbol, Any}}, Vector{Any}}  ## Generic Row
	const WhereClause = Function
	const HavingClause = Function
	const S_N = Union{Symbol, Nothing}
	# const PK = Pair{S_N, DataType}
	const D_SF = Dict{Symbol, Function}
	const U_IN = Union{Int, Nothing}
	const U_SN = Union{String, Nothing}
	const UDT = Union{Union, DataType, Type};
end

# ╔═╡ 2e8684fc-2ba6-48ec-81b3-f1c31cdde43a
md"""
Let us define our custom type for a relational table:
"""

# ╔═╡ 53accaea-92d4-11eb-0bc2-3d2ee10f9bfb
mutable struct Table
	columns::Vector{Symbol}
	types::Vector{UDT}
	pkey::Union{Symbol, Nothing}
	rows::Vector{Row}

  function Table(col_types::Vector{Pair{Symbol, DT}};
               pkey::Pair=(:id => Int)) where DT <: UDT
    # vector of pairs keeps implicit order
    @assert length(col_types) ≥ 1
    #
    cols = map(((k, _)=p) -> k, col_types)
    types = map(((_, v)=p) -> v, col_types)
    cols, types = check_pk(cols, types, pkey)
    new(cols, types, pkey.first, Vector{Row}[])
  end

  function Table(cols::Vector{Symbol}, types::Vector{DT};
               pkey::Pair=(:id => Int)) where DT <: UDT
    cols, types = check_pk(cols, types, pkey)
    new(cols, types, pkey.first, Vector{Row}[])
  end
end

# ╔═╡ 1b1b8d8c-3803-43fa-b3af-7a58d75bb13b
md"""
And let us start developing our API and utilities:
"""

# ╔═╡ 6f9c58d4-92d1-11eb-2c09-cb1ea5afcd6d
begin
	import Base: length
	##
	## API for Table
	##

	id(self::Table) = self.pkey

	length(self::Table) = length(self.rows)

	function gen_pkey_value(self::Table)
	  """
	  Assuming pkey is instance of a numeric type
	  """
	  length(self.rows) == 0 && (return 1)
	  (map(r -> r[id(self)], self.rows) |> maximum)  + 1
	end

	## utility
	function check_pk(cols, types, pkey)
	  if pkey.first !== nothing
		@assert pkey[2] <: Integer "Expecting an Integer, got: $(pkey[2])"
	  end

	  if pkey.first !== nothing && pkey.first ∉ cols
		cols  = [pkey[1], cols...]
		types = [pkey[2], types...]
	  end
	  (cols, types)
	end
end

# ╔═╡ 040a7354-92d2-11eb-2825-7bac34b9fdf9
YaUsers = Table([:name => String, :num_friends => Int];
	pkey=(:user_id => Int32))

# ╔═╡ 17d214a7-be15-4236-a614-3127b535ef66
AltUsers = Table([:name => String, :num_friends => Int]) ## Default pkey

# ╔═╡ 85211d8a-be50-454d-a6c5-cd820a20f8d4
## Default pkey and NULL values allowed
AltUsers₂ = Table([:name, :num_friends], [U_SN, U_IN])

# ╔═╡ adf67244-92d4-11eb-3004-41e62e906e32
begin
	##
	## Convention: by default all tables have a primary key called :id,
	## unless explicitly stated otherwise, something like:
	##   :pkey => (:user_id, Int)
	##
	Users = Table([:name => U_SN, :num_friends => U_IN];
		pkey=(:user_id => Int))
	@assert length(Users) == 0
	Users;
end

# ╔═╡ 650237ee-928e-11eb-0e43-17ab98a0cc11
begin
	##
	## API for Table (cont'ed)
	##

	function insert(self::Table, row::Vector{Pair{Symbol, DT}}) where DT <: Any
	  """
	  row is something like [:name => "FooBar", :num_friends => 3, id(Users) => 3]
	  """
	  res = filter(p -> p.first == id(self), row)
	  has_pkey = length(res) > 0
	  ##
	  ## NULL values permitted
	  #length(row) ≠ length(self.types) && has_pkey &&
	  #  throw(ArgumentError("Mismatch with expected number of columns"))
	  ##
	  check_value_dtype(self, row)
	  ## find the pair <pkey, value> pkey is given by id(self)
	  row_ids = filter(p -> p.first == id(self), row)
	  if length(row_ids) == 0
		pk_id =  gen_pkey_value(self)
		row = [row..., id(self) => pk_id]
	  else
		id_val = row_ids[1][2]  ## value assoc with this pkey
		row_already_inserted(self, id_val) && (return nothing)
	  end
	  push!(self.rows, Dict(row...))
	end

	function insert(self::Table, row::Vector{Any})
	  """
	  row is something like: [9, "Foo", 2 ] - just pure values
	  """
	  length(row) ≠ length(self.types) && id(self) ∈ map(p -> p[1], row) &&
		throw(ArgumentError("Mismatch with expected number of columns "))
	  #
	  check_value_dtype(self, row)
	  ## row[1] is the value assoc. with pkey
	  row_already_inserted(self, row[1]) && (return nothing)
	  push!(self.rows, Dict(zip(self.columns, row)))
	end

	function insert(self::Table, row::Row)
	  """
	  row is a  Dict, like: insert(Users, Dict(:name => "Foo", :num_friends => 5))
	  """
	  length(row) ≠ length(self.types) && haskey(row, id(self)) &&
		throw(ArgumentError("Mismatch with expected number of columns"))
	  check_value_dtype(self, row)
	  if haskey(row, id(self))
		row_already_inserted(self, row[id(self)]) && (return nothing)
	  else
		row[id(self)] = gen_pkey_value(self)
	  end
	  push!(self.rows, row)
	end

	insert(self::Table, rows::Vector{T}) where T <: GRow = insert.(Ref(self), rows)

	function coltype(self::Table, colname::Symbol)::UDT
	  ix = findfirst(col -> col == colname, self.columns)
	  ix === nothing && throw(ArgumentError("column $(colname) inexistent"))
	  self.types[ix]
	end

	# function Base.show(io::IO, self::Table)
	# 	s = "num. records: $(length(self)):\n"
	# 	for rd ∈ self.rows
	# 		s₁ = []
	# 		for (k, v) ∈ rd
	# 			push!(s₁, "$(k): $(v)")
	# 		end
	# 		s = string(s, "  <", join(s₁, ", "), ">\n")
	# 	end
	# 	print(io, "$(s)")
	# end

	##
	## Internal checkers
	##

	function check_value_dtype(self::Table, values::Vector{DT}) where DT <: Any
	  dtypes = self.types
	  for (v, dt) ∈ zip(values, dtypes)
		!(typeof(v) <: dt) && v !== nothing &&
		  throw("Expected type for <$(v)>: $(dt), got: $(typeof(v))")
	  end
	end

	function check_value_dtype(self::Table,
							   row::Union{Row, Vector{Pair{Symbol, DT}}}) where DT <: Any
	  (cols, types) = self.columns, self.types
	  col_type = Dict(zip(cols, types))
	  for (k, v) ∈ row
		(haskey(col_type, k) && typeof(v) <: col_type[k]) ||
		  throw("Expected type/2: $(col_type[k]), got: $(typeof(v))")
	  end
	end

	function row_already_inserted(self::Table, id_val::Any)::Bool
		id(self) === nothing && return false  ## no pkey => ignore check
		#
		res = filter(r -> r[id(self)] == id_val, self.rows)
		length(res) > 0    ## row already inserted id > 0
	end

end

# ╔═╡ d146a8ac-92cc-11eb-29b1-bb065a468477
md"""
Insert using an array/vector:
"""

# ╔═╡ 13a0940e-8ce4-11eb-231c-f331a607203c
begin
	insert(Users, [0, "Hero", 0])
	@assert length(Users) ≥ 1
end

# ╔═╡ e3458d34-92cc-11eb-2426-79f229abd908
md"""
Insert using a vector of pairs:
"""

# ╔═╡ 96c2a668-9290-11eb-074d-19ffeb68eba6
begin
	insert(Users, [id(Users) => 1, :name => "Dunn", :num_friends => 2])
	@assert length(Users) ≥ 2
end

# ╔═╡ f1357972-92cc-11eb-1f73-43aae5a43aa5
md"""
Insert using a dictionary
"""

# ╔═╡ 5a3a8fdc-929d-11eb-25e3-e309dffdb455
begin
	## Insert with a Dict
	insert(Users, Dict(id(Users) => 2, :name => "Sue", :num_friends => 3))
	@assert length(Users) ≥ 3
end

# ╔═╡ 0292e48e-92cd-11eb-1882-81a2faaad136
md"""
Insert using a dictionary, no primary key value specified:
"""

# ╔═╡ 8545d17c-929e-11eb-0989-ffd4749bbb25
begin
	## Insert with a Dict no pkey/1
	insert(Users, Dict(:name => "Ayumi", :num_friends => 5))
	@assert length(Users) ≥ 4
end

# ╔═╡ 95d2b19e-92a4-11eb-0f91-e343c9317983
begin
	## Insert with a Dict no pkey/2
	insert(Users, [:name => "PasMas", :num_friends => 5])
	@assert length(Users) ≥ 5
end

# ╔═╡ 8a945abc-012f-4ace-a67a-b7e0f520346b
begin
	## Insert with a Dict no pkey/2 + null value (encoded as nothing)
	# insert(Users,
	# 	[:name => "Foobar", :num_friends => nothing])
	insert(Users,
		[:name => "Foobar", :num_friends => nothing, :user_id => 17])
end

# ╔═╡ 28ba0bf6-92cd-11eb-0669-75eadb768518
md"""
Insert using a collection (vector of rows):
"""

# ╔═╡ d4562f7e-9290-11eb-2d8c-952f2e0edfed
begin
	## Insert a collection (Vector of) of Records
	insert(Users, [
		[5, "Chi", 3],
		[6, "Thor", 3],
		[7, "Clive", 2],
		[8, "Devin", 2],
		[9, "Kate", 2],
		[10, "Kaze", nothing]
	])
	@assert length(Users) ≥ 10
end

# ╔═╡ f877c3e6-929c-11eb-00cd-c15568f99627
Users

# ╔═╡ d5657a88-92c4-11eb-2a4f-c7dcaa97bcf6
begin
	@test_throws ArgumentError  coltype(Users, :foobar)
	@test coltype(Users, :num_friends) == U_IN
	@test coltype(Users, :name) == U_SN
end

# ╔═╡ 80ddf0f4-8ce2-11eb-3046-331119a0dc9b
html"""

"""

# ╔═╡ 80a71408-8ce2-11eb-3978-75a4a2df9116
md"""
#### Update
"""

# ╔═╡ d1ed1e36-8ce5-11eb-34ee-4db60fb0db8a
md"""
The key features for an update:
  - what table,
  - which fields/rows
  - what their new values will be
"""

# ╔═╡ 4ccdb1ce-92c2-11eb-18dd-c510e88474c6
##
## API (cont'ed)
##
function update(self::Table, updates::Row, pred::WhereClause=row -> true)
	## 1 - Make sure key/columns and values are consistent with column types
	for (col, val) ∈ updates
		col ∉ self.columns && throw(ArgumentError("invalid column $(col)"))
		col_type = coltype(self, col)
		!(typeof(val) <: col_type) && val !== nothing &&
			throw(ArgumentError("Expected $(col_type), got $(typeof(val))"))
	end

	## 2 - OK, update
	rows_to_update = filter(((ix, r)=t) -> pred(r), collect(enumerate(self.rows)))
	for (ix, row) ∈ rows_to_update
		# self.rows[ix] = Row(row..., updates...)
		for (k, v) ∈ updates
			self.rows[ix][k] = v
		end
	end
	nothing
end

# ╔═╡ 2af6dbec-92c6-11eb-1155-596ed44b04f8
begin
	update(Users, Dict{Symbol, Any}(:num_friends => 7),
		row -> row[:name] == "Ayumi")
	Users
end

# ╔═╡ 2abb0114-92c6-11eb-01d5-cd34aa6426b2
begin
	update(Users, Dict{Symbol, Any}(:num_friends => 4), row -> row[:num_friends] == 2)
	Users
end

# ╔═╡ bc1fb8d2-92ca-11eb-14e0-c9586461bc1e
begin
	update(Users, Dict{Symbol, Any}(:num_friends => 1))
	Users
end

# ╔═╡ d2b198b6-8ce2-11eb-170e-0f17904c9f2c
html"""

"""

# ╔═╡ d37f37b8-8ce8-11eb-2c00-3f98ca407f41
md"""
#### Delete
"""

# ╔═╡ 5bf7729e-8dfd-11eb-070f-9b7ec0746bd7
##
## API (cont'ed)
##
function delete(self::Table, pred::WhereClause=row -> true)|
	self.rows = [
		row for row ∈ self.rows if !pred(row)
	]
	nothing
end

# ╔═╡ 394093fa-8cea-11eb-0071-eff3045a012b
begin
	n = length(Users)
	delete(Users, row -> row[id(Users)] == 1)
	@assert length(Users) == n - 1
end

# ╔═╡ 435c4444-8dfd-11eb-24e8-5543f905f199
begin
	delete(Users)
	@assert length(Users) == 0
end

# ╔═╡ e95983ba-8ceb-11eb-38fd-ed92cdcf754c
html"""

"""

# ╔═╡ d3a749a2-8cec-11eb-1f06-b568f244b576
md"""
#### Select
"""

# ╔═╡ 31d3226a-8cf4-11eb-0897-39989ba76b58
md"""
We will give our Table struct a select method that returns a new Table . The
method accepts two optional arguments:
  - `keep_cols` which specifies the names of the columns we want to
keep in the result. if none supply, the result contains all the columns, and

  - `add_cols` is a dictionary whose keys are new column names and whose values are functions specifying how to compute the values of the new columns.
"""

# ╔═╡ 3ec9f5d5-6d6d-4586-b68b-84e3ee5aee27
## Utility function
##
function create_res_table(self::Table, new_cols, new_types;
                          pred::Function=_a -> true)
  ## keep pkey or not...
  pred(self) ?
	Table(new_cols, new_types; pkey=id(self) => coltype(self, id(self))) :
    Table(new_cols, new_types; pkey=nothing => Nothing)
end

# ╔═╡ a56253c8-92cd-11eb-31b9-51e1fd5ae026
function select(self::Table;
		keep_cols=Vector{Symbol}[], add_cols=D_SF())::Table
	##
	length(keep_cols) == 0 && (keep_cols = self.columns)

	## New column names and types
	new_cols = [keep_cols..., collect(keys(add_cols))...]
	keep_types = [coltype(self, col) for col in keep_cols]

	## collect the rows for result table
	n_rows = Vector{Any}[]
	add_types = Any[]
	for (ix, row) ∈ enumerate(self.rows)
		n_row = Any[row[col] for col ∈ keep_cols]
		## as we process the first row, we can get the return type...
		## ...of each function defined in add_cols
		## What if no row to process => tag with Any
		for (_col_name, fn) ∈ add_cols
			r = fn(row)
			ix == 1 && push!(add_types, typeof(r))
			push!(n_row, r)
		end
		push!(n_rows, n_row)
	end

	if length(add_cols) > 0 && length(add_types) == 0
		## add as many Any as column in add_cols
		push!(add_types, repeat(Any[Any], inner=length(add_cols)))
	end

	## Create result table
	new_types = UDT[keep_types..., add_types...]
	@assert(length(new_cols) == length(new_types),
		"length(new_cols) == length(new_types)")

	n_table = create_res_table(self, new_cols, new_types;
		pred=s -> id(s) ∈ keep_cols)

	insert(n_table, n_rows)
	n_table
end

# ╔═╡ e6348294-92cb-11eb-3bfd-09d80933b33a
begin
	delete(Users)
	insert(Users, [[0, "Hero", 0],
			[1, "Dunn", 2],
			[2, "Sue", 3],
			[3, "Chi", 3],
			[4, "Thor", 3],
			[5, "Clive", 2],
			[6, "Hicks", 3],
			[7, "Devin", 2],
			[8, "Kate", 2]
	])
	insert(Users, Dict(:name => "Ayumi", :num_friends => 5))
	insert(Users, [:name => "PasMas", :num_friends => 5])
end

# ╔═╡ 96386ab8-93ef-44c7-a524-1246379d7103
Table([:user_id, :name, :num_friends], UDT[Int64, Union{Nothing, String}, Union{Nothing, Int64}])

# ╔═╡ a5441766-92cd-11eb-1c8e-edcf3db0022e
begin
	## SELECT * FROM Users;
	n₀ = length(Users)
	all_users = select(Users)
	@assert length(all_users) == n₀
end

# ╔═╡ 3bfad78b-2ebd-4067-b306-b687f5a92a10
begin
	## SELECT id FROM Users;
	user_ids = select(Users, keep_cols=[:name])
end

# ╔═╡ cfed19ac-94e2-4c84-85b4-6134765cce0b
md"""
##### Limit
"""

# ╔═╡ 386a0591-92ce-4267-b1e6-73e293eb727c
##
## API (cont'ed)
##
function limit(self::Table, num_rows::Int=5)::Table
	"""
	Only the first num_rows are returned
	"""
	@assert 1 ≤ num_rows ≤ length(self.rows) "1 ≤ $(num_rows) ≤ $(length(self.rows))"

 	n_table = create_res_table(self, self.columns, self.types;
		pred=s -> id(s) !== nothing)

	## NOTE: mark vector as GRow type
	rows = Vector{Any}[]
	for row ∈ self.rows[1:num_rows]
		push!(rows, Any[v for (_, v) ∈ row])
	end
	insert(n_table, rows)
	n_table
end

# ╔═╡ 53c15148-1b04-4d85-aac6-874a6f750a00
user_names₁ = select(Users, keep_cols=[:name]) |> limit

# ╔═╡ 1ac32ed0-6e47-4606-8b70-524159660d93
user_names₂ = select(Users, keep_cols=[:name]) |>
  u -> limit(u, 2)

# ╔═╡ be425e03-c7b4-401d-9d29-fa2b07618ee5
md"""
###### Where
"""

# ╔═╡ e99eab12-08d9-4955-ba10-2eb46f74bd85
##
## API (cont'ed)
##
function where(self::Table, pred::WhereClause=row -> true)::Table
	"""
	Only the rows that satisfy pred are returned
	"""
	n_table = create_res_table(self, self.columns, self.types;
		pred=s -> id(s) !== nothing)

	n_rows = Vector{Any}[]
	for row ∈ self.rows
		pred(row) && (push!(n_rows, [row[col] for col ∈ n_table.columns]))
	end

	insert(n_table, n_rows)
	n_table
end

# ╔═╡ b08b9470-5d53-4420-bef6-ef1c0bb4414c
begin
	dunn_ids = where(Users, row -> row[:name] == "Dunn") |>
		u -> select(u, keep_cols=[id(Users)])
	# @assert length(dunn_ids) == 1
end

# ╔═╡ 955fe2a1-22b3-44db-b23d-d595d58bac3d
begin
	amp_ids = where(Users, row -> row[:name][1] ∈ ['A', 'P']) |>
		u -> select(u, keep_cols=[id(Users), :name])
end

# ╔═╡ 964507de-c48d-4f7c-ab05-9b7739ecd5f3
begin
	ncol = length(Users.columns)
	function name_len_fn(row)::Int
		length(row[:name])
	end

	name_lengths = select(Users;
		add_cols=D_SF(:name_length => name_len_fn))
end

# ╔═╡ 13c464eb-2a59-488b-b0ab-922ecda91bbd
begin
	name_lengths₂ = select(Users;
		add_cols=Dict{Symbol, Function}(
			:name_ini => row -> string(row[:name][1]),
			:len_name => row -> length(row)
		)
	)
end

# ╔═╡ 95ae9db4-f9ba-4fde-a021-0fc952550627
md"""
##### Aside on introspection
"""

# ╔═╡ 17b2e439-4f70-4db6-a01f-638ae2f88b6f
begin
	## with explicit return type
	##
	str_fn = """
function name_len_fn(row, bar, vargs...)::UInt16
	length(row[:name])
end
"""
	expr = Meta.parse(str_fn)
end

# ╔═╡ a57b3859-d026-4325-a5b2-7d4e7b4db408
begin
	## w/o explicit return type
	##
	str_fn₂ = """
function name_len_fn(row1, bar1; foo1=10)
	length(row[:name])
end
"""
	expr₂ = Meta.parse(str_fn₂)
end

# ╔═╡ 8e25eb57-b6f0-4c61-b657-4d815f1eb31f
with_terminal() do
	dump(expr)
end

# ╔═╡ 8ea2ba50-ad35-40b2-ab2d-1c34c9d90bd9
expr.args[1], expr₂.args[1]

# ╔═╡ b34d8aaa-7696-452f-9fd7-1002aa771547
## get return type
expr.args[1].args, length(expr.args[1].args), expr.args[1].args[end]

# ╔═╡ 4ff7e8e9-354e-4c99-b6f8-aedbf7ce1f3f
expr₂.args[1].args, length(expr₂.args[1].args)

# ╔═╡ 3a438295-3d30-41cc-95bd-db3b93e62d12
function get_fn_retype(str_fn::String)::DataType
	"""
	Return output type of a user function if such type is available...
	...which is the case iff expr.args[1].args has length of 2
	Otherwise fallback to Any
	"""
	expr = Meta.parse(str_fn)
	length(expr.args[1].args) == 2 ? eval(expr.args[1].args[end]) : Any
end

# ╔═╡ 9b24fead-ccfd-43ed-afa8-ed5ef6c4cf8f
get_fn_retype(str_fn), get_fn_retype(str_fn₂)

# ╔═╡ 6e7e7896-8cf8-11eb-062e-99492ec8cff8
html"""

"""

# ╔═╡ 70f707c0-8cf6-11eb-11c2-73d7b28f7a0c
md"""
#### Group by
"""

# ╔═╡ 70d84470-8cf6-11eb-23b6-c7ba506d8552
md"""
This function will take a list of columns we want to group by and a dictionary of aggragation functions plus an optional predicate parameter (having clause) that can operate on multiple rows.
"""

# ╔═╡ c451451a-8cf7-11eb-183b-e358c9d618e0
##
## API (cont'ed)
##

function group_by(self::Table; group_by_cols::Vector{Symbol}, agg::D_SF,
	having=HavingClause=row_gp -> true)::Table
	grouped_rows = Dict{}()

	## 1 - Populate groups
	for row ∈ self.rows
		key = Tuple(row[col] for col ∈ group_by_cols)
		ary = get(grouped_rows, key, Row[])
		push!(ary, row)
		grouped_rows[key] = ary
	end

	## 2 - populate rows
	n_rows = Vector{Any}[]
	agg_types = Any[]
	for (ix, (key, rows)) ∈ enumerate(grouped_rows)
		agg_row = Any[]
		if having(rows)
			agg_row = Any[agg_row..., key...]  ## Keep Any[]
			for (_, agg_fn) ∈ agg
				r = agg_fn(rows)
				ix == 1 && push!(agg_types, typeof(r))
				push!(agg_row, r)
			end
		end
		length(agg_row) > 0 && (push!(n_rows, agg_row))
	end

	## 3 - res. table consists of group_by columns and aggregates
	new_cols = [group_by_cols..., collect(keys(agg))...] |> unique
	gp_by_types = [coltype(self, col) for col ∈ group_by_cols]
	new_types = [gp_by_types..., agg_types...]

	n_table = create_res_table(self, new_cols, new_types;
		pred=s -> id(s) ∈ new_cols)
	insert(n_table, n_rows)
	n_table
end

# ╔═╡ 30bdc788-8d01-11eb-2f8f-9fd79593ddb8
begin
	"""
	-- find number of users and smallest user_id for each possible name length:
	SELECT LENGTH(name) as name_length, MIN(user_id) AS min_user_id,
		COUNT(*) AS num_users
	FROM users
	GROUP BY LENGTH(name);
	"""
	min_user_id = rows -> minimum(row[id(Users)] for row ∈ rows)
	num_rows = rows -> length(rows)

	stats_by_len = select(Users, add_cols=D_SF(:name_length => name_len_fn)) |>
		u -> group_by(u,
				group_by_cols=[:name_length],
				agg=D_SF(:min_user_id => min_user_id, :num_users => num_rows))
	#
end

# ╔═╡ c5bac3ae-c63b-44ba-9b10-75144410f456
Users

# ╔═╡ e696c2dd-b250-418c-9161-9fd87bfc378e
begin
	"""
	-- average number of friends for users whose names start with specific letters
	-- but see only the results for letters whose corresponding average is greater
	-- than 1
	SELECT SUBSTR(name, 1, 1) AS first_letter,
		AVG(num_friends) AS avg_num_friends
	FROM users
	GROUP BY SUBSTR(name, 1, 1)
	HAVING AVG(num_friends) > 1;
	"""
	first_letter_fn = row -> row[:name] !== nothing ? string(row[:name][1]) : ""

	avg_num_friends_fn =
		rows -> sum([row[:num_friends] for row in rows]) / length(rows)

	enough_friends_fn = rows -> avg_num_friends_fn(rows) > 1.

	avg_friends_by_letter =
		select(Users, add_cols=D_SF(:first_letter => first_letter_fn)) |>
		u -> group_by(u, group_by_cols=[:first_letter],
				agg=D_SF(:avg_num_friends => avg_num_friends_fn),
				having=enough_friends_fn)
end

# ╔═╡ 84946000-8d02-11eb-3160-571efee8fb0b
html"""

"""

# ╔═╡ d05f1c40-8d17-11eb-2724-4d989c4c6b92
md"""
#### Order by
"""

# ╔═╡ 848fd812-8d02-11eb-01d6-75965b08bcc5
##
## API (cont'ed)
##
function order_by(self::Table, order::Function)::Table
	n_table = select(self)
	sort!(n_table.rows, by=order)
	n_table
end

# ╔═╡ 6cde1330-e1bb-419e-87f3-cf73c0e11e3d
friendliest_letters =  avg_friends_by_letter |>
	u -> order_by(u, row -> -row[:avg_num_friends])
# NOTE: -row(...) to reverse the sort

# ╔═╡ 7d35f85e-8e07-11eb-228c-81d28a444b52
friendliest_letters₂ =  avg_friends_by_letter |>
	u -> order_by(u, row -> -row[:avg_num_friends]) |>
	u -> limit(u, 3)
# NOTE: -row(...) to reverse the sort

# ╔═╡ 40d46802-8e0f-11eb-2e47-53096e24dbd8
html"""

"""

# ╔═╡ 80dff352-8d02-11eb-06b1-5f5e325046f5
md"""
#### Join
"""

# ╔═╡ f1af4256-8e02-11eb-0238-2515d39a89cd
##
## API (cont'ed)
##
function join(self::Table, otable::Table; left_join=false)::Table
	join_on_cols = [c for c ∈ self.columns if c ∈ otable.columns]

	add_cols = [c for c ∈ otable.columns if c ∉ join_on_cols && c != id(otable)]
	new_cols = [self.columns..., add_cols...]
	new_types = UDT[self.types..., coltype.(Ref(otable), add_cols)...]

	n_table = Table(new_cols, new_types; pkey=nothing => Nothing)

	n_rows = Vector{Any}[]
	for row ∈ self.rows
		is_join = orow -> all(orow[c] == row[c] for c ∈ join_on_cols)

		o_rows = where(otable, is_join).rows
		for o_row ∈ o_rows
			push!(n_rows, Any[[row[c] for c ∈ self.columns]...,
					[o_row[c] for c ∈ add_cols]...])
		end

		if left_join && length(o_rows) == 0
			push!(n_rows, Any[[row[c] for c ∈ self.columns]...,
					[nothing for _ ∈ add_cols]...])
		end
	end
	@show "join/4", n_rows, length(n_rows)
	insert(n_table, n_rows)
	n_table
end

# ╔═╡ 26571096-f30c-476a-ad0e-ef133ba2562f
begin
	User_Interests = Table([:user_id => Int, :interest => String]) ## Default pkey
	insert(User_Interests, [
			[1, 0, "SQL"],
			[2, 0, "NoSQL"],
			[3, 2, "SQL"],
			[4, 2, "MySQL"],
			[5, 7, "PostgreSQL"],
			[6, 7, "SQL"]
	])
end

# ╔═╡ a1a91e04-ffdb-49b1-a999-e1cf1a730cb4
begin
	sql_users = join(Users, User_Interests) |>
		u -> where(u, r -> r[:interest] == "SQL") |>
		u -> select(u, keep_cols=[:name])
end

# ╔═╡ 43a52692-8e10-11eb-2043-8f0be195f58a
begin
	sql_users₂ = join(Users, User_Interests; left_join=true) |>
		u -> where(u, r -> r[:interest] !== nothing &&
							occursin(r"SQL\z", r[:interest])) |>
		u -> select(u, keep_cols=[:name])
	@test length(sql_users₂.rows) == 6
end

# ╔═╡ Cell order:
# ╟─8c80e072-8b59-11eb-3c21-a18fe43c4536
# ╠═de098278-8e74-11eb-13e3-49c1b6e06e7d
# ╠═ac463e7a-8b59-11eb-229e-db560e17c5f5
# ╟─e7373726-8b59-11eb-2a2b-b5138e4f5268
# ╟─f5ee64b2-8b59-11eb-2751-0778efd589cd
# ╟─81290d1c-8ce2-11eb-3340-337957fd81b7
# ╟─8ff1bb20-8ce2-11eb-1de6-fd84daec8930
# ╟─d3ee2138-8ce2-11eb-0b29-659a3be01512
# ╠═6d759180-928b-11eb-1cc1-0593c7f3b0c2
# ╟─2e8684fc-2ba6-48ec-81b3-f1c31cdde43a
# ╠═53accaea-92d4-11eb-0bc2-3d2ee10f9bfb
# ╟─1b1b8d8c-3803-43fa-b3af-7a58d75bb13b
# ╠═6f9c58d4-92d1-11eb-2c09-cb1ea5afcd6d
# ╠═040a7354-92d2-11eb-2825-7bac34b9fdf9
# ╠═17d214a7-be15-4236-a614-3127b535ef66
# ╠═85211d8a-be50-454d-a6c5-cd820a20f8d4
# ╠═adf67244-92d4-11eb-3004-41e62e906e32
# ╠═650237ee-928e-11eb-0e43-17ab98a0cc11
# ╟─d146a8ac-92cc-11eb-29b1-bb065a468477
# ╠═13a0940e-8ce4-11eb-231c-f331a607203c
# ╟─e3458d34-92cc-11eb-2426-79f229abd908
# ╠═96c2a668-9290-11eb-074d-19ffeb68eba6
# ╟─f1357972-92cc-11eb-1f73-43aae5a43aa5
# ╠═5a3a8fdc-929d-11eb-25e3-e309dffdb455
# ╟─0292e48e-92cd-11eb-1882-81a2faaad136
# ╠═8545d17c-929e-11eb-0989-ffd4749bbb25
# ╠═95d2b19e-92a4-11eb-0f91-e343c9317983
# ╠═8a945abc-012f-4ace-a67a-b7e0f520346b
# ╟─28ba0bf6-92cd-11eb-0669-75eadb768518
# ╠═d4562f7e-9290-11eb-2d8c-952f2e0edfed
# ╠═f877c3e6-929c-11eb-00cd-c15568f99627
# ╠═d5657a88-92c4-11eb-2a4f-c7dcaa97bcf6
# ╟─80ddf0f4-8ce2-11eb-3046-331119a0dc9b
# ╟─80a71408-8ce2-11eb-3978-75a4a2df9116
# ╟─d1ed1e36-8ce5-11eb-34ee-4db60fb0db8a
# ╠═4ccdb1ce-92c2-11eb-18dd-c510e88474c6
# ╠═2af6dbec-92c6-11eb-1155-596ed44b04f8
# ╠═2abb0114-92c6-11eb-01d5-cd34aa6426b2
# ╠═bc1fb8d2-92ca-11eb-14e0-c9586461bc1e
# ╟─d2b198b6-8ce2-11eb-170e-0f17904c9f2c
# ╟─d37f37b8-8ce8-11eb-2c00-3f98ca407f41
# ╠═5bf7729e-8dfd-11eb-070f-9b7ec0746bd7
# ╠═394093fa-8cea-11eb-0071-eff3045a012b
# ╠═435c4444-8dfd-11eb-24e8-5543f905f199
# ╟─e95983ba-8ceb-11eb-38fd-ed92cdcf754c
# ╟─d3a749a2-8cec-11eb-1f06-b568f244b576
# ╟─31d3226a-8cf4-11eb-0897-39989ba76b58
# ╠═3ec9f5d5-6d6d-4586-b68b-84e3ee5aee27
# ╠═a56253c8-92cd-11eb-31b9-51e1fd5ae026
# ╠═e6348294-92cb-11eb-3bfd-09d80933b33a
# ╠═96386ab8-93ef-44c7-a524-1246379d7103
# ╠═a5441766-92cd-11eb-1c8e-edcf3db0022e
# ╠═3bfad78b-2ebd-4067-b306-b687f5a92a10
# ╟─cfed19ac-94e2-4c84-85b4-6134765cce0b
# ╠═386a0591-92ce-4267-b1e6-73e293eb727c
# ╠═53c15148-1b04-4d85-aac6-874a6f750a00
# ╠═1ac32ed0-6e47-4606-8b70-524159660d93
# ╟─be425e03-c7b4-401d-9d29-fa2b07618ee5
# ╠═e99eab12-08d9-4955-ba10-2eb46f74bd85
# ╠═b08b9470-5d53-4420-bef6-ef1c0bb4414c
# ╠═955fe2a1-22b3-44db-b23d-d595d58bac3d
# ╠═964507de-c48d-4f7c-ab05-9b7739ecd5f3
# ╠═13c464eb-2a59-488b-b0ab-922ecda91bbd
# ╟─95ae9db4-f9ba-4fde-a021-0fc952550627
# ╠═17b2e439-4f70-4db6-a01f-638ae2f88b6f
# ╠═a57b3859-d026-4325-a5b2-7d4e7b4db408
# ╠═8e25eb57-b6f0-4c61-b657-4d815f1eb31f
# ╠═8ea2ba50-ad35-40b2-ab2d-1c34c9d90bd9
# ╠═b34d8aaa-7696-452f-9fd7-1002aa771547
# ╠═4ff7e8e9-354e-4c99-b6f8-aedbf7ce1f3f
# ╠═3a438295-3d30-41cc-95bd-db3b93e62d12
# ╠═9b24fead-ccfd-43ed-afa8-ed5ef6c4cf8f
# ╟─6e7e7896-8cf8-11eb-062e-99492ec8cff8
# ╟─70f707c0-8cf6-11eb-11c2-73d7b28f7a0c
# ╟─70d84470-8cf6-11eb-23b6-c7ba506d8552
# ╠═c451451a-8cf7-11eb-183b-e358c9d618e0
# ╠═30bdc788-8d01-11eb-2f8f-9fd79593ddb8
# ╠═c5bac3ae-c63b-44ba-9b10-75144410f456
# ╠═e696c2dd-b250-418c-9161-9fd87bfc378e
# ╟─84946000-8d02-11eb-3160-571efee8fb0b
# ╟─d05f1c40-8d17-11eb-2724-4d989c4c6b92
# ╠═848fd812-8d02-11eb-01d6-75965b08bcc5
# ╠═6cde1330-e1bb-419e-87f3-cf73c0e11e3d
# ╠═7d35f85e-8e07-11eb-228c-81d28a444b52
# ╟─40d46802-8e0f-11eb-2e47-53096e24dbd8
# ╟─80dff352-8d02-11eb-06b1-5f5e325046f5
# ╠═f1af4256-8e02-11eb-0238-2515d39a89cd
# ╠═26571096-f30c-476a-ad0e-ef133ba2562f
# ╠═a1a91e04-ffdb-49b1-a999-e1cf1a730cb4
# ╠═43a52692-8e10-11eb-2043-8f0be195f58a

max_stars_repo_path: Data Science From Scratch/24_DB_SQL.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/24_DB_SQL.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/24_DB_SQL.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
# This file was generated, do not modify it.

using MLJ, DataFrames, Random
@load RidgeRegressor pkg=MultivariateStats

Random.seed!(5) # for reproducibility
x1 = rand(300)
x2 = rand(300)
x3 = rand(300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(300))
X = DataFrame(x1=x1, x2=x2, x3=x3)

test, train = partition(eachindex(y), 0.8);

Xs = source(X)
ys = source(y, kind=:target)

std_model = Standardizer()
stand = machine(std_model, Xs)
W = transform(stand, Xs)

box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, ys)
z = transform(box, ys)

ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

ŷ = inverse_transform(box, ẑ)

@from_network CompositeModel(std=std_model, box=box_model,
                             ridge=ridge_model) <= ŷ;

cm = machine(CompositeModel(), X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8),
                measure=rms)
round(res.measurement[1], sigdigits=3)

mutable struct CompositeModel2 <: DeterministicNetwork
    std_model::Standardizer
    box_model::UnivariateBoxCoxTransformer
    ridge_model::RidgeRegressor
end

function MLJ.fit(m::CompositeModel2, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y, kind=:target)
    W = transform(machine(m.std_model, Xs), Xs)
    box = machine(m.box_model, ys)
    z = transform(box, ys)
    ẑ = predict(machine(m.ridge_model, W, z), W)
    ŷ = inverse_transform(box, ẑ)
    fit!(ŷ, verbosity=0)
    return fitresults(ŷ)
end

mdl = CompositeModel2(Standardizer(), UnivariateBoxCoxTransformer(),
                      RidgeRegressor(lambda=0.1))
cm = machine(mdl, X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8), measure=rms)
round(res.measurement[1], sigdigits=3)


text:
using Test
using DataDeps
using DelimitedFiles

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

@testset "Pi" begin
    register(DataDep(
     "Pi",
     "There is no real reason to download Pi, unlike say lists of prime numbers, it is always faster to compute than it is to download. No matter how many digits you want.",
     "https://www.angio.net/pi/digits/10000.txt",
     sha2_256
    ))

    pi_string = read(datadep"Pi/10000.txt", String)
    @test parse(Float64, pi_string) ≈ π
    @test parse(BigFloat, pi_string) ≈ π

end

@testset "TrumpTweets" begin
    # This tests demostrates how the `post_fetch_method` can be used to synthesize new files

    register(DataDep("TrumpTweets",
    """
    Tweets from 538's article:
    [The World’s Favorite Donald Trump Tweets](https://fivethirtyeight.com/features/the-worlds-favorite-donald-trump-tweets/)

    Includes a filtered view that is the tweats filtered to remove any tweets that @mention anyone,
    so no conversations etc, just announcements of opinions/thoughts.

    Used under Creative Commons Attribution 4.0 International License.
    """,
    "https://raw.githack.com/fivethirtyeight/data/master/trump-twitter/realDonaldTrump_poll_tweets.csv",
    "5a63b6cb2503a20517b5d41bd73e821ffbfdddd5cdc1977a547f1c925790bb15",
    post_fetch_method = function(in_fn) # Multiline anon function.
        out_fn = "nonmentions_"*basename(in_fn)
        print(out_fn)
        open(out_fn, "w") do out_fh
            for line in eachline(in_fn)
                if '@' ∉ line
                    println(out_fh, line)
                end
            end
        end

    end
    ))

    # Read the original file
    all_tweets = Set(eachline(datadep"TrumpTweets/realDonaldTrump_poll_tweets.csv"))
    # Read the file that we are generating
    nonmentions_tweets = Set(eachline(datadep"TrumpTweets/nonmentions_realDonaldTrump_poll_tweets.csv"))
    # Use them both
    mentions_tweets = setdiff(all_tweets, nonmentions_tweets)
    @test length(mentions_tweets) > 0
    @test all(Ref('@') .∈ collect(mentions_tweets))
end

@testset "MNIST" begin

    register(DataDep(
        "MNIST train",
        """
        Dataset: THE MNIST DATABASE of handwritten digits, (training subset)
        Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
        Website: http://yann.lecun.com/exdb/mnist/
        [LeCun et al., 1998a]
            Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.
            "Gradient-based learning applied to document recognition."
            Proceedings of the IEEE, 86(11):2278-2324, November 1998
        The files are available for download at the offical
        website linked above. We can download these files for you
        if you wish, but that doesn't free you from the burden of
        using the data responsibly and respect copyright. The
        authors of MNIST aren't really explicit about any terms
        of use, so please read the website to make sure you want
        to download the dataset.
        """,
        "http://yann.lecun.com/exdb/mnist/".*["train-images-idx3-ubyte.gz", "train-labels-idx1-ubyte.gz"];
        # Not providing a checksum at all so can check it gives output
        # TODO : automate this test with new 0.7 stuff
    ))

    register(DataDep(
        "MNIST",
        """
        Dataset: THE MNIST DATABASE of handwritten digits
        Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
        Website: http://yann.lecun.com/exdb/mnist/
        [LeCun et al., 1998a]
            Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.
            "Gradient-based learning applied to document recognition."
            Proceedings of the IEEE, 86(11):2278-2324, November 1998
        The files are available for download at the offical
        website linked above. We can download these files for you
        if you wish, but that doesn't free you from the burden of
        using the data responsibly and respect copyright. The
        authors of MNIST aren't really explicit about any terms
        of use, so please read the website to make sure you want
        to download the dataset.
        """,
        "http://yann.lecun.com/exdb/mnist/".*["train-images-idx3-ubyte.gz", "train-labels-idx1-ubyte.gz", "t10k-images-idx3-ubyte.gz", "t10k-labels-idx1-ubyte.gz"],
        "0bb1d5775d852fc5bb32c76ca15a7eb4e9a3b1514a2493f7edfcf49b639d7975"
    ))
read(datadep"MNIST"*"/train-labels-idx1-ubyte.gz")
    @test read(datadep"MNIST"*"/train-labels-idx1-ubyte.gz") == read(datadep"MNIST train"*"/train-labels-idx1-ubyte.gz")
end

@testset "UCI Banking" begin
    register(DataDep(
        "UCI Banking",
        """
        Dataset: Bank Marketing Data Set
        Authors: S. Moro, P. Cortez and P. Rita.
        Website: https://archive.ics.uci.edu/ml/datasets/bank+marketing

        This dataset is public available for research. The details are described in [Moro et al., 2014].
        Please include this citation if you plan to use this database:
        [Moro et al., 2014] S. Moro, P. Cortez and P. Rita. A Data-Driven Approach to Predict the Success of Bank Telemarketing. Decision Support Systems, Elsevier, 62:22-31, June 2014
        """,
        [
        "https://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank.zip",
        "https://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank-additional.zip"
        ],
        [(SHA.sha1, "785118991cd7d7ee7d8bf75ea58b6fae969ac185"),
         (SHA.sha3_224, "01b53f5b69d0b169070219b4391c623d84ab17d4cea8c8895cbf951d")];

         post_fetch_method = unpack
    ))

    data, header = readdlm(datadep"UCI Banking/bank.csv", ';', header=true)
    @test size(header) == (1,17)
    @test size(data) == (4521,17)

end

@testset "UCI Adult, Hierarchical checksums" begin
    # This is an example of using hierachacy in the remote URLs,
    # and similar (partially matching up to depth) hierachacy in the checksums
    # for processing some groups of elements differently to others.
    # Doing this with checksums is not particularly useful
    # But the same thing applies to `fetch_method` and `post_fetch_method`.
    # So for example the
    register(DataDep(
        "UCI Adult",
        """
        Dataset: Adult Data Set UCI ML Repository
        Website: https://archive.ics.uci.edu/ml/datasets/Adult
        Abstract : Predict whether income exceeds \$50K/yr based on census data.  Also known as "Census Income" dataset.

        If you make use of this data it is requested that you cite:
        - Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science.
        """,
        [
            [
                "https://archive.ics.uci.edu/ml/datasets/../machine-learning-databases/adult/adult.data",
                "https://archive.ics.uci.edu/ml/datasets/../machine-learning-databases/adult/adult.test"
            ],

            [
                "https://archive.ics.uci.edu/ml/datasets/../machine-learning-databases/adult/Index",
                [
                    "https://archive.ics.uci.edu/ml/datasets/../machine-learning-databases/adult/adult.names"
                    "https://archive.ics.uci.edu/ml/datasets/../machine-learning-databases/adult/old.adult.names"
                ]
             ]
        ],
        [
            "f9a9220df6bc5d9848bf450fd9ad45b9496503551af387d4a1bbe38ce1f8fc38", #adult.data ⊻ adult.test
            [
             "c53c35ce8a0eb10c12dd4b73830f3c94ae212bb388389d3763bce63e8d6bc684", #Index
             "818481d320861c4b623626ff6fab3426ad93dae4434b7f54ca5a0f357169c362" # adult.names ⊻ old.adult.names
            ]
        ]
    ))

    @test length(collect(eachline(datadep"UCI Adult/adult.names"))) == 110

end

@testset "gzipped source code" begin
    register(DataDep("DataDeps Source v0.5.0",
        """
        This is the source code of DataDeps.jl v0.5.0
        This test checked we can unpack gzipped tarballs.
        """,
        "https://github.com/oxinabox/DataDeps.jl/archive/v0.5.0.tar.gz",
        "cd1fc3e58b4272ec559d1c5bcda5e4f0339647dab709baa2d507a73f3a89168d";
        post_fetch_method=DataDeps.unpack
    ));
    @test length(readdir(datadep"DataDeps Source v0.5.0")) == 2
end


text:
# Type handling

"Returns the signature of the matching method, specializing types as necessary."
function signature(fn_type::Type, arg_types::Type...)
    meta = IRTools.meta(Tuple{fn_type, arg_types...})
    if isnothing(meta) return nothing end
    sig, typevars = unparameterize(meta.method.sig)
    sig_types = collect(sig.parameters[2:end])
    for (i, a) in enumerate(arg_types)
        if sig_types[i] isa TypeVar
            sig_types[i] = a
        else
            sig_types[i] = reparameterize(sig_types[i], typevars)
        end
    end
    return sig_types
end

"Unwraps UnionAll types, returning the innermost body with a list of TypeVars."
function unparameterize(@nospecialize(T::Type))
    vars = TypeVar[]
    while T isa UnionAll
        pushfirst!(vars, T.var)
        T = T.body
    end
    return T, vars
end

"Rewraps parametric types with the provided TypeVars."
function reparameterize(@nospecialize(T::Type), vars::Vector{TypeVar})
    if T isa UnionAll
        T, origvars = unparameterize(T)
        vars = unique!([origvars; vars])
    end
    parameters = Set(T.parameters)
    for v in vars
        if (v in parameters) T = UnionAll(v, T) end
    end
    return T
end

"Unwrap `GlobalRef`s, returning the unqualified name."
unwrap(ref::GlobalRef) = ref.name
unwrap(val) = val

# Loop detection.

struct NaturalLoop
    header::Int # Unique entry point to loop
    latch::Int # In-loop block that has backedge to header
    body::Set{Int} # The set of all blocks in the loop
    backedges::Set{Int} # All blocks in CFG with backedges to the header
end

function Base.display(nl::NaturalLoop)
    println("--- Natural loop ---")
    println(" Header : $(nl.header)")
    println(" Blocks : $(nl.body)")
    println(" Latch : $(nl.latch)")
    println(" Backedges : $(nl.backedges)")
    println()
end

Base.display(ls::AbstractVector{NaturalLoop}) = display.(ls)

IRTools.dominators(ir::IR) = dominators(IRTools.CFG(ir))

"Returns all edges in the control flow graph."
function edges(ir::IR; reverse=false)
    neighbors = reverse ? predecessors : successors
    return Dict(s.id => Set(t.id for t in neighbors(s)) for s in blocks(ir))
end

"Prunes unreachable blocks."
function pruneblocks!(ir::IR)
    unreachable = [b.id for b in blocks(ir) if isempty(predecessors(b))][2:end]
    map(i -> deleteblock!(ir, i), reverse(unreachable))
    return ir
end

"Builds a single loop from the IR, header block, and final block (latch)."
function buildloop(ir::IR, header::Int, latch::Int, backedges::Set{Int})
    body = Set{Int}([header])
    stack = [latch]
    while !isempty(stack)
        i = pop!(stack)
        if i in body continue end
        push!(body, i)
        append!(stack, [b.id for b in predecessors(block(ir, i))])
    end
    return NaturalLoop(header, latch, body, backedges)
end

"Detect the set of natural loops for a piece of IR."
function detectloops(ir::IR)
    loops = NaturalLoop[]
    domdict = dominators(ir)
    backedges = Dict{Int,Set{Int}}()
    # Find all backedges
    for (id, doms) in domdict, header in successors(block(ir, id))
        if !(header.id in doms) continue end
        push!(get!(backedges, header.id, Set{Int}()), id)
    end
    # Create loop for each backedge
    loops = [buildloop(ir, header, latch, edges)
             for (header, edges) in backedges for latch in edges]
    # Sort outermost to innermost, then by header index
    sort!(loops, lt=(l1, l2) -> l2.body in l1.body || l1.header < l2.header)
    return loops
end

"Detect the set of natural loops for a piece of IR, pruning unreachable blocks."
detectloops!(ir::IR) = detectloops(pruneblocks!(IRTools.explicitbranch!(ir)))

"Insert preheader for a block in a piece of IR."
function preheader!(ir::IR, i::Int, backedges::AbstractVector{Int}=Int[])
    # Insert preheader, add arguments
    preheader = block!(ir, i)
    header = block(ir, i+1)
    args = [argument!(preheader, insert = false) for a in arguments(header)]
    # Rewire branches from entering blocks to preheader
    backedges = map(e -> e >= i ? e+1 : e, backedges)
    entrances = [p for p in predecessors(header) if !(p.id in backedges)]
    for e in entrances, (bi, br) in enumerate(branches(e))
        if br.block != header.id continue end
        branches(e)[bi] = Branch(br, block=preheader.id)
    end
    # Branch from preheader to header
    branch!(preheader, header, args...)
    return preheader
end

"Insert preheader for a loop."
preheader!(ir::IR, loop::NaturalLoop) =
    preheader!(ir, loop.header, collect(loop.backedges))

"Insert preheaders for every loop in a piece of IR."
function preheaders!(ir::IR, loops::AbstractVector{NaturalLoop})
    preheaders, inserts = Block[], Int[]
    for (n, loop) in enumerate(sort!(loops, by=l->l.header, rev=true))
        backedges = map(e -> e + sum(inserts .<= e), collect(loop.backedges))
        preheader!(ir, loop.header, backedges)
        pushfirst!(preheaders, block(ir, loop.header + length(loops) - n))
        push!(inserts, loop.header) # Track which indices have been inserted
    end
    return preheaders
end

"Add loop counters to each natural loop in a piece of IR."
function loopcounts!(ir::IR)
    # Prune unreachable blocks, detect loops, and add preheaders
    preheaders!(ir, detectloops!(ir))
    # Iterate over loops in topological order
    prev_header = -1
    loops, countvars = NaturalLoop[], Variable[]
    for loop in detectloops(ir)
        # Skip nested loops with the same header
        if loop.header == prev_header continue end
        # Create loop count variables, and pass them along blocks
        preheader, header = block(ir, loop.header-1), block(ir, loop.header)
        push!(arguments(branches(preheader)[1]), 0) # Initialize loop count
        count = argument!(header, insert=false) # Pass to header
        count = pushfirst!(header, xcall(Base, :+, count, 1)) # Increment count
        for b in (block(ir, i) for i in loop.backedges)
            for (bi, br) in enumerate(branches(b)) # Pass count along backedges
                if br.block != loop.header continue end
                push!(arguments(branches(b)[bi]), count)
            end
        end
        push!(loops, loop); push!(countvars, count)
        prev_header = loop.header
    end
    return loops, countvars
end


text:
using Oceananigans.Operators

struct EnergyConservingScheme end
struct EnstrophyConservingScheme end

struct VectorInvariant{S}
    scheme :: S
end

VectorInvariant(; scheme::S = EnstrophyConservingScheme()) where S = VectorInvariant{S}(scheme)

const VectorInvariantSchemes = Union{VectorInvariant, WENOVectorInvariant}

const VectorInvariantEnergyConserving = VectorInvariant{<:EnergyConservingScheme}
const VectorInvariantEnstrophyConserving = VectorInvariant{<:EnstrophyConservingScheme}

######
###### Horizontally-vector-invariant formulation of momentum scheme
######
###### Follows https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#vector-invariant-momentum-equations
######

@inline U_dot_∇u(i, j, k, grid, scheme::VectorInvariantSchemes, U) = (
    + vertical_vorticity_U(i, j, k, grid, scheme, U.u, U.v)  # Vertical relative vorticity term
    + vertical_advection_U(i, j, k, grid, scheme, U.u, U.w)  # Horizontal vorticity / vertical advection term
    + bernoulli_head_U(i, j, k, grid, scheme, U.u, U.v))     # Bernoulli head term

@inline U_dot_∇v(i, j, k, grid, scheme::VectorInvariantSchemes, U) = (
    + vertical_vorticity_V(i, j, k, grid, scheme, U.u, U.v)  # Vertical relative vorticity term
    + vertical_advection_V(i, j, k, grid, scheme, U.v, U.w)  # Horizontal vorticity / vertical advection term
    + bernoulli_head_V(i, j, k, grid, scheme, U.u, U.v))     # Bernoulli head term

####
#### Bernoulli head terms
####

@inline bernoulli_head_U(i, j, k, grid, scheme::VectorInvariantSchemes, u, v) = ∂xᶠᶜᶜ(i, j, k, grid, Khᶜᶜᶜ, scheme, u, v)
@inline bernoulli_head_V(i, j, k, grid, scheme::VectorInvariantSchemes, u, v) = ∂yᶜᶠᶜ(i, j, k, grid, Khᶜᶜᶜ, scheme, u, v)

@inline ϕ²(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2
@inline Khᶜᶜᶜ(i, j, k, grid, ::VectorInvariantSchemes, u, v) = (ℑxᶜᵃᵃ(i, j, k, grid, ϕ², u) + ℑyᵃᶜᵃ(i, j, k, grid, ϕ², v)) / 2

####
#### Horizontal advection terms
####

@inline ζ_ℑx_vᶠᶠᵃ(i, j, k, grid, u, v) = ζ₃ᶠᶠᶜ(i, j, k, grid, u, v) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
@inline ζ_ℑy_uᶠᶠᵃ(i, j, k, grid, u, v) = ζ₃ᶠᶠᶜ(i, j, k, grid, u, v) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)

@inline vertical_vorticity_U(i, j, k, grid, ::VectorInvariantEnergyConserving, u, v) = - ℑyᵃᶜᵃ(i, j, k, grid, ζ_ℑx_vᶠᶠᵃ, u, v) / Δxᶠᶜᶜ(i, j, k, grid)
@inline vertical_vorticity_V(i, j, k, grid, ::VectorInvariantEnergyConserving, u, v) = + ℑxᶜᵃᵃ(i, j, k, grid, ζ_ℑy_uᶠᶠᵃ, u, v) / Δyᶜᶠᶜ(i, j, k, grid)

@inline vertical_vorticity_U(i, j, k, grid, ::VectorInvariantEnstrophyConserving, u, v) = - ℑyᵃᶜᵃ(i, j, k, grid, ζ₃ᶠᶠᶜ, u, v) * ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_qᶜᶠᶜ, v) / Δxᶠᶜᶜ(i, j, k, grid)
@inline vertical_vorticity_V(i, j, k, grid, ::VectorInvariantEnstrophyConserving, u, v) = + ℑxᶜᵃᵃ(i, j, k, grid, ζ₃ᶠᶠᶜ, u, v) * ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_qᶠᶜᶜ, u) / Δyᶜᶠᶜ(i, j, k, grid)

@inline function vertical_vorticity_U(i, j, k, grid, scheme::WENOVectorInvariant{FT, XT, YT, ZT, XS, YS, ZS, VI}, u, v) where {FT, XT, YT, ZT, XS, YS, ZS, VI}
    v̂  =  ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_qᶜᶠᶜ, v) / Δxᶠᶜᶜ(i, j, k, grid)
    ζᴸ =  _left_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme, ζ₃ᶠᶠᶜ, VI, u, v)
    ζᴿ = _right_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme, ζ₃ᶠᶠᶜ, VI, u, v)
    return - upwind_biased_product(v̂, ζᴸ, ζᴿ)
end

@inline function vertical_vorticity_V(i, j, k, grid, scheme::WENOVectorInvariant{FT, XT, YT, ZT, XS, YS, ZS, VI}, u, v) where {FT, XT, YT, ZT, XS, YS, ZS, VI}
    û  =  ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_qᶠᶜᶜ, u) / Δyᶜᶠᶜ(i, j, k, grid)
    ζᴸ =  _left_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, ζ₃ᶠᶠᶜ, VI, u, v)
    ζᴿ = _right_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, ζ₃ᶠᶠᶜ, VI, u, v)
    return + upwind_biased_product(û, ζᴸ, ζᴿ)
end

####
#### Vertical advection terms
####

@inbounds ζ₂wᶠᶜᶠ(i, j, k, grid, u, w) = ℑxᶠᵃᵃ(i, j, k, grid, Az_qᶜᶜᶠ, w) * ∂zᶠᶜᶠ(i, j, k, grid, u)
@inbounds ζ₁wᶜᶠᶠ(i, j, k, grid, v, w) = ℑyᵃᶠᵃ(i, j, k, grid, Az_qᶜᶜᶠ, w) * ∂zᶜᶠᶠ(i, j, k, grid, v)

@inline vertical_advection_U(i, j, k, grid, ::VectorInvariant, u, w) =  ℑzᵃᵃᶜ(i, j, k, grid, ζ₂wᶠᶜᶠ, u, w) / Azᶠᶜᶜ(i, j, k, grid)
@inline vertical_advection_V(i, j, k, grid, ::VectorInvariant, v, w) =  ℑzᵃᵃᶜ(i, j, k, grid, ζ₁wᶜᶠᶠ, v, w) / Azᶜᶠᶜ(i, j, k, grid)

@inline function vertical_advection_U(i, j, k, grid, scheme::WENOVectorInvariant{FT, XT, YT, ZT, XS, YS, ZS, VI}, u, w) where {FT, XT, YT, ZT, XS, YS, ZS, VI}
    ŵ  =  ℑxᶠᵃᵃ(i, j, k, grid, ℑzᵃᵃᶜ, Δx_qᶜᶜᶠ, w) / Δxᶠᶜᶜ(i, j, k, grid)
    ζᴸ =  _left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, ∂zᶠᶜᶠ, VI, u)
    ζᴿ = _right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, ∂zᶠᶜᶠ, VI, u)
    return upwind_biased_product(ŵ, ζᴸ, ζᴿ)
end

@inline function vertical_advection_V(i, j, k, grid, scheme::WENOVectorInvariant{FT, XT, YT, ZT, XS, YS, ZS, VI}, v, w) where {FT, XT, YT, ZT, XS, YS, ZS, VI}
    ŵ  =  ℑyᵃᶠᵃ(i, j, k, grid, ℑzᵃᵃᶜ, Δy_qᶜᶜᶠ, w) / Δyᶜᶠᶜ(i, j, k, grid)
    ζᴸ =  _left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, ∂zᶜᶠᶠ, VI, v)
    ζᴿ = _right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, ∂zᶜᶠᶠ, VI, v)
    return upwind_biased_product(ŵ, ζᴸ, ζᴿ)
end

######
###### Conservative formulation of momentum advection
######

@inline U_dot_∇u(i, j, k, grid, scheme::AbstractAdvectionScheme, U) = div_𝐯u(i, j, k, grid, scheme, U, U.u)
@inline U_dot_∇v(i, j, k, grid, scheme::AbstractAdvectionScheme, U) = div_𝐯v(i, j, k, grid, scheme, U, U.v)

######
###### No advection
######

@inline U_dot_∇u(i, j, k, grid::AbstractGrid{FT}, scheme::Nothing, U) where FT = zero(FT)
@inline U_dot_∇v(i, j, k, grid::AbstractGrid{FT}, scheme::Nothing, U) where FT = zero(FT)


text:
@inline function U!(k, U)
  @inbounds for r = 1:k
    U[1,r] = -r
    for j = 2:k
      U[j,r] = U[j-1,r] * ((j-1) - r)/j
    end
  end
  nothing
end

function R!(k, ρ, cache)
  @unpack R = cache
  @inbounds for r = 1:k
    R[1,r] = -r * ρ
    for j = 2:k
      R[j,r] = R[j-1,r] * ((j-1) - r * ρ)/j
    end
  end
  nothing
end

# This functions takes help of D2 array to create backward differences array D
# Ith row of D2 keeps Ith order backward differences (∇ⁱyₙ)

function backward_diff!(cache::OrdinaryDiffEqMutableCache, D, D2, k, flag=true)
  flag && copyto!(D[1], D2[1,1])
  for i = 2:k
    for j = 1:(k-i+1)
      @.. D2[i,j] = D2[i-1,j] - D2[i-1,j+1]
    end
    flag && copyto!(D[i], D2[i,1])
  end
end

function backward_diff!(cache::OrdinaryDiffEqConstantCache, D, D2, k, flag=true)
  flag && (D[1] = D2[1,1])
  for i = 2:k
    for j = 1:(k-i+1)
      D2[i,j] = D2[i-1,j] - D2[i-1,j+1]
    end
    flag && (D[i] = D2[i,1])
  end
end

# this function updates backward difference array D when stepsize gets change
# Formula -> D = D * (R * U)
# and it is taken from the paper -
# Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s paper
# E. Alberdi Celaya, J. J. Anza Aguirrezabala, and P. Chatzipantelidis
function reinterpolate_history!(cache::OrdinaryDiffEqMutableCache, D, R, k)
  @unpack tmp = cache.nlsolver
  fill!(tmp,zero(eltype(D[1])))
  for j = 1:k
    for k = 1:k
      @. tmp += D[k] * R[k,j]
    end
    D[j] .= tmp
    fill!(tmp, zero(eltype(tmp)))
  end
end

function reinterpolate_history!(cache::OrdinaryDiffEqConstantCache, D, R, k)
  tmp = zero(D[1])
  for j = 1:k
    for k = 1:k
      tmp += D[k] * R[k,j]
    end
    D[j] = tmp
  end
end

function calc_R(ρ, k, ::Val{N}) where {N}
  R = zero(MMatrix{N,N,typeof(ρ)})
  @inbounds for r = 1:k
    R[1,r] = -r * ρ
    for j = 2:k
      R[j,r] = R[j-1,r] * ((j-1) - r * ρ)/j
    end
  end
  SArray(R)
end

function update_D!(D, dd, k)
  dd = _vec(dd)
  @views @.. D[:,k+2] = dd - D[:,k+1]
  @views @.. D[:,k+1] = dd
  for i in k:-1:1
    @views @.. D[:,i] = D[:,i] + D[:,i+1]
  end
  return nothing
end

const γₖ = @SVector[sum(1//j for j in 1:k) for k in 1:6]

error_constant(integrator, order) = error_constant(integrator, integrator.alg, order)
function error_constant(integrator, alg::QNDF, k)
  @unpack γₖ = integrator.cache
  κ = alg.kappa[k]
  κ*γₖ[k] + inv(k+1)
end

function compute_weights!(ts,k,weights)
  for i = 1:k+1
      weights[i] = one(eltype(weights))
      for j = 1:i-1
          weights[i] *= ts[i] - ts[j]
      end
      for j = i+1:k+1
          weights[i] *= ts[i] - ts[j]
      end
      weights[i] = 1/weights[i]
  end
end

function calc_Lagrange_interp(k,weights,t,ts,u_history,u::Number)
  if t in ts
    i = searchsortedfirst(ts,t,rev=true)
    return u_history[i]
  else
    for i in 1:k+1
      u += weights[i]/(t-ts[i])*u_history[i]
    end
    for i in 1:k+1
      u *= t-ts[i]
    end
  end
  u
end

function calc_Lagrange_interp(k,weights,t,ts,u_history,u)
  vc = _vec(u)
  if t in ts
    i = searchsortedfirst(ts,t,rev=true)
    @.. @views vc = u_history[:,i]
    return u
  else
    for i in 1:k+1
      @.. @views vc += weights[i]/(t-ts[i])*u_history[:,i]
    end
    for i in 1:k+1
      @.. u *= t-ts[i]
    end
  end
  u
end

function calc_Lagrange_interp!(k,weights,t,ts,u_history,u)
  vc = _vec(u)
  if t in ts
    i = searchsortedfirst(ts,t,rev=true)
    @.. @views vc = u_history[:,i]
  else
    for i in 1:k+1
      @.. @views vc += weights[i]/(t-ts[i])*u_history[:,i]
    end
    for i in 1:k+1
      @.. u *= t-ts[i]
    end
  end
end

#This code refers to https://epubs.siam.org/doi/abs/10.1137/S0036144596322507

function calc_finite_difference_weights(ts,t,order,::Val{N}) where {N}
  max_order = N
  c = zero(MMatrix{max_order+1,max_order+1,eltype(ts)})
  c1 = one(t)
  c4 = ts[1] - t
  c[1,1] = one(t)
  for i in 2:order+1
    c2 = one(t)
    c5 = c4
    c4 = ts[i] - t
    @inbounds for j = 1:i-1
      c3 = ts[i] - ts[j]
      c2 *= c3
      if j == i-1
        for k in i:-1:2
          c[i,k] = c1*((k-1)*c[i-1,k-1]-c5*c[i-1,k])/c2
        end
        c[i,1] = zero(t)
      end
      for k in i:-1:2
        c[j,k] = (c4*c[j,k] - (k-1)*c[j,k-1])/c3
      end
      c[j,1] = zero(t)
    end
    c1 = c2
  end
  return SArray(c)
end


text:
using DFTK
using Test
using LinearAlgebra: norm

include("testcases.jl")

function test_chi0(;symmetry=false, use_symmetry=false, temperature=0,
                   spin_polarization=:none,
                   kgrid=[3, 1, 1], fft_size=[10, 1, 10], Ecut=3)
    testcase = silicon
    n_bands  = 12
    tol      = 1e-14
    ε        = 1e-8
    testtol  = 2e-6

    collinear = spin_polarization == :collinear
    label = [
        symmetry        ? "   symm" : "no symm",
        use_symmetry    ? "   use" : "no use",
        temperature > 0 ? "temp" : "  0K",
        collinear       ? "coll" : "none",
    ]
    @testset "Computing χ0 ($(join(label, ", ")))" begin
        spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
        magnetic_moments = collinear ? [spec => 2rand(2)] : []
        model_kwargs = (temperature=temperature, symmetries=symmetry,
                        magnetic_moments=magnetic_moments,
                        spin_polarization=spin_polarization)
        basis_kwargs = (kgrid=kgrid, fft_size=fft_size, use_symmetry=use_symmetry)
        model = model_LDA(testcase.lattice, [spec => testcase.positions]; model_kwargs...)
        basis = PlaneWaveBasis(model, Ecut; basis_kwargs...)

        ρ0     = guess_density(basis)
        ρspin0 = guess_spin_density(basis, magnetic_moments)
        energies, ham0 = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0, ρspin=ρspin0)
        res = DFTK.next_density(ham0, tol=tol, eigensolver=diag_full, n_bands=n_bands)
        ρ1  = res.ρout
        ρspin1 = res.ρ_spin_out

        # Now we make the same model, but add an artificial external potential ε * dV
        n_spin = model.n_spin_components
        dV = [from_real(basis, randn(eltype(basis), basis.fft_size)) for _ in 1:n_spin]
        dV_sym = DFTK.symmetrize.(dV)
        if symmetry
            dV = dV_sym
        else
            for σ in 1:n_spin
                @test dV_sym[σ].real ≈ dV[σ].real
            end
        end

        εdV = cat((ε .* dV[σ].real for σ in 1:n_spin)...; dims=4)
        term_builder = basis -> DFTK.TermExternal(basis, εdV)
        model = model_LDA(testcase.lattice, [spec => testcase.positions];
                          model_kwargs..., extra_terms=[term_builder])
        basis = PlaneWaveBasis(model, Ecut; basis_kwargs...)
        energies, ham = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0, ρspin=ρspin0)
        res = DFTK.next_density(ham, tol=tol, eigensolver=diag_full, n_bands=n_bands)
        ρ2     = res.ρout
        ρspin2 = res.ρ_spin_out
        diff_findiff = [(ρ2 - ρ1) / ε]
        !isnothing(ρspin1) && push!(diff_findiff, (ρspin2 - ρspin1) / ε)

        EVs = [eigen(Hermitian(Array(Hk))) for Hk in ham0.blocks]
        Es = [EV.values[1:n_bands] for EV in EVs]
        Vs = [EV.vectors[:, 1:n_bands] for EV in EVs]
        occ, εF = DFTK.compute_occupation(basis, Es)

        # Test apply_χ0 and compare against finite differences
        diff_applied_χ0 = apply_χ0(ham0, Vs, εF, Es, dV...; droptol=0)
        for σ in 1:n_spin
            @test norm(diff_findiff[σ].real - diff_applied_χ0[σ].real) < testtol
        end

        if !symmetry
            dV_vec = vcat(vec.(getproperty.(dV, :real))...)

            # Test compute_χ0 against finite differences
            χ0 = compute_χ0(ham0)
            diff_computed_χ0 = real(reshape(χ0 * dV_vec, basis.fft_size..., n_spin))
            for σ in 1:n_spin
                @test norm(diff_findiff[σ].real - diff_computed_χ0[:, :, :, σ]) < testtol
            end

            # Test that apply_χ0 is self-adjoint
            dV1 = [from_real(basis, randn(eltype(basis), basis.fft_size)) for _ in 1:n_spin]
            dV2 = [from_real(basis, randn(eltype(basis), basis.fft_size)) for _ in 1:n_spin]
            χ0dV1 = apply_χ0(ham0, Vs, εF, Es, dV1...)
            χ0dV2 = apply_χ0(ham0, Vs, εF, Es, dV2...)
            if n_spin == 1
                @test abs(  dot(dV1[1].real, χ0dV2[1].real)
                          - dot(dV2[1].real, χ0dV1[1].real)) < testtol
            else
                # Form linear combination to get density changes in α and β
                χ0dV1αβ = [(χ0dV1[1].real + χ0dV1[2].real) / 2,  # (ρtot + ρspin) / 2
                           (χ0dV1[1].real - χ0dV1[2].real) / 2]  # (ρtot - ρspin) / 2
                χ0dV2αβ = [(χ0dV2[1].real + χ0dV2[2].real) / 2,
                           (χ0dV2[1].real - χ0dV2[2].real) / 2]
                V1χ0V2 = dot(χ0dV1αβ[1], dV2[1].real) + dot(χ0dV1αβ[2], dV2[2].real)
                V2χ0V1 = dot(χ0dV2αβ[1], dV1[1].real) + dot(χ0dV2αβ[2], dV1[2].real)
                @test abs(V1χ0V2 - V2χ0V1) < testtol
            end

            # Test the diagonal_only option
            χ0_diag = compute_χ0(ham0; droptol=Inf)
            diff_diag_1 = real(reshape(χ0_diag * dV_vec, basis.fft_size..., n_spin))
            diff_diag_2 = apply_χ0(ham0, Vs, εF, Es, dV...; droptol=Inf,
                                   sternheimer_contribution=false)
            for σ in 1:n_spin
                @test norm(diff_diag_1[:, :, :, σ] - diff_diag_2[σ].real) < testtol
            end
        end
    end
end

for temperature in (0, 0.03), spin_polarization in (:none, :collinear)
    for use_symmetry in (false, true), symmetry in (false, true)
        test_chi0(;symmetry=symmetry, use_symmetry=use_symmetry,
                  temperature=temperature, spin_polarization=spin_polarization)
    end
end


text:
using Statistics
using Oceananigans.Buoyancy: g_Earth
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface, FreeSurface,
      compute_volume_scaled_divergence!, add_previous_free_surface_contribution,
      compute_vertically_integrated_volume_flux!

function run_implicit_free_surface_solver_tests(arch, grid)
    ### Nx = 360
    ### Ny = 360
    # A spherical domain with bounded toplogy
    ### grid = RegularLatitudeLongitudeGrid(size = (Nx, Ny, 1),
    ###                                 longitude = (-30, 30),
    ###                                 latitude = (15, 75),
    ###                                 z = (-4000, 0))

    # A boring grid
    ## RegularRectilinearGrid(FT, size=(3, 1, 4), extent=(3, 1, 4))
    ### grid = RegularRectilinearGrid(size = (128, 1, 1),
    ###                         x = (0, 1000kilometers), y = (0, 1), z = (-400, 0),
    ###                         topology = (Bounded, Periodic, Bounded))
    Nx = grid.Nx
    Ny = grid.Ny
    Δt = 900

    free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth)

    # Create fields with default bounded and zeroflux boundaries.
    velocities = VelocityFields(arch, grid)
    η          = CenterField(arch, grid, TracerBoundaryConditions(grid))
    @. η.data  = 0

    # Initialize the solver
    free_surface = FreeSurface(free_surface, velocities, arch, grid)

    # Create a divergent velocity
    u, v, w = velocities
    @. u.data = 0
    @. v.data = 0
    @. w.data = 0
    imid = Int(floor(grid.Nx / 2)) + 1
    jmid = Int(floor(grid.Ny / 2)) + 1
    CUDA.@allowscalar u.data[imid, jmid, 1] = 1

    fill_halo_regions!(u.data, u.boundary_conditions, arch, grid, nothing, nothing)
    fill_halo_regions!(v.data, v.boundary_conditions, arch, grid, nothing, nothing)
    fill_halo_regions!(w.data, w.boundary_conditions, arch, grid, nothing, nothing)

    # Create a fake model
    model = (architecture=arch,grid=grid,free_surface=free_surface, Δt=Δt, velocities=velocities)

    ## We need vertically integrated U,V
    event = compute_vertically_integrated_volume_flux!(free_surface, model)
    wait(device(model.architecture), event)
    u=free_surface.barotropic_volume_flux.u
    v=free_surface.barotropic_volume_flux.v
    fill_halo_regions!(u.data ,u.boundary_conditions, arch, grid, nothing, nothing )
    fill_halo_regions!(v.data ,v.boundary_conditions, arch, grid, nothing, nothing )

    ### We don't need the halo below, its just here for some debugging
    Ax=free_surface.vertically_integrated_lateral_face_areas.Ax
    Ay=free_surface.vertically_integrated_lateral_face_areas.Ay
    fill_halo_regions!(Ax.data ,Ax.boundary_conditions, arch, grid, nothing, nothing )
    fill_halo_regions!(Ay.data ,Ay.boundary_conditions, arch, grid, nothing, nothing )

    # Calculate the vertically ingrated divergence term
    event = compute_volume_scaled_divergence!(free_surface, model)
    wait(device(model.architecture), event)

    # Scale and add in η term
    RHS = free_surface.implicit_step_solver.solver.settings.RHS
    RHS .= RHS/(free_surface.gravitational_acceleration*Δt)
    event = add_previous_free_surface_contribution(free_surface, model, Δt )
    wait(device(model.architecture), event)
    fill_halo_regions!(RHS   , η.boundary_conditions, model.architecture, model.grid)

    x  = free_surface.implicit_step_solver.solver.settings.x
    x .= η.data
    fill_halo_regions!(x ,η.boundary_conditions, model.architecture, model.grid)
    solve_poisson_equation!(free_surface.implicit_step_solver.solver, RHS, x; Δt=Δt, g=free_surface.gravitational_acceleration)
    ## exit()
    fill_halo_regions!(x ,η.boundary_conditions, model.architecture, model.grid)
    free_surface.η.data .= x

    # Amatrix_function!(result, x, arch, grid, bcs; args...)
    result = free_surface.implicit_step_solver.solver.settings.A(x;Δt=Δt,g=free_surface.gravitational_acceleration)

    CUDA.@allowscalar begin
     @test abs(minimum(result[1:Nx, 1:Ny, 1] .- RHS[1:Nx, 1:Ny, 1])) < 1e-11
     @test abs(maximum(result[1:Nx, 1:Ny, 1] .- RHS[1:Nx, 1:Ny, 1])) < 1e-11
     @test std(result[1:Nx, 1:Ny, 1] .- RHS[1:Nx, 1:Ny, 1]) < 1e-13
    end

    return CUDA.@allowscalar result[1:Nx, 1:Ny, 1] ≈ RHS[1:Nx, 1:Ny, 1]
end

@testset "Implicit free surface solver tests" begin
    for arch in archs
        @info "Testing implicit free surface solver [$(typeof(arch))]..."
        # A spherical domain with bounded toplogy
        Nx = 360
        Ny = 360
        grid = RegularLatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                    longitude = (-30, 30),
                                    latitude = (15, 75),
                                    z = (-4000, 0))
        @test run_implicit_free_surface_solver_tests(arch, grid)
        grid = RegularRectilinearGrid(size = (128, 1, 1),
                            x = (0, 1000kilometers), y = (0, 1), z = (-400, 0),
                            topology = (Bounded, Periodic, Bounded))
        @test run_implicit_free_surface_solver_tests(arch, grid)

    end
end


text:
import DataFrames.DataFrame
import GZip
import ProgressMeter

# IEX pcap file intro constants
const IEX_VERSION = 0x01
const IEX_RESERVED = 0x00
const CHANNEL_ID = 0x00000001
const TOPS_PROTOCOL_ID_1_5 = 0x8002
const TOPS_PROTOCOL_ID_1_6 = 0x8003
const DEEP_PROTOCOL_ID_1_0 = 0x8004

# data type for stock tickers
const EightASCIICharBytes = NTuple{8,UInt8}

struct TradeReportMessage
    sale_condition_flags::UInt8
    timestamp::Int64
    symbol::EightASCIICharBytes
    size::UInt32
    price::Int64
    trade_id::Int64
end

# TODO: Parse more message types like QuoteUpdateMessage and TradeBreakMessage
struct QuoteUpdateMessage
    flags::UInt8
    timestamp::Int64
    symbol::EightASCIICharBytes
    bid_size::UInt32
    bid_price::Int64
    ask_price::Int64
    ask_size::UInt32
end

struct TradeBreakMessage
    sale_condition_flags::UInt8
    timestamp::Int64
    symbol::EightASCIICharBytes
    size::UInt32
    price::Int64
    trade_id::Int64
end

"""
Header message format, after the version, reserved, protocol, and channel information
(version, reserved bit, protocol, and channel are fixed for a given feed and version)
"""
struct TPHeaderVaryingHalf
    session_id::UInt32
    payload_length::UInt16
    message_count::UInt16
    stream_offset::Int64
    first_message_sequence_number::Int64
    send_time::Int64
end

struct TradeMessageIter
    pcap_file::AbstractString
    show_progress_bar::Bool
    protocol_magic_bytes::UInt16
end
TradeMessageIter(pcap_file::AbstractString; show_progress_bar::Bool=true, protocol_magic_bytes::UInt16=TOPS_PROTOCOL_ID_1_6) = TradeMessageIter(pcap_file, show_progress_bar, protocol_magic_bytes)

struct TradeMessageIterState
    io::Union{Nothing,IO}
    messages_left::Int
    n_trades::Int
    progress_bar::Union{Nothing,ProgressMeter.ProgressUnknown}
end
TradeMessageIterState() = TradeMessageIterState(nothing, 0, 0, nothing)

function Base.iterate(
    iter::TradeMessageIter,
    state::TradeMessageIterState=TradeMessageIterState()
)::Union{Tuple{TradeReportMessage,TradeMessageIterState},Nothing}
    io = state.io
    if io === nothing
        io = GZip.open(iter.pcap_file, "r")
    end

    n_trades = state.n_trades
    messages_left = state.messages_left
    progress_bar = state.progress_bar
    if iter.show_progress_bar && progress_bar === nothing
        progress_bar = ProgressMeter.ProgressUnknown()
    end

    while !eof(io)
        # seek for next header while no messages left
        while messages_left == 0
            tp_header = seek_header(io, iter.protocol_magic_bytes)
            if tp_header === nothing
                break
            end
            messages_left = tp_header.message_count
        end

        # if seeking got to the end of the file, break out
        if eof(io)
            break
        end

        # read messsage
        message_size = read(io, UInt16)
        message_type = Char(read(io, 1)[1])
        message = nothing
        if message_type == 'T'  # 'T' means TradeReportMessage
            if !(message_size - 1 == sum(sizeof(x) for x in fieldtypes(TradeReportMessage)))
                throw("WTF!")
            end
            message = read_struct(io, TradeReportMessage)
            n_trades += 1
        else
            _ = read(io, message_size - 1)
        end
        messages_left -= 1
        if progress_bar !== nothing
            ProgressMeter.next!(progress_bar, showvalues=[("Trades", n_trades)])
        end
        if message !== nothing
            new_state = TradeMessageIterState(io, messages_left, n_trades, progress_bar)
            return message, new_state
        end
    end
end
Base.IteratorSize(iter::TradeMessageIter) = Base.SizeUnknown()
Base.IteratorEltype(iter::TradeMessageIter) = Base.HasEltype()
Base.eltype(iter::TradeMessageIter) = TradeReportMessage

"""Equivalent of the `read` function but for structs of fixed-size types"""
function read_struct(io::IO, struct_type::Type)
    field_gen = (
        reinterpret(field_type, read(io, sizeof(field_type)))[1]
        for field_type in fieldtypes(struct_type)
    )
    return struct_type(field_gen...)
end

"""
Reads PCAP file contents until header encountered.
Returns the varying portion of the header.
"""
function seek_header(io, protocol_id)::Union{TPHeaderVaryingHalf,Nothing}
    match_bytes = [
        reinterpret(UInt8, [IEX_VERSION]);
        reinterpret(UInt8, [IEX_RESERVED]);
        reinterpret(UInt8, [protocol_id]);
        reinterpret(UInt8, [CHANNEL_ID]);
    ]
    i = 1
    read_bytes = 0
    while !eof(io)
        next_byte = read(io, 1)[1]
        read_bytes += 1
        if next_byte == match_bytes[i]
            i += 1
            # if we match the first part of the header, read and return the second half
            if i == length(match_bytes) + 1
                remaining_header = read_struct(io, TPHeaderVaryingHalf)
                return remaining_header
            end
        else
            i = 1
        end
    end
    # if we hit EOF, return nothing
    return nothing
end

"""Read all TradeReportMessage messages in a PCAP file."""
function read_trade_report_messages(
    gzipped_pcap_filepath::String;
    show_progress_bar::Bool=true,
    protocol_magic_bytes::UInt16=TOPS_PROTOCOL_ID_1_6
)::Vector{TradeReportMessage}
    trade_messages = TradeReportMessage[]
    trade_message_iter = TradeMessageIter(
        gzipped_pcap_filepath;
        show_progress_bar=show_progress_bar,
        protocol_magic_bytes=protocol_magic_bytes
    )
    for message in trade_message_iter
        push!(trade_messages, message)
    end
    return trade_messages
end

"""Convert a list of TradeReportMessage objects into a DataFrame."""
function assemble_trade_report_frame(trade_report_messages::Vector{TradeReportMessage})::DataFrame
    # extract message structs into a columnar format
    cols = Dict([
        (name, Vector{type}(undef, length(trade_report_messages)))
        for (name, type) in
        zip(fieldnames(TradeReportMessage),
            fieldtypes(TradeReportMessage))
    ])
    for (i, msg) in enumerate(trade_report_messages)
        for field in fieldnames(typeof(msg))
            cols[field][i] = getfield(msg, field)
        end
    end

    # convert types
    cols[:symbol] = [strip(transcode(String, [x...])) for x in cols[:symbol]]
    cols[:sale_condition_flags] = Int64.(cols[:sale_condition_flags])
    cols[:size] = Int64.(cols[:size])

    # wrap as DataTable
    df = DataFrame(;cols...)
    return df
end