text:
if PROC=="GPU"

    import Base.sum!
    function sum!(A::CudaArray{Float64},ind::Int,beta::Float64,out::CudaArray{Float64})
        if ind==1
            tmp=CudaArray(Float64,(1,size(A,1))); fill!(tmp,1.0)
            CUBLAS.gemm!('N','N',1.0,tmp,A,beta,out)
            free(tmp)
        elseif ind==2
            tmp=CudaArray(Float64,(size(A,2),1)); fill!(tmp,1.0)
            CUBLAS.gemm!('N','N',1.0,A,tmp,beta,out)
            free(tmp)
        else
            error("sum!(A::CudaArray,ind,beta,out) only defined for ind=1 or 2")
        end
    end

    function sum!(A::CudaArray{Float32},ind::Int,beta::Real,out::CudaArray{Float32})
        if ind==1
            tmp=CudaArray(Float32,(1,size(A,1))); fill!(tmp,1.0)
            CUBLAS.gemm!('N','N',Float32(1.0),tmp,A,Float32(beta),out)
            free(tmp)
        elseif ind==2
            tmp=CudaArray(Float32,(size(A,2),1)); fill!(tmp,1.0)
            CUBLAS.gemm!('N','N',Float32(1.0),A,tmp,Float32(beta),out)
            free(tmp)
        else
            error("sum!(A::CudaArray,ind,beta,out) only defined for ind=1 or 2")
        end
    end

    export sum!

    function ArrayToCudaArray!(nodes,net)
        nds=union(nodes)
        for n in nds
            net.value[n]=CudaArray(net.value[n])
        end
    end

    export ArrayToCudaArray!

    scalval(x::CudaArray)=sqrt(CUBLAS.dot(x,x))
    export scalval

    import CUDArt.to_host
    to_host(A::Array)=A
    export to_host

    import Base.println
    println(A::CudaArray)=println(to_host(A))
    export println
end

function cArray(eltype,sz::Tuple)
    global PROC
     if PROC=="GPU"
        return CudaArray(eltype,sz)
    else
        return zeros(eltype,sz)
    end
end

function cArray(sz::Tuple)
    global PROC
    if PROC=="GPU"
        return CudaArray(Float64,sz)
    else
        return Array(Float64,sz)
    end
end

function cArray(A::Array)
    global PROC
    if PROC=="GPU"
        return CudaArray(A)
    else
        return A
    end
end

function cArray(proc::ASCIIString,A::Array)
    if proc=="GPU"
        return CudaArray(A)
    else
        return A
    end
end

function cArray(proc::ASCIIString,eltype,A::Array)
    if proc=="GPU"
        return CudaArray(eltype,A)
    else
        return map(eltype,A)
    end
end

function cArray(GPUproc::Bool,eltype,A::Array)
    if GPUproc==true
        return CudaArray(eltype,A)
    else
        return map(eltype,A)
    end
end

export cArray

function cArray(gpu::Bool,eltype,sz::Tuple)
    if gpu==true
       return CudaArray(eltype,sz)
    else
        return zeros(eltype,sz)
    end
end

function cArray(gpu,sz::Tuple)
    if gpu==true
        return CudaArray(Float64,sz)
    else
        return Array(Float64,sz)
    end
end
function cArray(gpu,A::Array)
    if gpu==true
        return CudaArray(A)
    else
        return A
    end
end
export cArray

function CudaArray(eltype,A::Array)
    CudaArray(map(eltype,A))
end
export CudaArray

import Base.zeros
@gpu zeros(A::CudaArray)=CudaArray(zeros(size(A)))
zeros(A::Float64)=collect(0.0)
export zeros

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Documenter, AssignOnce

makedocs(
    modules = [AssignOnce],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "AssignOnce.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/goretkin/AssignOnce.jl.git",
)

ext: jl


text:
using Compat
using Compat.Test
using MathOptInterface
const MOI = MathOptInterface

@testset "Functions" begin
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    @testset "Broadcast" begin
        xf = MOI.SingleVariable(x)
        yf = MOI.SingleVariable(y)
        zf = MOI.SingleVariable(z)
        function sum_indices(sv1::MOI.SingleVariable, sv2::MOI.SingleVariable)
            return sv1.variable.value + sv2.variable.value
        end
        @test sum_indices.(xf, [yf, zf]) == [3, 4]
    end
    @testset "Copy" begin
        @testset "VectorOfVariables" begin
            f = MOI.VectorOfVariables([x, y])
            f_copy = copy(f)
            f_copy.variables[2] = z
            @test f.variables[2] == y
        end
    end
end


text:
# Simulation
# Assuming all bodies connected in a tree by joints

# include("RigidBody.jl")
# include("OrientationConversion.jl")
# include("Joint.jl")
# include("Force.jl")
# include("extF.jl")

using LinearAlgebra
using Revise
using DifferentialEquations
using StaticArrays
using BenchmarkTools

global GravityInInertial = [0.0,0.0,-9.806];

struct solSim
    # Structure to store results of Simulation for each rigid body
    r::Matrix{Float64}# where T<:Real
    v::Matrix{Float64}# where T<:Real
    β::Matrix{Float64}# where T<:Real
    βdot::Matrix{Float64}# where T<:Real
end

"""
simulate(tEnd,tSpan,j;g,extFVec)

tEnd: Simulation End time. (Float64) \\
tSpan: Result saved at every specified tSpan. (Float64) \\
j: Tuple of Joints. Length need not be specified. \\
g: Acceleration due to gravity. Default value: [0.0;0.0;-9.806] (Vector{Float64}) \\
extFVec: Vector of extForces acting on each rigid body. (Vector{extForces}) \\

Tolerance of ode solver set. RelTol = 1e-10, AbsTol = 1e-10. \\
Tsit5() solver used. \\
Returns a tuple containing tSim and vector of solutions in solSim form.
"""
function simulate(tEnd::Float64,tInt::Float64,j::Vector{Joint};
                  g::Vector{Float64}=[0.0;0.0;-9.806]);# where T<:Real#,extFVec::Vector{extForces}=Vector{extForces}(undef,1))
    # Ellipsis (...) to facilitate supplying variable number of arguments to the function
    # Initial condition (all bodies connected)
    global GravityInInertial = g
    # global extFList = extFVec
    X0 = Vector{Float64}(j[1].RB1.x)
    for k = 1:length(j)
        append!(X0,j[k].RB2.x)
    end
    # println("X0 type = ", typeof(X0))
    # println("tEnd type = ", typeof(tEnd))
    # Declaring the ODE Problem as per DifferentialEquations convention
    tSpanODE = (0.0, tEnd);

    prob = ODEProblem(mainDynODEDiff!,X0,tSpanODE,j)
    # sol = solve(prob, saveat = tSpan, RK4(), reltol=1e-10, abstol=1e-10)
    # sol = solve(prob, Tsit5(), reltol=1e-10, abstol=1e-10)#, saveat = 0:tSpan:tEnd)
    sol = solve(prob, Tsit5(), saveat=tInt, reltol=1e-10, abstol=1e-10)
    # sol = solve(prob, DP5(), reltol=1e-10, abstol=1e-10)

    solFinal = Vector{solSim}(undef,length(j))
    for i=1:length(j)
        rSol = Matrix{Float64}(undef,(length(sol.t),3))
        vSol = Matrix{Float64}(undef,(length(sol.t),3))
        βSol = Matrix{Float64}(undef,(length(sol.t),4))
        βdotSol = Matrix{Float64}(undef,(length(sol.t),4))
        for t=1:length(sol.t)
            rSol[t,:] = (sol.u[t][14*i+1:14*i+3])
            vSol[t,:] = (sol.u[t][14*i+8:14*i+10])
            βSol[t,:] = (sol.u[t][14*i+4:14*i+7])
            βdotSol[t,:] = (sol.u[t][14*i+11:14*i+14])
        end
    solFinal[i] = solSim(rSol,vSol,βSol,βdotSol)
    end
    return (sol.t, solFinal)
end

function mainDynODEDiff!(dX,X, j::Vector{Joint}, t::Float64)# where T<:Real
    dX[:] = mainDynODEDiff(X,j,t);
end

function mainDynODEDiff(X::Vector{T}, j::Vector{Joint}, t::Float64)::Vector{T} where T<:Real
    # ODE function to be used as per DifferentialEquations convention
    # Create extForcesList storing extForces for each rigid body
    # Create ForceConstraints Array storing constraint forces acting on each rigid body
    global GravityInInertial
    # Update RigidBodies
    updateRigidBody!(j[1].RB1,X[1:14]) # update inertial frame (not moving)
    for k=1:length(j)
        updateRigidBody!(j[k].RB2,X[14*k+1:14*(k+1)])
    end

    extFListCopy = externalFTotal(t,j); # Generates all the external forces explicitly specified, through joint actions or directly.

    U = genU(extFListCopy) # Construct an input matrix from external forces and torques supplied.

    ## Inplace
    dX = zeros(length(X));
    fxdot!(dX, X, U,j,GravityInInertial);

    ## Not inplace
    # dX = fxdot(X,U,j,GravityInInertial);

    # println("t = $t")
    return dX
end

function fxdot!(dX, X::Vector{T}, U::Matrix{S}, j::Vector{Joint}, GravityInInertial::Vector{Float64}) where {T<:Real,S<:Real}
    dX[:] = fxdot(X,U,j,GravityInInertial);
end

function fxdot(X::Vector{T},U::Matrix{S},j::Vector{Joint},GravityInInertial::Vector{Float64}) where {T<:Real,S<:Real}
    nJ = length(j); # Number of joints

    unconstrF, constrF = Constraint(X, U, j, GravityInInertial)

    xdot = Vector{Union{T,S}}(undef,14*(nJ+1));
    xdot[1:14] = zeros(14);

    for k=1:nJ
        xb = X[14*k+1:14*(k+1)]

        unconstrainedF_rb = unconstrF[:,j[k].RB2.bodyID]
        Fc_rb = constrF[:,j[k].RB2.bodyID]
        TotalForce = unconstrainedF_rb[1:3] + Fc_rb[1:3]
        TotalMoment = unconstrainedF_rb[4:7] + Fc_rb[4:7]
        genF = [TotalForce; TotalMoment]; # generalized force acting on RB
        genM = Array(genMatM(X, j[k].RB2)) # generalized mass matrix of RB
        xdot[14*k+1:14*k+7] = X[14*k+8:14*k+14]#xb[8:10]
        xdot[14*k+8:14*k+14] = inv(genM)*genF;
    end

    return xdot
end

function externalFTotal(t::T, j::Vector{Joint})::Vector{extForces} where T<:Real
    # Returns externally applied forces in total, not including gravity. Includes Forces from joints and other explicitly specified external forces.
    extFList = extF(t,j)

    # Generate forces from actuated joints on each body
    # ForceJoints = Matrix{Float64}(undef,6,2*length(j))
    # for k=1:length(j)
    #     ForceJoints[:,2*k-1], ForceJoints[:,2*k] = genJointF(t,j[k])
    # end

    # Add ForceJoints to extFList
    extFListCopy = deepcopy(extFList)
    # for k=1:length(j)
    #     extFListCopy[j[k].RB1.bodyID].Forces =
    #     vcat(extFListCopy[j[k].RB1.bodyID].Forces, reshape(ForceJoints[1:3,2*k-1],(1,3)))
    #     extFListCopy[j[k].RB1.bodyID].Positions =
    #     vcat(extFListCopy[j[k].RB1.bodyID].Positions, reshape(j[k].pos1,(1,3)))
    #     extFListCopy[j[k].RB1.bodyID].Torques =
    #     vcat(extFListCopy[j[k].RB1.bodyID].Torques, reshape(ForceJoints[4:6,2*k-1],(1,3)))
    #
    #     extFListCopy[j[k].RB2.bodyID].Forces =
    #     vcat(extFListCopy[j[k].RB2.bodyID].Forces, reshape(ForceJoints[1:3,2*k],(1,3)))
    #     extFListCopy[j[k].RB2.bodyID].Positions =
    #     vcat(extFListCopy[j[k].RB2.bodyID].Positions, reshape(j[k].pos2,(1,3)))
    #     extFListCopy[j[k].RB2.bodyID].Torques =
    #     vcat(extFListCopy[j[k].RB2.bodyID].Torques, reshape(ForceJoints[4:6,2*k],(1,3)))
    # end
    return extFListCopy
end

function genU(extFList::Vector{extForces})::Matrix{Float64}
    nB = length(extFList); # number of bodies
    u = zeros(6,nB); # Each column in u contains a vector of [F;τ] acting on the corresponding rigid body
    for b=1:nB
        extF_b = extFList[b]
        u[1:3,b] = sum(extF_b.Forces,dims = 1) # Sum up forces acting on body
        for n in 1:size(extF_b.Forces,1) # Sum up moments acting on body
            u[4:6,b] += cross(extF_b.Positions[n,:], extF_b.Forces[n,:])
        end
        u[4:6,b] += sum(extF_b.Torques, dims = 1)[:] # Sum up torques acting on body
    end
    return u
end

function getXU_0(j::Vector{Joint})
    X0 = Vector{Float64}(j[1].RB1.x)
    for k = 1:length(j)
        append!(X0,j[k].RB2.x)
    end

    extFListCopy = externalFTotal(0.0,j); # Generates all the external forces explicitly specified, through joint actions or directly.

    U0 = genU(extFListCopy) # Construct an input matrix from external forces and torques supplied.

    return X0, U0
end

function checkRevJoint(solQuad::solSim, solProp1::solSim, rjCube1::Array{Float64}, rjProp1::Array{Float64})
    ## Function that outputs relevant quantities constrained under a rotational joint
    tLen = size(solQuad.r,1)
    ωCube = Matrix{Float64}(undef,tLen,3);
    ωProp1 = Matrix{Float64}(undef,tLen,3);
    ωProp1InCube = Matrix{Float64}(undef,tLen,3);
    ωRel = Matrix{Float64}(undef,tLen,3);
    jointLoc1 = Matrix{Float64}(undef,length(tSim),3);
    for i=1:tLen
        dcm = quat2dcm(solQuad.β[i,:])
        dcm2 = quat2dcm(solProp1.β[i,:])
        ωCube[i,:] = angVel(solQuad.β[i,:],solQuad.βdot[i,:])
        ωProp1[i,:] = angVel(solProp1.β[i,:],solProp1.βdot[i,:])
        ωProp1InCube[i,:] = quat2dcm(solQuad.β[i,:])*transpose(quat2dcm(solProp1.β[i,:]))*ωProp1[i,:]
        ωRel[i,:] = ωProp1InCube[i,:] - ωCube[i,:]
        jointLoc1[i,:] = solQuad.r[i,:] + transpose(dcm)*rjCube1 - solProp1.r[i,:] - transpose(dcm2)*rjProp1;
        # jointLoc1[i,:] = solProp1.r[i,:] + transpose(dcm2)*rjProp1;
    end
    return ωCube, ωProp1, ωRel, jointLoc1
end

function checkRevJointIn(solPend1::solSim, rj2::Array{Float64})
    ## Function that outputs relevant quantities constrained under a rotational joint when first body is inertial frame
    tLen = size(solPend1.r,1)
    jointLoc1 = Matrix{Float64}(undef,tLen,3)
    ωSol = Matrix{Float64}(undef,tLen,3)
    for i=1:tLen
        tmp_dcm = quat2dcm(solPend1.β[i,:])
        ωSol[i,:] = angVel(solPend1.β[i,:],solPend1.βdot[i,:])
        jointLoc1[i,:] = solPend1.r[i,:] + transpose(tmp_dcm)*rj2
    end
    return ωSol, jointLoc1
end

function joinBodies!(j1::Vector{Joint}, j2::Vector{Joint}, b1_id::Integer, b2_id::Integer, rj1::Vector{T}, rj2::Vector{T}; type::String="Weld", axis::Vector{T}=[0.0;0.0;1.0], k::T=0.0, rL::T=0.0, jointForce::Vector{T} = zeros(T,3), jointTorque::Vector{T} = zeros(T,3))  where T<:Real
    # join the bodies b1 given by b1_id in j1 and b2_id in j2 at rj1 and rj2 respectively
    # Create new joint as a weld joint, unless specified
    len_j1 = length(j1); # old length of j1
    RB1 = body_in_jointVec(j1,b1_id);
    RB2 = body_in_jointVec(j2,b2_id);

    newJ = Joint(RB1, RB2, rj1, rj2, type = type, axis = axis, k = k, rL = rL, jointForce = jointForce, jointTorque = jointTorque);
    push!(j1,newJ)

     # Update bodyIDs of all bodies present in joint tree j2
    for i=1:length(j2)
        j2[i].RB2.bodyID += len_j1;
        push!(j1,j2[i])
    end
    filterInertialJoint!(j1); # This removes the connection of the body in joint tree j2 that was connected to the inertial frame (removes loop inside joint trees)
end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
#=
Implementation of shootf, the Numerical Recipes
algorithm for solving a boundary value problem
by shooting to a fixed point.
=#

# provides c0, G, m_H, k_B, sigma_sb, a_rad, Msun, Rsun, Lsun
include("constants.jl")
include("physics.jl")
include("opacity.jl")
include("derivs.jl")
include("nr.jl")

function load1(m, Pc, Tc, X, Y)
    #=
    Subroutine for shootf.
    Loads the (r, l, P, T) center (m = 0) values for starting the outward
    integration.
    =#
    mu = mu_from_composition(X, Y)
    # m = 0.01 * Msun # needs tuning
    rhoc = rho_ideal(Pc, Tc, mu)
    rc = (3 * m / (4 * pi * rhoc)) ^ (1 / 3)
    lc = dldm(m, 0, 0, 0, Tc, rhoc, X, Y) * m
    # @debug("load1: center values = ", [rc, lc, Pc, Tc])
    return [rc, lc, Pc, Tc]
end # load1

function load2(M, Rs, Ls, X, Y; max_iterations=1000)
    #=
    Subroutine for shootf.
    Loads the (r, l, P, T) surface (m = M) values for starting the inward
    integration.
    =#
    mu = mu_from_composition(X, Y)
    converged = false
    rho_guess = 1e-8 # need a less arbitrary guess?
    Ts = try
        T_surface(Rs, Ls)
    catch
        @warn("load2: Rs = ", Rs)
        @warn("load2: Ls = ", Ls)
    end
    Ps = 0
    count = 1
    while ~converged & (count < 100)
        kappa = opacity(logkappa_spl, rho_guess, Ts)
        Ps = P_surface(M, Rs, kappa)
        rho = rho_ideal(Ps, Ts, mu)
        relerr = (rho - rho_guess) / rho
        if abs(relerr) < 1e-6 # need less arbitrary convergence condition
            converged = true
        else
            rho_guess = rho_guess + relerr * rho / 6
            count += 1
        end # if
    end # while
    if count >= max_iterations
        @warn("load2: Exceeded ", max_iterations, " iterations!")
    end # if
    # @debug("load2: surface values = ", [Rs, Ls, Ps, Ts])
    return [Rs, Ls, Ps, Ts]
end # load2

function init_guess(M, X, Y)
    #=
    Gives a reasonable guess for Pc, Tc, R, and L for a star of mass M
    =#
    yhalf = [Rsun * 10 ^ 0.062544715341505616, Lsun * 10 ^ -0.36630780815983022, 10 ^ 19.700356657315954, 10 ^ 7.2863105891013289]
    y2 = [Rsun * 10 ^ 0.36072469919806871, Lsun * 10 ^ 1.3021335874028637, 10 ^ 17.186879800538229, 10 ^ 7.3557279873577759]
    y5 = [Rsun * 10 ^ 0.59322648100945885, Lsun * 10 ^ 2.8743140385067893, 10 ^ 16.752833342638592, 10 ^ 7.4547058342916337]
    y4 = [Rsun * 10 ^ 0.52962615387538248, Lsun * 10 ^ 2.50275963393719, 10 ^ 16.855407680647915, 10 ^ 7.429616124439673]
    y10 = [Rsun * 10 ^ 0.68191099104456998, Lsun * 10 ^ 3.8447263619458569, 10 ^ 16.517451455410988, 10 ^ 7.5012272142698366]
    if M == Msun
        return [Rsun, Lsun, 2.4e17, 1.6e7]
    elseif M == 0.5 * Msun
        @info("init_guess: from MESA-Web")
        return yhalf
    elseif M == 2 * Msun
        @info("init_guess: from MESA-Web")
        return y2
    elseif (M >= 4 * Msun) & (M <= 5 * Msun)
        @info("init_guess: from MESA-Web")
        dm = M - 4 * Msun
        dydm = (y5 - y4) / Msun
        return y4 + dydm * dm
    elseif M == 10 * Msun
        @info("init_guess: from MESA-Web")
        return y10
    end
    mu = mu_from_composition(X, Y)

    # mass-radius relation
    if M / Msun > 1.3 # then CNO burning dominates
        z1 = 0.8
        z2 = 0.64
    else # then pp chain dominates
        z1 = 0.5
        z2 = 0.1
    end # if
    R = Rsun * (M / Msun) ^ z1 * (mu / 0.61) ^ z2
    # mass-luminosity relation
    L = Lsun * (M / Msun) ^ 3 * (mu / 0.61) ^ 4
    # Pc, Tc guesses, essentially dimensional analysis with fudges
    Pc = 2 * G * M ^ 2 / (pi * R ^ 4)
    Tc = 8 * 0.02 * mu * m_H * G * M / (k_B * R)
    return [R, L, Pc, Tc]
end # init_guess

function profiles(m, M, Rs, Ls, Pc, Tc, X, Y, mf; n=1000)
    #=
    Integrates outward from m to mf and inward from M to mf to construct
    r, l, P, and T profiles from the given initial conditions.
    =#
    D(x, y) = deriv(x, y[1], y[2], y[3], y[4], X, Y)
    # load the center BC guess
    yc = load1(m, Pc, Tc, X, Y)
    # load the surface BC guess
    ys = load2(M, Rs, Ls, X, Y)
    # create the mass grid
    mass_out = collect(linspace(m, mf, n))
    mass_in = collect(linspace(M, mf, 100 * n))
    # integrate out from center to fixed point
    x1, y1 = rk(D, yc, mass_out)
    # integrate in from surface to fixed point
    x2, y2 = rk(D, ys, mass_in)
    return x1, y1, x2, y2
end

function profiles(star::Star, Rs, Ls, Pc, Tc; n=1000)
    return profiles(star.m, star.M, Rs, Ls, Pc, Tc, star.X, star.Y, star.mf)
end

function score(m, M, Rs, Ls, Pc, Tc, X, Y, mf; n=Int(1e3))
    #=
    Function to zero.  Integrates out and in.
    =#
    x1, y1, x2, y2 = profiles(m, M, Rs, Ls, Pc, Tc, X, Y, mf)
    yf1 = y1[end, 1:end]
    yf2 = y2[end, 1:end]
    return (yf1 - yf2) ./ yf1
end # score

function shootf(m, M, X, Y; fixed_point=0.8)
    #=
    Shooting to a fixed point
    Go from initial guesses to a solution to the two-point boundary value
    problem.
    =#
    y0 = init_guess(M, X, Y)
    @info("shootf: Initial guess for M = ", M, ":")
    @info("shootf:     y0 = ", y0)
    mf = fixed_point * M
    f(y) = score(m, M, y[1], y[2], y[3], y[4], X, Y, mf)
    Rs, Ls, Pc, Tc = newton(f, y0)
    return [Rs, Ls, Pc, Tc]
end # shootf

function shootf(star::Star)
    return shootf(star.m, star.M, star.X, star.Y)
end

function shootf!(star::Star)
    bc = shootf(star)
    x1, y1, x2, y2 = profiles(star, bc...)
    star.grid = vcat(x1, reverse(x2))
    star.r = vcat(y1[1:end, 1], reverse(y2[1:end, 1]))
    star.l = vcat(y1[1:end, 2], reverse(y2[1:end, 2]))
    star.P = vcat(y1[1:end, 3], reverse(y2[1:end, 3]))
    star.T = vcat(y1[1:end, 4], reverse(y2[1:end, 4]))
    return x1, y1, x2, y2
end


text:
function test_create()
    create(
        TEST_CONFIG;
        dir=joinpath(TEST_TMP_DIR, "test-playground"),
        name="myproject",
        reqs_file=joinpath(TEST_DIR, "../REQUIRE"),
        julia="julia-bin"
    )
    @test ispath(joinpath(TEST_TMP_DIR, "test-playground"))
    @test isdir(joinpath(TEST_TMP_DIR, "test-playground"))
    @test ispath(joinpath(TEST_CONFIG.dir.store, "myproject"))
    @test islink(joinpath(TEST_CONFIG.dir.store, "myproject"))

    create(
        TEST_CONFIG;
        name="otherproject",
        reqs_file=joinpath(TEST_DIR, "DECLARE"),
        reqs_type=:DECLARE,
        julia="julia-nightly-dir"
    )

    @test ispath(joinpath(TEST_CONFIG.dir.store, "otherproject"))
    @test isdir(joinpath(TEST_CONFIG.dir.store, "otherproject"))

    create(TEST_CONFIG)
    @test ispath(joinpath(TEST_DIR, ".playground"))
    @test isdir(joinpath(TEST_DIR, ".playground"))
end

test_create()


text:
# ------------------------------------------------------------------------------
#                           BLAS LEVEL 2
# ------------------------------------------------------------------------------

if debug println("MKL : BLAS LEVEL 2") end
for fname in (  :mkl_gemv!,)

                if debug println("export : ", fname) end
                @eval export $fname
end

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

## ?gemv
# fname = :mkl_gemv!
# for (fmkl, elty) in ((:dgemv,:Float64),
#                       (:sgemv,:Float32))
#
#     @eval begin
#         function $fname(layout::Cint, trans::Char,
#                         m::Integer, n::Integer,
#                         alpha::Float64, A::Union{Ptr{Float64},DenseArray{Float64}}, lda::Integer, X::Union{Ptr{Float64},DenseArray{Float64}}, incx::Integer,
#                         beta::(Float64), Y::Union{Ptr{Float64},DenseArray{Float64}}, incy::Integer)
#
#             # m,n = size(A,1),size(A,2)
#             # lda = max(1,stride(A,2))
#             #
#             # if trans == 'N' && (length(X) != n || length(Y) != m)
#             #     throw(DimensionMismatch("A has dimensions $(size(A)), X has length $(length(X)) and Y has length $(length(Y))"))
#             # elseif trans == 'C' && (length(X) != m || length(Y) != n)
#             #     throw(DimensionMismatch("A' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
#             # elseif trans == 'T' && (length(X) != m || length(Y) != n)
#             #     throw(DimensionMismatch("A.' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
#             # end
#
#             ccall((:dgemv, libmkl), Void,
#                     (Ptr{Cint}, Ptr{Int},
#                     Ptr{Int}, Ptr{Int},
#                     Ptr{Float64}, Ptr{Float64}, Ptr{Int}, Ptr{Float64}, Ptr{Int},
#                     Ptr{Float64}, Ptr{Float64}, Ptr{Int}),
#                     &layout, &trans,
#                     &m, &n,
#                     &alpha, A, &lda, X, incx,
#                     &beta, Y, &incy)
#             Y
#         end
#     end
# end


text:
abstract type AbstractInterferometryDatum{T} end

abstract type AbstractVisibilityDatum{T} <: AbstractInterferometryDatum{T} end
abstract type AbstractLinearPolDatum{S<:AbstractVisibilityDatum, T} <: AbstractInterferometryDatum{T} end
abstract type AbstractCrossPolDatum{S,T} <: AbstractInterferometryDatum{T} end

abstract type ClosureProducts{T} <: AbstractInterferometryDatum{T} end

abstract type Observation{T} end

using DelimitedFiles
using AstroTime: modified_julian

export uvpositions, stations, getdata, arrayconfig,
       getuv, baselines

"""
    $(TYPEDEF)
This defined the abstract type for an array configuration. Namely, baseline
times, SEFD's, bandwidth, observation frequencies, etc.
"""
abstract type ArrayConfiguration end

struct EHTArrayConfiguration{S,F,T<:AbstractArray} <: ArrayConfiguration
    stations::S
    frequency::F
    bandwidth::F
    uvsamples::T
end

function getuv(ac::ArrayConfiguration)
    return ac.uvsamples.u, ac.uvsamples.v
end

function uvtimefreq(ac::EHTArrayConfiguration)
    u,v = getuv(ac)
    t = ac.uvsamples.time
    ν = ac.frequency
    return u, v, t, fill(ν, length(u))
end

struct ArrayBaselineDatum{T}
    time::T
    u::T
    v::T
    baseline::Tuple{Symbol, Symbol}
    error_real::T
    error_imag::T
end

const ArrayQuadrangleDatum = NTuple{4, ArrayBaselineDatum{T}} where {T}
const ArrayTriangleDatum = NTuple{3, ArrayBaselineDatum{T}} where {T}

uvpositions(D::AbstractVisibilityDatum) = D.u, D.v

Base.@kwdef struct EHTObservation{F,T<:AbstractInterferometryDatum{F},S<:StructArray{T}, N} <: Observation{F}
    data::S
    mjd::N
    ra::F
    dec::F
    bandwidth::F
    frequency::F
    source::Symbol
    timetype::Symbol = :UTC
end

Base.getindex(data::EHTObservation, i::Int) = data.data[i]
Base.length(data::EHTObservation) = length(data.data)

function stations(d::EHTObservation{T,A}) where {T,A<:AbstractInterferometryDatum}
    bl = getdata(d, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    return unique([s1..., s2...])
end

function Base.show(io::IO, d::EHTObservation{F,D}) where {F,D}
    println(io, "EHTObservation{$F,$D, ...}")
    println(io, "  source: ", d.source)
    println(io, "  mjd: ", d.mjd)
    println(io, "  frequency: ", d.frequency)
    println(io, "  bandwidth: ", d.bandwidth)
    println(io, "  stations: ", stations(d))
    println(io, "  nsamples: ", length(d))
end

getdata(obs::Observation, s::Symbol) = getproperty(obs.data, s)
Base.getindex(obs::Observation, i) = getindex(obs.data, i)

Base.@kwdef struct EHTVisibilityDatum{T<:Number} <: AbstractVisibilityDatum{T}
    visr::T
    visi::T
    error::T
    u::T
    v::T
    time::T
    baseline::NTuple{2,Symbol}
end

function amplitude(D::EHTVisibilityDatum)
    amp = hypot(D.visr,D.visi)
    return EHTVisibilityAmplitudeDatum(amp, D.error,
                                       D.u, D.v,
                                       D.time,
                                       D.frequency,
                                       D.bandwidth,
                                       D.baseline
                                    )
end

Base.@kwdef struct EHTVisibilityAmplitudeDatum{T<:Number} <: AbstractVisibilityDatum{T}
    amp::T
    error::T
    u::T
    v::T
    time::T
    baseline::NTuple{2,Symbol}
end

function checktriangle(D1::EHTVisibilityDatum,
                       D2::EHTVisibilityDatum,
                       D3::EHTVisibilityDatum)
    b1 = D1.baseline
    b2 = D2.baseline
    b3 = D3.baseline
    l = length(unique([b1..., b2..., b3...]))
    @assert l == 3 "For a valid closure phase you need 3 unique stations not $l"
    @assert (D1.time == D2.time == D3.time) "For a valid closure phase the times need to match"

end

@inline function visibility(D::EHTVisibilityDatum{T}) where {T}
    return Complex{T}(D.visr, D.visi)
end

@inline function amplitude(D::EHTVisibilityAmplitudeDatum{T}) where {T}
    return D.amp
end

@inline function bispectrum(D1::EHTVisibilityDatum, D2::EHTVisibilityDatum, D3::EHTVisibilityDatum)
    visibility(D1)*visibility(D2)*visibility(D3)
end

"""
    $(TYPEDEF)

A Datum for a single closure phase. Note in the future this may get replaced
with the fully covariant formalism from Blackburn et al. (2020).
"""
Base.@kwdef struct EHTClosurePhaseDatum{T<:Number} <: ClosureProducts{T}
    phase::T
    error::T
    u1::T
    v1::T
    u2::T
    v2::T
    u3::T
    v3::T
    time::T
    triangle::NTuple{3,Symbol}
end

function stations(d::EHTObservation{T,A}) where {T,A<:EHTClosurePhaseDatum}
    bl = getdata(d, :triangle)
    return unique(vcat(collect.(bl)...))
end

"""
    $(SIGNATURES)
Computes the closure phase of the three visibility datums.

# Notes
We currently use the high SNR Gaussian error approximation for the closure phase.
In the future we may use the moment matching from Monte Carlo sampling.
"""
function closure_phase(D1::EHTVisibilityDatum,
                      D2::EHTVisibilityDatum,
                      D3::EHTVisibilityDatum)

    checktriangle(D1,D2,D3)

    amp1 = amplitude(D1).amp
    amp2 = amplitude(D2).amp
    amp3 = amplitude(D3).amp
    u1,v1 = uvpositions(D1)
    u2,v2 = uvpositions(D2)
    u3,v3 = uvpositions(D3)

    bis = bispectrum(D1, D2, D3)
    s12 = unique([D1.baseline..., D2.baseline...])
    s123 = unique([s12..., D3.baseline...])
    #Use the Gaussian approximation TODO hook this into Measurements.jl?
    err = sqrt((D1.error/amp1)^2 + (D2.error/amp2)^2 + (D3.error/amp3)^2)
    return EHTClosurePhaseDatum(angle(bis), err,
                                u1, v1, u2, v2, u3, v3,
                                time, s123)
end

function baselines(CP::EHTClosurePhaseDatum)
    tri = CP.triangle
    return ((tri[1],tri[2]), (tri[2], tri[3]), (tri[3], tri[1]))
end

uvpositions(datum::EHTClosurePhaseDatum) = (datum.u1, datum.v1, datum.u2, datum.v2, datum.u3, datum.v3)

Base.@kwdef struct EHTLogClosureAmplitudeDatum{T<:Number} <: ClosureProducts{T}
    amp::T
    error::T
    u1::T
    v1::T
    u2::T
    v2::T
    u3::T
    v3::T
    u4::T
    v4::T
    time::T
    quadrangle::NTuple{4,Symbol}
end

function baselines(CP::EHTLogClosureAmplitudeDatum)
    quad = CP.quadrangle
    return ((quad[1],quad[2]), (quad[3], quad[4]), (quad[1], quad[3]), (quad[2], quad[4]))
end

function stations(d::EHTObservation{T,A}) where {T,A<:EHTLogClosureAmplitudeDatum}
    bl = getdata(d, :quadrangle)
    return unique(vcat(collect.(bl)...))
end

uvpositions(datum::EHTLogClosureAmplitudeDatum) = (datum.u1, datum.v1, datum.u2, datum.v2, datum.u3, datum.v3, datum.u4, datum.v4)

# @doc raw"""
#     $(SIGNATURES)
# Compute the log-closure amplitude of 4 visibility datums.
# This uses the same ordering as eht-imaging i.e.

# ```math
# C_A = \log \frac{|V_{1}| |V_{2}|}{|V_{3}| |V_{4}|}
# ```

# # Notes

# Currently we use the high SNR Gaussian approximation to the thermal error.
# In the future we may use the moment matching approach using Monte Carlo sampling.
# """
# function logclosure_amplitude(D1::EHTVisibilityDatum,
#                               D2::EHTVisibilityDatum,
#                               D3::EHTVisibilityDatum,
#                               D4::EHTVisibilityDatum)
#     check_quadrangle(D1,D2,D3,D4)

#     amp1 = amplitude(D1)
#     amp2 = amplitude(D2)
#     amp3 = amplitude(D3)
#     amp4 = amplitude(D4)

#     u1,v1 = uvpositions(D1)
#     u2,v2 = uvpositions(D2)
#     u3,v3 = uvpositions(D3)
#     u4,v4 = uvpositions(D4)

#     #Construct the quadrangle
#     s12 = unique([D1.baseline..., D2.baseline...])
#     s123 = unique([s12..., D3.baseline...])
#     s1234 = unique([s123..., D4.baseline...])

#     # Now do the error propogation
#     err = hypot(D1.err/map1, D2.err/amp2, D3.err/amp3, D4.err/amp4)

#     #construct the lcamp
#     lcamp = log(amp1*amp2/(amp3*amp4))

#     return EHTLogClosureAmplitudeDatum(lcamp, err,
#                                        u1, v1,
#                                        u2, v2,
#                                        u3, v3,
#                                        u4, v4,
#                                        time, s1234)
# end

function arrayconfig(vis::EHTObservation{F,A}) where {F,A<:Union{EHTVisibilityDatum, EHTVisibilityAmplitudeDatum}}
    u = getdata(vis, :u)
    v = getdata(vis, :v)
    st = stations(vis)
    times = getdata(vis, :time)
    bandwidth = vis.bandwidth
    frequency = vis.frequency
    error = getdata(vis, :error)
    baseline = getdata(vis, :baseline)
    uvsamples = StructArray{ArrayBaselineDatum}(time=times,
                                        u=u,
                                        v=v,
                                        baseline=baseline,
                                        error_real=error,
                                        error_imag=error
                                    )
    return EHTArrayConfiguration(st, frequency, bandwidth, uvsamples)
end

include(joinpath(@__DIR__, "io.jl"))
include(joinpath(@__DIR__, "ehtim.jl"))


text:
# kron.jl
# test

using LinearMapsAA: LinearMapAA, LinearMapAM, LinearMapAO
using LinearAlgebra: I, Diagonal
using Test: @test, @testset

M1 = rand(ComplexF64, 3, 3)
M2 = rand(ComplexF64, 2, 2)

M12 = kron(M1, M2)
A1 = LinearMapAA(M1)
A2 = LinearMapAA(M2)
@testset "kron basics" begin
    @test kron(A1, A1) isa LinearMapAM
    for pair in ((A1,A2), (A1,M2), (M1,A2)) # all combinations
        AAk = kron(pair[1], pair[2])
        @test AAk isa LinearMapAM
        @test Matrix(AAk) == M12
        @test Matrix(AAk)' == Matrix(AAk')
    end
end

A3 = LinearMapAA(M1 ; odim=(1,size(M1,1)))
A4 = LinearMapAA(M2 ; idim=(size(M2,2),1))
@testset "kron(I,A)" begin
    K = kron(I(2), A4)
    M = kron(I(2), M2)
    @test K isa LinearMapAO
       @test Matrix(K) == M
       @test Matrix(K') == Matrix(K)'
end

@testset "kron(A,I)" begin
    K = kron(A4, I(2))
    M = kron(M2, I(2))
    @test K isa LinearMapAO
       @test Matrix(K) == M
       @test Matrix(K') == Matrix(K)'
end

@testset "kron(Ao,M)" begin
    @test kron(A3, M2) isa LinearMapAO
    @test kron(M1, A4) isa LinearMapAO
end

@testset "kron(Ao,Bo)" begin
    K = kron(A3, A4)
    M = kron(M1, M2)
    @test K isa LinearMapAO
       @test Matrix(K) == M
       @test Matrix(K') == Matrix(K)'
end

@testset "kron(A,D)" begin
    D = Diagonal(1:3)
    K = kron(A1, D)
    @test Matrix(K) == kron(M1, D)
    K = kron(D, A1)
    @test Matrix(K) == kron(D, M1)
end


text:
# Example julia lambda function
# Couting words in a file in S3 path

# See following resources for S3 in julia
# - https://juliacloud.github.io/AWSCore.jl/build/AWSS3.html
# - https://github.com/JuliaCloud/AWSS3.jl

module word_count

using AWSLambdaJuliaRuntime
using AWSS3
using AWSCore

#= The handler function =#
function handler(event_data::InvocationRequest)
	# download some file from S3 path
	# count number of unique words
	# put result back to S3 path or send in response

	println("Hello World : $(event_data.payload)")
	run(`ls`)

	i_am_good = true
	if i_am_good
		return success_invocation_response("""{"msg": "Me Rockz!"}""", "application/json")
	else
		return failure_invocation_response("Me Suckz!", "some_error_type")
	end
end

end # module word_count

max_stars_repo_name: FugroRoames/aws-lambda-julia-runtime

max_issues_repo_name: FugroRoames/aws-lambda-julia-runtime

max_forks_repo_name: FugroRoames/aws-lambda-julia-runtime


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
##

using DataFrames
using CSV
using Plots
using StatsPlots
using MCPhylo
using ProgressMeter
##

fpairs = CSV.read("../data/fpairs.txt", DataFrame, header=false)[:,1]

##

theme(:solarized_light)
#
upscale = 1 #8x upscaling in resolution
fntsm = Plots.font("sans-serif", pointsize=round(20.0*upscale))
fntlg = Plots.font("sans-serif", pointsize=round(25.0*upscale))
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
default(size=(800*upscale,600*upscale)) #Plot canvas size

##

d = CSV.read("ppp/empiricalStats.csv", DataFrame)

d[:,:char] = 1:28
d[:,:method] .= "empirical"
d[:,:model] .= "nature"

for charNum in 1:28
    global d
    cdf = CSV.read("ppp/$(lpad(charNum, 2, "0")).csv", DataFrame)
    cdf[:,:char] .= charNum
    d = vcat(d, cdf)
end

##
n=28
##

tvPlots = []

for charNum in 1:n
    cDF = filter(x -> x.char == charNum, d)
    cl = filter(x -> x.model=="lineageSpecific", cDF).totalVariance
    cu = filter(x -> x.model=="universal", cDF).totalVariance
    eu = filter(x -> x.model=="nature", cDF).totalVariance
    p = plot([1,1], hpd(cu),
        size=(4*600, 400*7),
        linewidth=5,
        xlim=[0.5,2.5],
        legend=false,
        color=:black,
        xrotation=0,
        xticks=([1,2], ["", ""]),
        ylim=(0.27,0.85)
    )
    if charNum%4 != 2
        yticks!(p,[false])
    end
    plot!([1,1], hpd(cu, alpha=0.5), linewidth=15, color=:black, alpha=0.5)
    plot!([2,2], hpd(cl), linewidth=5, color=:black)
    plot!([2,2], hpd(cl, alpha=0.5), linewidth=15, color=:black, alpha=0.5)
    plot!(p,[.9,2.1],[eu[1],eu[1]], linewidth=2)
    title!(fpairs[charNum])
    push!(tvPlots, p)
end
##
xticks!(tvPlots[end-3], ([1,2], ["universal", "lineage-specific"]))
xticks!(tvPlots[end-2], ([1,2], ["universal", "lineage-specific"]))
xticks!(tvPlots[end-1], ([1,2], ["universal", "lineage-specific"]))
xticks!(tvPlots[end], ([1,2], ["universal", "lineage-specific"]))
p1 = plot(tvPlots...,
    layout=grid(7,4, heights=repeat([1/7],n)))

savefig("../data/img/pppTV.pdf")
##

lvPlots = []

for charNum in 1:n
    cDF = filter(x -> x.char == charNum, d)
    cl = filter(x -> x.model=="lineageSpecific", cDF).lineagewiseVariance
    cu = filter(x -> x.model=="universal", cDF).lineagewiseVariance
    eu = filter(x -> x.model=="nature", cDF).lineagewiseVariance
    p = plot([1,1], hpd(cu),
        size=(4*600, 400*7),
        linewidth=5,
        xlim=[0.5,2.5],
        legend=false,
        color=:black,
        xrotation=0,
        xticks=([1,2], ["", ""]),
        ylim=(0.03,0.16),
        )
    if charNum%4 != 2
        yticks!(p,[false])
    end
    plot!([1,1], hpd(cu, alpha=0.5), linewidth=15, color=:black, alpha=0.5)
    plot!([2,2], hpd(cl), linewidth=5, color=:black)
    plot!([2,2], hpd(cl, alpha=0.5), linewidth=15, color=:black, alpha=0.5)
    plot!(p,[.9,2.1],[eu[1],eu[1]], linewidth=2)
    title!(fpairs[charNum])
    push!(lvPlots, p)
end
##
xticks!(lvPlots[end-3], ([1,2], ["universal", "lineage-specific"]))
xticks!(lvPlots[end-2], ([1,2], ["universal", "lineage-specific"]))
xticks!(lvPlots[end-1], ([1,2], ["universal", "lineage-specific"]))
xticks!(lvPlots[end], ([1,2], ["universal", "lineage-specific"]))
p1 = plot(lvPlots...,
    layout=grid(7,4, heights=repeat([1/7],n)))

savefig("../data/img/pppLV.pdf")
##

clvPlots = []

for charNum = 1:n
    cDF = filter(x -> x.char == charNum, d)
    cl = filter(x -> x.model == "lineageSpecific", cDF).crossLineageVariance
    cu = filter(x -> x.model == "universal", cDF).crossLineageVariance
    eu = filter(x -> x.model == "nature", cDF).crossLineageVariance
    p = plot(
        [1, 1],
        hpd(cu),
        size=(4*600, 400*7),
        linewidth = 5,
        xlim = [0.5, 2.5],
        legend = false,
        color = :black,
        xrotation = 0,
        xticks = ([1, 2], ["", ""]),
        ylim = (0.28, 0.7),
    )
    if charNum % 4 != 2
        yticks!(p, [false])
    end
    plot!([1, 1], hpd(cu, alpha = 0.5), linewidth = 15, color = :black, alpha = 0.5)
    plot!([2, 2], hpd(cl), linewidth = 5, color = :black)
    plot!([2, 2], hpd(cl, alpha = 0.5), linewidth = 15, color = :black, alpha = 0.5)
    plot!(p, [0.9, 2.1], [eu[1], eu[1]], linewidth = 2)
    title!(fpairs[charNum])
    push!(clvPlots, p)
end
##
xticks!(clvPlots[end-3], ([1,2], ["universal", "lineage-specific"]))
xticks!(clvPlots[end-2], ([1,2], ["universal", "lineage-specific"]))
xticks!(clvPlots[end-1], ([1,2], ["universal", "lineage-specific"]))
xticks!(clvPlots[end], ([1,2], ["universal", "lineage-specific"]))
p1 = plot(clvPlots...,
    layout=grid(7,4, heights=repeat([1/7],n)))

savefig("../data/img/pppCLV.pdf")
##

function getStatistics(charNum, stats)
    df = select(filter(x -> x.char == charNum, d), [stats, :model])
    r = [filter(x -> x.model == "nature", df)[1,stats]]
    append!(r, hpd(filter(x -> x.model == "universal", df)[:,stats]))
    append!(r, hpd(filter(x -> x.model == "lineageSpecific", df)[:,stats]))
end

##

tv = DataFrame(Tuple.(getStatistics.(1:28, :totalVariance)))
rename!(
    tv,
    [:empirical, :universal_lower, :universal_upper, :lineage_lower, :lineage_upper],
)

insertcols!(tv, 1, :fpair => fpairs)

e = Array(rpad.(round.(tv[:,2], digits = 3), 5, "0"))

uHPD = mapslices(
    x -> "(" * join(x, ", ") * ")",
    Array(rpad.(round.(tv[:, 3:4], digits = 3), 5, "0")),
    dims = 2,
)

lHPD = mapslices(
    x -> "(" * join(x, ", ") * ")",
    Array(rpad.(round.(tv[:, 5:6], digits = 3), 5, "0")),
    dims = 2,
)

output = DataFrame(
    "feature pair" => tv.fpair,
    "empirical" => e,
    "universal model (HPD)" => uHPD[:,1],
    "lineage-specific model (HPD)" => lHPD[:,1],
)

CSV.write("../data/tables/tv.tex", output, newline="\\\\\n", delim="&")

##

lv = DataFrame(Tuple.(getStatistics.(1:28, :lineagewiseVariance)))
rename!(
    lv,
    [:empirical, :universal_lower, :universal_upper, :lineage_lower, :lineage_upper],
)

insertcols!(lv, 1, :fpair => fpairs)

e = Array(rpad.(round.(lv[:,2], digits = 3), 5, "0"))

uHPD = mapslices(
    x -> "(" * join(x, ", ") * ")",
    Array(rpad.(round.(lv[:, 3:4], digits = 3), 5, "0")),
    dims = 2,
)

lHPD = mapslices(
    x -> "(" * join(x, ", ") * ")",
    Array(rpad.(round.(lv[:, 5:6], digits = 3), 5, "0")),
    dims = 2,
)

output = DataFrame(
    "feature pair" => lv.fpair,
    "empirical" => e,
    "universal model (HPD)" => uHPD[:,1],
    "lineage-specific model (HPD)" => lHPD[:,1],
)

CSV.write("../data/tables/lv.tex", output, newline="\\\\\n", delim="&")
##

clv = DataFrame(Tuple.(getStatistics.(1:28, :crossLineageVariance)))
rename!(
    clv,
    [:empirical, :universal_lower, :universal_upper, :lineage_lower, :lineage_upper],
)

insertcols!(clv, 1, :fpair => fpairs)

e = Array(rpad.(round.(clv[:,2], digits = 3), 5, "0"))

uHPD = mapslices(
    x -> "(" * join(x, ", ") * ")",
    Array(rpad.(round.(clv[:, 3:4], digits = 3), 5, "0")),
    dims = 2,
)

lHPD = mapslices(
    x -> "(" * join(x, ", ") * ")",
    Array(rpad.(round.(clv[:, 5:6], digits = 3), 5, "0")),
    dims = 2,
)

output = DataFrame(
    "feature pair" => clv.fpair,
    "empirical" => e,
    "universal model (HPD)" => uHPD[:,1],
    "lineage-specific model (HPD)" => lHPD[:,1],
)

CSV.write("../data/tables/clv.tex", output, newline="\\\\\n", delim="&")

##

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
module BenchKMeans
using Random
using ParallelKMeans
using BenchmarkTools

suite = BenchmarkGroup()

Random.seed!(2020)
X = rand(10, 100_000)

centroids3 = ParallelKMeans.smart_init(X, 3, 1, init="kmeans++").centroids
centroids10 = ParallelKMeans.smart_init(X, 10, 1, init="kmeans++").centroids

suite["10x100_000x3x1     Lloyd"] = @benchmarkable kmeans($X, 3, init = $centroids3, n_threads = 1, verbose = false, tol = 1e-6, max_iters = 1000)
suite["10x100_000x3x1  Hammerly"] = @benchmarkable kmeans(Hamerly(), $X, 3, init = $centroids3, n_threads = 1, verbose = false, tol = 1e-6, max_iters = 1000)

suite["10x100_000x3x2     Lloyd"] = @benchmarkable kmeans($X, 3, init = $centroids3, n_threads = 2, verbose = false, tol = 1e-6, max_iters = 1000)
suite["10x100_000x3x2  Hammerly"] = @benchmarkable kmeans(Hamerly(), $X, 3, init = $centroids3, n_threads = 2, verbose = false, tol = 1e-6, max_iters = 1000)

suite["10x100_000x10x1    Lloyd"] = @benchmarkable kmeans($X, 10, init = $centroids10, n_threads = 1, verbose = false, tol = 1e-6, max_iters = 1000)
suite["10x100_000x10x1 Hammerly"] = @benchmarkable kmeans(Hamerly(), $X, 10, init = $centroids10, n_threads = 1, verbose = false, tol = 1e-6, max_iters = 1000)

suite["10x100_000x10x2    Lloyd"] = @benchmarkable kmeans($X, 10, init = $centroids10, n_threads = 2, verbose = false, tol = 1e-6, max_iters = 1000)
suite["10x100_000x10x2 Hammerly"] = @benchmarkable kmeans(Hamerly(), $X, 10, init = $centroids10, n_threads = 2, verbose = false, tol = 1e-6, max_iters = 1000)

end # module

BenchKMeans.suite


text:
# Desired:
#   using Test
#
# In order that Test need not to be defined in the user environment, it is included via Modia:
import Modia

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

@testset "RunSimulations" begin

    Modia.ModiaLogging.resetTestStatus()

    include("models/TestVariableTypes.jl")
    include("models/TestEquations.jl")

    include("models/TestElectrical.jl")
    include("models/TestFilter.jl")
    include("models/TestArrayOfComponents.jl")
    include("models/TestConditionalComponents.jl")
    include("models/TestConditionalEquations.jl")

    include("models/TestSpatialDiscretization.jl")

    include("models/MergingModifiers.jl")

    #include("models/TestFluid.jl")

    #include("models/TestCoupledInertias.jl")
    #include("models/TestPendulum.jl")
    #include("models/TestStateSelection.jl")
    #include("models/TestRetranslationOmega.jl")

    Modia.ModiaLogging.printTestStatus()

end
nothing


text:
import Base: size, length, isless, getindex

export
    BitVar,
    BitVarArray,
    size,
    indexes,
    name,
    getindex,
    BitAssignment

name(v::AbstractVariable) = v.name

struct BitVar <: AbstractBitVariable
    name::String
    function BitVar(name::String)
        @assert length(name) > 0
        new(name)
    end
end

size(::BitVar) = 1 # BitVar always has one element
length(::BitVar) = 1 # BitVar always has one element

struct BitVarArray  <: AbstractBitVariable
    name::String
    size::Vector{Int}
    function BitVarArray(name::String,
                         size::Vector{Int})
        @assert name != ""
        @assert length(size) > 0
        @assert all(0 .< size)
        new(name, size)
    end
end

size(v::BitVarArray) = v.size
length(v::BitVarArray) = prod(size(v))

struct BitVarSubArray  <: AbstractBitVariable
    var::BitVarArray
    indexes::Vector{Int}
    function BitVarSubArray(var::BitVarArray,
                            indexes::Vector{Int})
        @assert length(indexes) == length(size(var))
        @assert all(0 .<= indexes .<= size(var)) # 0 means "all"
        new(var, indexes)
    end
end

name(v::BitVarSubArray) = name(v.var)
size(v::BitVarSubArray) = size(v.var)
length(v::BitVarSubArray) = prod(size(v))
indexes(v::BitVarSubArray) = v.indexes

# for syntax:
# var = BitVarArray("mine", [3,2,3])
# var2 = var[:,2,:]
# TODO syntax var[:,:,[1,3]] or var[:,:,1:2]
function Base.getindex(var::BitVarArray, I::Vararg)
    I = map(x-> ifelse(x == Colon(), 0, x), I)
    BitVarSubArray(var, Int[I...])
end

isless(v::AbstractBitVariable, w::AbstractBitVariable) = (name(v) < name(w))


text:
module G2
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    1.19041545916754e-6*l1.^5 + 8.22231271886465e-6*l1.^4.*l2 - 1.41790652493433e-5*l1.^4.*l3 - 3.71236760805448e-8*l1.^4 + 1.87237365187196e-5*l1.^3.*l2.^2 - 7.03670987705473e-5*l1.^3.*l2.*l3 - 9.63717125914612e-8*l1.^3.*l2 + 6.3566217782067e-5*l1.^3.*l3.^2 + 2.45015054239551e-7*l1.^3.*l3 - 1.40800538973871e-10*l1.^3 + 1.60903292357735e-5*l1.^2.*l2.^3 - 0.000104503275318621*l1.^2.*l2.^2.*l3 + 4.83183384860215e-8*l1.^2.*l2.^2 + 0.000210160563816057*l1.^2.*l2.*l3.^2 + 1.92814428599642e-7*l1.^2.*l2.*l3 - 4.2295310009243e-10*l1.^2.*l2 - 0.000133679823732731*l1.^2.*l3.^3 - 4.64320735546665e-7*l1.^2.*l3.^2 + 8.45817472956316e-10*l1.^2.*l3 - 1.55805030397592e-13*l1.^2 + 5.86039594070883e-6*l1.*l2.^4 - 5.56603856920193e-5*l1.*l2.^3.*l3 + 8.31924615763953e-8*l1.*l2.^3 + 0.000188112147098247*l1.*l2.^2.*l3.^2 - 3.46504541199643e-7*l1.*l2.^2.*l3 - 6.72525431864125e-11*l1.*l2.^2 - 0.000265665096331678*l1.*l2.*l3.^3 + 1.53644624301111e-7*l1.*l2.*l3.^2 + 9.81678669090356e-10*l1.*l2.*l3 - 4.6644143835908e-13*l1.*l2 + 0.000133323717219143*l1.*l3.^4 + 2.58608071622046e-7*l1.*l3.^3 - 1.33776521557279e-9*l1.*l3.^2 + 7.78418394428497e-13*l1.*l3 + 1.14158920155552e-5*l1 + 7.69886693605312e-7*l2.^5 - 9.71720191286928e-6*l2.^4.*l3 + 2.04003115047211e-8*l2.^4 + 4.72984134088345e-5*l2.^3.*l3.^2 - 1.64978151138711e-7*l2.^3.*l3 + 3.67917494575423e-11*l2.^3 - 0.000110075710259202*l2.^2.*l3.^3 + 4.21141402372243e-7*l2.^2.*l3.^2 - 4.32020406856046e-11*l2.^2.*l3 - 8.45245878638712e-14*l2.^2 + 0.000121528242051865*l2.*l3.^4 - 3.32241916958472e-7*l2.*l3.^3 - 4.48196555086032e-10*l2.*l3.^2 + 6.36432165987811e-13*l2.*l3 + 4.68195293026799e-6*l2 - 5.09987206760422e-5*l3.^5 + 1.84060408137646e-8*l3.^4 + 5.95882149068005e-10*l3.^3 - 7.08086905029087e-13*l3.^2 - 1.61426797114482e-5*l3 + 6.92950794258085e-9
end
end


text:
@testset "Sequence" begin
    𝒯 = Taylor(1)
    ℱ = Fourier(1, 1.0)
    𝒞 = Chebyshev(2)
    𝑇 = 𝒯 ⊗ ℱ ⊗ 𝒞
    𝕊 = 𝑇^1
    coeffs = vec([inv(prod((2, 2, 2) .^ abs.(α))) for α ∈ indices(𝑇)])
    a = Sequence(𝕊, coeffs)
    #
    @test space(a) == 𝕊
    @test coefficients(a) == coeffs
    @test order(a) == order(𝕊)
    @test order(a, 1) == order(𝕊, 1)
    @test_throws MethodError frequency(a)
    @test frequency(component(a, 1), 2) == frequency(𝑇, 2)
    @test firstindex(a) == RadiiPolynomial._firstindex(𝕊)
    @test lastindex(a) == RadiiPolynomial._lastindex(𝕊)
    @test length(a) == length(coeffs)
    @test size(a) == size(coeffs)
    @test iterate(a) == iterate(coeffs)
    @test iterate(a, 2) == iterate(coeffs, 2)
    @test eltype(a) == eltype(typeof(a)) == eltype(coeffs)
    @test a[1:3] == coeffs[1:3]
    @test view(a, 1:3) == view(coeffs, 1:3)
    @test setindex!(copy(a), 0.0, 1) == setindex!(copy(coeffs), 0.0, 1)
    @test a == a ≈ a
    @test !iszero(a)
    @test copy(a) == a
    @test iszero(zero(a))
    @test float(a) == Sequence(𝕊, float(coeffs))
    @test complex(a) == Sequence(𝕊, complex(coeffs))
    @test real(a) == Sequence(𝕊, real(coeffs))
    @test imag(a) == Sequence(𝕊, imag(coeffs))
    @test conj(a) == Sequence(𝕊, conj(coeffs))
    @test conj!(copy(a)) == Sequence(𝕊, conj!(copy(coeffs)))
    #
    a_𝑇 = component(a, 1)
    @test coefficients(a_𝑇) == coeffs
    @test selectdim(a_𝑇, 2, 0) == selectdim(reshape(coeffs, 2, 3, 3), 2, 2)
end


text:
# # [Density and collision frequency as interpolating functions](@id interpolating-functions)
#
# The ionospheric [`Species`](@ref) type used when defining a [`HomogeneousWaveguide`](@ref)
# requires `numberdensity` and `collisionfrequency` fields to be callable functions of
# altitude `z` in meters that return `Float64` values of number density in m⁻³ and
# species-neutral collision frequency in s⁻¹.
#
# Sometimes a simple analytical function can be provided, such as [`waitprofile`](@ref).
# Other times, the profile is constructed from tabular data at discrete heights which must
# be wrapped in an interpolator before being passed to `Species`.
#
# The choice of interpolating function not only impacts the accuracy of the propagation
# result (how well the interpolation resembles the true profile), but it can also impact
# the runtime of [`propagate`](@ref) and related functions.
#
# In this example we'll compare a few different interpolating functions.
#
# ## Profiles and interpolators
#
# We'll use the [FIRI-2018](https://doi.org/10.1029/2018JA025437) profiles from
# [FIRITools.jl](https://github.com/fgasdia/FIRITools.jl) which are slightly more
# complicated than a pure exponential profile. FIRITools isn't registered, but can be
# installed from the Pkg by copy-pasting the entire url.
#
# From [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) we'll use
# a linear interpolator, a cubic spline interpolator, and monotonic interpolators.
# We'll use [NormalHermiteSplines.jl](https://github.com/IgorKohan/NormalHermiteSplines.jl)
# to construct Hermite splines.

using Printf, Statistics
using LongwaveModePropagator
using LongwaveModePropagator: QE, ME
using Plots, Distances

using FIRITools

using Interpolations, NormalHermiteSplines
nothing #hide

# Here's the discrete profile data. Some interpolators will benefit from random sample
# points, but most density data will be on a grid. The profile uses an exponential
# extrapolation of the base of FIRI from about 60 km altitude down to the ground.

zs = 0:1e3:110e3
Ne = FIRITools.extrapolate(firi(50, 30), zs);

# Let's construct the interpolators.

linear_itp = LinearInterpolation(zs, Ne)
cubic_itp = CubicSplineInterpolation(zs, Ne);

# There's not great documentation on the monotonic interpolators of Interpolations.jl as of
# `v0.13`, but several are supported.

fb_itp = Interpolations.interpolate(zs, Ne, FritschButlandMonotonicInterpolation())
fc_itp = Interpolations.interpolate(zs, Ne, FritschCarlsonMonotonicInterpolation())
s_itp = Interpolations.interpolate(zs, Ne, SteffenMonotonicInterpolation());

# And the Hermite splines

spline = prepare(collect(zs), RK_H1())
spline = construct(spline, Ne)
hermite_itp(z) = evaluate_one(spline, z);

# (a spline built with RK_H0 kernel is a continuous function,
#  a spline built with RK_H1 kernel is a continuously differentiable function,
#  a spline built with RK_H2 kernel is a twice continuously differentiable function).

zs_fine = 40e3:100:110e3
Ne_fine = FIRITools.extrapolate(firi(50, 30), zs_fine);

linear_fine = linear_itp.(zs_fine)
cubic_fine = cubic_itp.(zs_fine)
fb_fine = fb_itp.(zs_fine)
fc_fine = fc_itp.(zs_fine)
s_fine = s_itp.(zs_fine)
hermite_fine = hermite_itp.(zs_fine);

# The profiles are compared using percentage difference relative to the true profile.

cmp(a,b) = (a - b)/b*100

dNe = cmp.(Ne_fine, Ne_fine)
dlinear = cmp.(linear_fine, Ne_fine)
dcubic = cmp.(cubic_fine, Ne_fine)
dfb = cmp.(fb_fine, Ne_fine)
dfc = cmp.(fc_fine, Ne_fine)
ds = cmp.(s_fine, Ne_fine)
dhermite = cmp.(hermite_fine, Ne_fine);

# To plot densities with a log scale, we set values less than 0.1 to NaN.

cl(x) = replace(v->v <= 0.1 ? NaN : v, x)
lc(x) = replace(x, NaN => 0)

p1 = plot(cl([Ne_fine linear_fine cubic_fine hermite_fine fb_fine fc_fine s_fine]),
    zs_fine/1000, xscale=:log10, xlabel="Ne (m⁻³)", ylabel="Altitude (km)",
    legend=:topleft, labels=["Truth" "Linear" "Cubic" "Hermite" "FritschButland" "FritschCarlson" "Steffen"])
p2 = plot(lc([dNe dlinear dcubic dhermite dfb dfc ds]),
    zs_fine/1000, xlabel="% difference", legend=false, xlims=(-1, 1))
plot(p1, p2, layout=(1,2), size=(800,400), margin=3Plots.mm)
#md savefig("interpolatingfunctions_profiles.png"); nothing # hide
#md # ![](interpolatingfunctions_profiles.png)

# Unsurprisingly, the error is highest at the cutoff altitude of 40 km where the densities
# below are 0.
# The linear interpolation has positive-biased errors at all other heights because linear
# interpolation does not accurately capture the true exponential profile.
# To avoid this, the interpolator could have been constructed with a finer initial grid, but
# that is not always possible.
#
# Let's compute the total absolute difference between each interpolator and the truth and
# also compute the average percentage difference for each:

for (n, v) in ("linear"=>(linear_fine, dlinear), "cubic"=>(cubic_fine, dcubic),
    "hermite"=>(hermite_fine, dhermite), "FritschButland"=>(fb_fine, dfb),
    "FritschCarlson"=>(fc_fine, dfc), "Steffen"=>(s_fine, ds))
    @printf("%s:  %.3g  %.3g\n", n, cityblock(Ne_fine, v[1]), mean(abs, v[2]))
end

# Somewhat surprisingly, the linear interpolation has the lowest total absolute difference
# from the truth despite the fact that (as seen in the plot) the percentile difference from
# the linear interpolation to the truth is very high.
#
# ## Propagation results
#
# Let's compare propagation results and performance with each of these interpolators.
# Unfortunately, as of `v0.5.1` of NormalHermiteSplines.jl, the implementation of the
# package isn't efficient (see #3)[https://github.com/IgorKohan/NormalHermiteSplines.jl/issues/3]
# and because `hermite_itp` would be called millions of times in `propagate`, it is
# prohibitively slow.
#
# Here are the `Species` for each of the other interpolators.
# They will all use the analytic [`electroncollisionfrequency`](@ref).
#
# There's no functional form of the truth FIRI profile, so we'll build a fine FritschButland
# interpolator and call it the truth.

interpolators = (
    "truth" => Interpolations.interpolate(0:100:110e3,
        FIRITools.extrapolate(firi(50, 30), 0:100:110e3),
        FritschButlandMonotonicInterpolation()),
    "linear"=> linear_itp,
    "cubic" => cubic_itp,
    "FritschButland" => fb_itp,
    "FritschCarlson" => fc_itp,
    "Steffen" => s_itp
)

function propagateitp(interpolators)
    bfield = BField(50e-6, deg2rad(68), deg2rad(111))
    ground = GROUND[5]

    tx = Transmitter(24e3)
    rx = GroundSampler(0:5e3:3000e3, Fields.Ez)

    results = Dict{String,Tuple{Float64,Vector{Float64},Vector{Float64}}}()
    for (n, itp) in interpolators
        species = Species(QE, ME, itp, electroncollisionfrequency)
        waveguide = HomogeneousWaveguide(bfield, species, ground)
        t0 = time()
        _, amp, phase = propagate(waveguide, tx, rx)
        runtime = time() - t0
        results[n] = (runtime, amp, phase)
    end
    return results
end

propagateitp(interpolators);  # warmup
results = propagateitp(interpolators);

# First let's evaluate the propagation results.

d = 0:5:3000
p1 = plot(ylabel="Amplitude (dB μV/m)")
p2 = plot(xlabel="Range (km)", ylims=(-0.02, 0.02), ylabel="Δ", legend=false)
for (n, v) in results
    plot!(p1, d, v[2], label=n)
    plot!(p2, d, v[2]-results["truth"][2])
end
plot(p1, p2, layout=grid(2,1,heights=[0.7, 0.3]))
#md savefig("interpolatingfunctions_amp.png"); nothing # hide
#md # ![](interpolatingfunctions_amp.png)

p1 = plot(ylabel="Phase (deg)")
p2 = plot(xlabel="Range (km)", ylims=(-0.4, 0.4), ylabel="Δ", legend=false)
for (n, v) in results
    plot!(p1, d, rad2deg.(v[3]), label=n)
    plot!(p2, d, rad2deg.(v[3])-rad2deg.(results["truth"][3]))
end
plot(p1, p2, layout=grid(2,1,heights=[0.7, 0.3]))
#md savefig("interpolatingfunctions_phase.png"); nothing # hide
#md # ![](interpolatingfunctions_phase.png)

# Here is the mean absolute amplitude difference between each technique and the true profile
# amplitude:

for n in ("linear", "cubic", "FritschButland", "FritschCarlson", "Steffen")
    @printf("%s: %.3e\n",n, meanad(results[n][2], results["truth"][2]))
end

# The amplitude and phase for each of the interpolators matches the true exponential profile
# extremely closely relative to the typical noise of real VLF measurements.
#
# Timing can vary when run by GitHub to build the documentation, so results here are from a
# local run:
#
# |  Interpolator  | Runtime relative to `truth` |
# | -------------- | --------------------------- |
# | truth          |              1              |
# | linear         |             0.99            |
# | cubic          |             0.84            |
# | FritschButland |             0.94            |
# | FritschCarlson |             0.96            |
# | Steffen        |             0.95            |
#
# Really, any of these except for the linear interpolation could be used to interpolate a
# discrete profile. The cubic interpolation is the fastest, but the FritschButland is a
# little more true to the actual profile. It's for that reason that FritschButland is used
# to interpolate [`TableInput`](@ref) types in LongwaveModePropagator.


text:
"""
The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?
"""
# Running isprime will have to test prime each time it is called
# It might be better to generate a prime prime_sieve till the sqrt and check from reverse if the number is divisible.

function prime_sieve(n::Int)
    isprime = trues(n)
    isprime[1] = false

    for i in 2:isqrt(n)
        if isprime[i]
            for j in i^2:i:n
                isprime[j] = false
            end
        end
    end
    return filter(x -> isprime[x],1:n)
end

function largest_prime_factor(n)
    primes = prime_sieve(Int(trunc(√(n))))
    for prime in Iterators.reverse(primes)
        if n%prime == 0
            return prime
        end
    end
    return -1
end

@assert largest_prime_factor(13195) == 29
@show largest_prime_factor(600851475143)


text:
module Model

using Gen
using ..MangroveBase: ReviewInfo, Rating, RATINGS

@gen function subset(items::Vector)::Vector
  filter(item -> @trace(bernoulli(0.5), item[1]), items)
end

@dist subject_quality() = uniform_discrete(1, RATINGS)

@enum Bias neutral owner

@dist reviewer_bias()::Bias = Bias(categorical([0.8, 0.2]) - 1)

@gen function triangle_categorical(peak::Real)::Float64
  raw = [(5 - abs(x - peak)) for x in 1:5]
  @trace(categorical(raw ./ sum(raw)), :rating)
end

@gen function rating(quality, bias)::Rating
  if bias == neutral
    @trace(triangle_categorical(quality))
  elseif bias == owner
    if bernoulli(0.3)
      # This place is his, rate it high.
      @trace(triangle_categorical(RATINGS))
    else
      # Someone else's place, rate it low.
      @trace(triangle_categorical(1))
    end
  end
end

@gen function mangrove_model(subjects::Set{String}, reviewers::Set{String})
  biases = [id => @trace(reviewer_bias(), :bias => id) for id in reviewers]
  for sub in subjects
    true_quality = @trace(subject_quality(), :quality => sub)
    for (id, bias) in @trace(subset(biases), :selected => sub)
      @trace(rating(true_quality, bias), ReviewInfo(sub, id))
    end
  end
end

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
@with_kw struct ContinuousBanditMDP <: MDP{Float32, Float32}
    goal::Float32 = 0
    σ::Float32 = 1
end

POMDPs.reward(mdp::ContinuousBanditMDP, s, a) = exp(-(a[1] - mdp.goal)^2 / mdp.σ) - 0.1*abs(a[1] - mdp.goal)
POMDPs.gen(mdp::ContinuousBanditMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict()) = (sp = a[1], r = reward(mdp, s, a))
POMDPs.initialstate(mdp::ContinuousBanditMDP) = Deterministic(Float32(randn()))
POMDPs.isterminal(mdp::ContinuousBanditMDP, s) = true
isfailure(mdp::ContinuousBanditMDP, s) = s[1] > mdp.goal
POMDPs.discount(mdp::ContinuousBanditMDP) = 0.99
POMDPs.convert_s(::Type{V}, s::Float32, mdp::ContinuousBanditMDP) where {V<:AbstractArray} = [s]


text:
include("header.jl")

#=
# The Rayleigh problem and self-similarity

This notebook reviews the **Rayleigh problem**, the flow that develops in
a viscous fluid when an infinite flat wall is impulsively started in motion
to constant velocity $U$ tangent to the wall.

Here, we review the concepts of
- Viscous diffusion
- Self-similar flow
- Vorticity creation
=#

# ### Set up the module
using MAE250B
#-
using Plots
using SpecialFunctions

#=
To solve this flow, we assume that the flow is incompressible and two-dimensional.
There is also no $x$ dependence. From continuity and the no-flow-through condition,
the velocity normal to the wall ($v$) is exactly zero.

The Navier-Stokes equations for this problem reduce to a linear partial differential
equation for $u(y,t)$, the velocity in the $x$ direction, parallel to the wall.

$$ \dfrac{\partial u}{\partial t} = \nu \dfrac{\partial^2 u}{\partial y^2} $$

subject to the no-slip boundary condition on the wall:

$$ u(0,t) = U $$

and another condition at $\infty$, where the flow is always at rest (i.e., it
never feels the wall's motion):

$$ u(\infty,t) = 0 $$

The initial condition is that the flow is at rest:

$$ u(y,0) = 0 $$

Obviously this is at odds with the boundary condition at the wall. In fact,
this mismatch is what leads to all of the behavior that follows.
=#

#=
Before we present the solution, let's reason this out physically: as the wall
starts moving, the fluid right next to the wall will start moving, too. Because
of viscosity, this fluid, in turn, makes fluid above it start moving, and so on.
Over time, diffusion will allow "information" about this wall's motion to
spread upward.

We expect that the thickness $\delta$ of the layer of fluid that is in motion must be
dependent on time $t$ and on kinematic viscosity $\nu$. By dimensional analysis,
this thickness must be $\delta(t) \propto \sqrt{\nu t}$. It does not actually
matter what the constant of proportionality is. By convention, we
define the thickness as $\delta(t) = \sqrt{4 \nu t}$
=#

δ(t;ν = 1) = sqrt(4*ν*t)

#=
The exact solution for this problem can be derived using a *similarity transform*.
Basically, we hypothesize that the $y$ dependence of $u$ should be measured
in units of $\delta(t)$, because there is no other length scale in the problem.
We will make sense of this later. For now, we express this idea mathematically. We
propose that the solution should take the form

$$ u(y,t) = U F(\eta) $$

where $\eta = y/\delta(t)$. This $\eta$ is a dimensionless $y$ coordinate.

The $U$ factor makes sense: if we double the wall speed, we expect the
velocity in the fluid to double, as well. Note that we are also proposing that
$y$ and $t$ do not affect the velocity separately, but instead, in a mixed way.
We will explore what this means physically a little bit later.
But mathematically, we are proposing that there is only *one* independent variable $\eta$, not
two ($y$ and $t$). If we are right about this, then, by substituting this into the
original PDE for $u(y,t)$, we should get an ODE for $F(\eta)$.
=#

#=
In fact,
after doing this, applying the chain rule, and simplifying, we do get an ODE:

$$ F'' + 2 \eta F' = 0$$

subject to two conditions, $F(0) = 1$ (the "dimensionless" no-slip condition)
and $F(\infty) = 0$. The solution is

$$ F(\eta) = 1 - \mathrm{erf}(\eta) $$

where $\mathrm{erf}$ is the error function.
=#

F(η) = 1 - erf(η)
#-

#=
Let's plot this. It's a dimensionless profile, so we put $\eta$ on the
vertical axis:
=#
η = range(0.0,5.0,length=501)
plot(F.(η),η,xlims=(0,2),ylims=(0,5),xguide=L"F(\eta)",yguide=L"\eta")

#=
### Back to the actual velocity
But what does this mean for the velocity? First, let's assemble the velocity,
now that we have $F$

$$ u(y,t) = U (1 - \mathrm{erf}(y/\sqrt{4\nu t})) $$
=#
u(y,t;U=1,ν=1) = U*F(y/δ(t,ν=ν))

#=
The velocity is time dependent, so we need to plot this at different times.
In fact, let's make it into a movie, so we can see how the profile changes:
=#
y = range(0,5,length=501)
timerange = range(0.001,2,length=81)

@gif for t in timerange
    plot(u.(y,t),y,xlim=(0,2),ylim=(0,Inf),xlabel=L"u/U",ylabel=L"y",label="Profile",title=string("t = ",round(t,digits=2)))
    hline!([δ(t)],linecolor=:black,linewidth=1,label=L"\delta")
    vline!([F(1)],label="",linecolor=:black,linewidth=0.5,linestyle=:dash,legend=:true)
end

#=
This movie shows the layer getting thicker with time. However, the
shape of the velocity profile is always the same -- it just gets stretched out
as the layer grows!

Another way to understand this. Let's evaluate velocity at $y = 3$ at $t = 1$.
It is very small, about 3 percent of the wall velocity ($U=1$)
=#
u(3,1)

#=
Now let's evaluate it at $t = 10$. It's now about half of the wall speed ($U=1$):
=#
u(3,10)

#=
### Vorticity
We hinted before that something funny happens at $t = 0$. The velocity
is equal to $U$ at $y = 0$ and $0$ right above it. This is a discontinuous
velocity, which viscosity immediately smooths over.

The discontinuous velocity means an *infinite* vorticity at $y = 0$ and 0 above it.
As time proceeds, this spike in vorticity (a Dirac delta function) gets
smoothed out into a Gaussian (bell curve) shape, maximum at the wall. In
other words, vorticity is *generated* at $t = 0$ at the wall, and then
spread upward by diffusion.
=#
Fp(η) = 2/sqrt(π)*exp(-η^2)
ω(y,t;U=1,ν=1) = U/δ(t,ν=ν)*Fp(y/δ(t,ν=ν))

#-
@gif for t in timerange
    plot(ω.(y,t),y,xlim=(0,2),ylim=(0,Inf),xlabel=L"\omega",ylabel=L"y",label="Profile",title=string("t = ",round(t,digits=2)))
    hline!([δ(t)],linecolor=:black,linewidth=1,label=L"\delta")
end

#=
Though the vorticity spreads out and gets weaker, the integral under this curve
does not actually change at all! This is the *circulation* (per unit length of
the wall). The circulation remains constant over time.
=#