text:
#####
##### 3D modes on pyramids
#####

"""
    basis(elem::Pyr,N,r,s,t,tol=1e-12)

Computes orthonormal semi-nodal basis on the biunit pyramid element.

Warning: nodal derivative matrices may contain errors for nodes at t = 1.
A way to avoid this is to use weak differentiation matrices computed using
quadrature rules with only interior nodes.
"""
function basis(elem::Pyr, N, r, s, t, tol = 1e3*eps())

    # convert to abc
    a = @. 2*(r+1)/(1-t)-1
    b = @. 2*(s+1)/(1-t)-1
    c = t

    # fix top point
    ids = @. abs(1 - t) < tol
    t[ids] .= 1 - tol # WARNING: hacky (note factor of (1-c) = (1-t) in "pc" below)
    a[ids] .= -1
    b[ids] .= -1

    # change of vars from a to b
    dadr = @. 2 / (1-t)
    dbds = @. 2 / (1-t)
    dadt = @. (2 * r + 2) / (1 - t)^2
    dbdt = @. (2 * s + 2) / (1 - t)^2
    dcdt = 1

    Np = (N + 1) * (N + 2) * (2 * N + 3) ÷ 6
    V, Vr, Vs, Vt = ntuple(x->zeros(length(r), Np), 4)

    ind = 1
    for k = 0:N
        # make nodal quad basis
        bq, aq, wab = quad_nodes(Quad(), k)
        VDM, _ = basis(Quad(), k, aq, bq)
        VDMab, Va, Vb = basis(Quad(), k, a, b)
        Vab, DVa, DVb = map(A->A / VDM, (VDMab, Va, Vb))

        CNk = (N+2) / (2^(2*k+2) * (2*k+3))
        pNk = jacobiP(c, 2*k+3, 0, N-k)
        pc = @. ((1-c) / 2)^k * pNk
        dpNk = grad_jacobiP(c, 2*k+3, 0, N-k)
        dpc = @. ((1-c)/2)^k * dpNk - (k/2^k) * ((1 - c)^(k - 1)) * pNk

        for ij = 1:(k+1)^2
            scale = 1 / sqrt(CNk * wab[ij]) # normalization
            @. V[:,ind]  = scale * Vab[:,ij] * pc
            @. Vr[:,ind] = scale * (DVa[:,ij] * dadr) * pc
            @. Vs[:,ind] = scale * (DVb[:,ij] * dbds) * pc
            @. Vt[:,ind] = scale * ((DVa[:,ij] * dadt + DVb[:,ij] * dbdt) * pc + Vab[:,ij] * dpc * dcdt)
            ind += 1
        end
    end
    return V, Vr, Vs, Vt
end

"""
    abctorst(elem::Pyr,a,b,c)

Converts from Stroud coordinates (a,b,c) on [-1,1]^3 to reference element
coordinates (r,s,t).
"""
function abctorst(elem::Pyr,a,b,c)
    r = @. .5 * (1+a) * (1-c) - 1
    s = @. .5 * (1+b) * (1-c) - 1
    t = @. c
    return r, s, t
end

"""
    nodes(elem::Pyr,N)

Computes interpolation nodes of degree N. Edge nodes coincide with (N+1)-point Lobatto
points. Triangular face nodes coincide with Tri.nodes(N), quadrilateral face nodes
coincide with tensor product (N+1)-point Lobatto points.
"""
function nodes(elem::Pyr, N)

    if N == 1
        return equi_nodes(Pyr(), N)
    end

    # append quad face nodes
    r1D_equi = equi_nodes(Line(), N)
    rst_equi = interp_1D_to_edges(Pyr(), r1D_equi)
    quad_face_pts_equi = (vec.(meshgrid(r1D_equi[2:end-1]))..., -ones((N - 1) * (N - 1)))
    append!.(rst_equi, quad_face_pts_equi)

    # append quad face nodes
    r1D = nodes(Line(), N)
    rst_lobatto = interp_1D_to_edges(Pyr(), r1D)
    quad_face_pts = (vec.(meshgrid(r1D[2:end-1]))..., -ones((N - 1) * (N - 1)))
    append!.(rst_lobatto, quad_face_pts)

    # append quad face "bubble" basis functions
    function edge_face_pyr_basis(N, rst)
        V_edge = edge_basis(Pyr(), N, rst...)
        # WARNING: assumes first 4 vertices of Pyr = quad vertices
        # also assumes first 5 cols of V_edge = vertex functions for Pyr
        V_quad_face = zeros(length(first(rst)), (N - 1) * (N - 1))
        V_quad = vandermonde(Quad(), N-2, rst[1], rst[2])
        V_quad_bubble = @. V_edge[:,1] * V_edge[:,2] * V_edge[:,3] * V_edge[:,4]
        for i = 1:size(V_quad_face, 2)
            @. V_quad_face[:,i] = V_quad_bubble * V_quad[:,i]
        end
        return hcat(V_edge, V_quad_face)
    end

    # same interp problem as before
    V_edge_face = edge_face_pyr_basis(N, rst_equi)
    c = (x->V_edge_face \ x).(rst_lobatto) # should be lsq solve
    return (x->edge_face_pyr_basis(N, equi_nodes(elem, N)) * x).(c)
end

function equi_nodes(elem::Pyr, N)
    # Np = (N+1)*(N+2)*(2*N+3)/6
    a, b, c = ntuple(x->Float64[],3)
    c1D = LinRange(-1, 1, N+1)
    for k = N:-1:0
        if k==0
            r1D = [.5]
        else
            r1D = LinRange(-1, 1, k + 1)
        end
        ak, bk = vec.(meshgrid(r1D, r1D))
        ck = c1D[N + 1 - k]*one.(ak)
        append!.((a, b, c),(ak, bk, ck))
    end
    return abctorst(elem, a, b, c)
end

quad_nodes(elem::Pyr, N) = stroud_quad_nodes(elem, N)

function stroud_quad_nodes(elem::Pyr, N)

    a1D, w1D = gauss_quad(0, 0, N)
    c1D, wc1D = gauss_quad(2, 0, N)

    a, c, b = vec.(meshgrid(a1D, c1D, a1D))
    wa, wc, wb = vec.(meshgrid(w1D, wc1D, w1D))
    w = @. wa * wb * wc
    w = (8/3) * w ./ sum(w) # scale by b*h/3 = volume of pyr. b = 4, h = 2 for biunit

    return abctorst(elem, a, b, c)..., w
end


text:
using Random

const rng = MersenneTwister()
const letters = 'A':'Z'
const digits = '0':'9'

const possiblenames = [c1 * c2 * d1 * d2 * d3 for c1 in letters for c2 in letters for d1 in digits for d2 in digits for d3 in digits]
shuffle!(rng, possiblenames)

mutable struct Robot
    name::String
    Robot() = new(generatename!())
end

function reset!(instance::Robot)
    instance.name = generatename!()
    return instance
end

function generatename!()
    isempty(possiblenames) ? error("there are no more available robot names") : pop!(possiblenames)
end

function name(instance::Robot)
    instance.name
end

ext: jl

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
module Tests

include("utilities.jl")

include("term_test.jl")
include("bivariate_quadratic_test.jl")
include("concave_then_convex_cubic_test.jl")
include("horn_test.jl")

const linear_tests = Dict(
    "dsos_term" => dsos_term_test,
    "dsos_bivariate_quadratic" => dsos_bivariate_quadratic_test,
    "dsos_concave_then_convex_cubic" => dsos_concave_then_convex_cubic_test,
    "dsos_horn_test" => dsos_horn_test
)

const soc_tests = Dict(
    "sdsos_term" => sdsos_term_test,
    "sdsos_bivariate_quadratic" => sdsos_bivariate_quadratic_test,
    "sdsos_concave_then_convex_cubic" => sdsos_concave_then_convex_cubic_test,
    "sdsos_horn_test" => sdsos_horn_test
)

const sd_tests = Dict(
    "sos_term" => sos_term_test,
    "sos_bivariate_quadratic" => sos_bivariate_quadratic_test,
    "sos_concave_then_convex_cubic" => sos_concave_then_convex_cubic_test,
    "sos_horn_test" => sos_horn_test
)

@test_suite linear
@test_suite soc
@test_suite sd

end


text:
using Caesar, RoME, IncrementalInference
using RoMEPlotting, Gadfly, KernelDensityEstimatePlotting
using KernelDensityEstimate
using DocStringExtensions
using DelimitedFiles, JLD

include(joinpath(@__DIR__,"slamUtils.jl"))
include(joinpath(@__DIR__,"plotSASUtils.jl"))

expnum = 2;
if expnum == 1
    topdir = joinpath(ENV["HOME"],"data","kayaks","07_18")
    trialstr = "07_18_parsed_set1";
    datadir = joinpath(topdir,trialstr);
    ijldname = datadir * "/exp1" * ".jld"
elseif expnum==2
    topdir = joinpath(ENV["HOME"],"data","kayaks","07_18")
    trialstr = "07_18_parsed_set2";
    datadir = joinpath(topdir,trialstr);
    ijldname = datadir * "/exp1" * ".jld"
end

cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml");
chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");
cfgd=loadConfigFile(cfgFile)
cfgd["range_snr_floor"]=0.9
iload = load(ijldname)
nrange = iload["nrange"]
irange = iload["irange"];
itemp = readdlm(datadir*"/inav.csv",',',Float64,'\n', skipstart=1);

wstart = 190; wlen = 25;
dataframes = collect(nrange[wstart]:nrange[wstart]+wlen);

saswindow = 5;
sasstart = 1;
nsasfac = 5;
sasdataframes = [zeros(Int,saswindow) for _ in 1:nsasfac];
sasposes = [];
irangenew = irange[wstart] : 5 : irange[wstart+nsasfac]
igt = [itemp[irangenew,2] itemp[irangenew,4]];

poses = [Symbol("x$i") for i in 1:wlen]
posData = importdata_nav(dataframes,datadir=datadir);
navchecked, errorind = sanitycheck_nav(posData)

# dposData = deepcopy(posData)
# cumulativeDrift!(dposData,[0.0;0],[0.1,0.1])

# Check the trajectory
posPl = Gadfly.plot(layer(x=posData[:,1],y=posData[:,2], Geom.point, Geom.path, Theme(default_color=colorant"green",highlight_width = 0pt)),layer(x=igt[:,1],y=igt[:,2], Geom.point, Geom.path, Theme(default_color=colorant"red" ,highlight_width = 0pt))); posPl |> PDF("/tmp/test.pdf");

for i=1:nsasfac
    framestart = nrange[wstart]+sasstart-1+(i-1)*saswindow;
    sasdataframes[i] = collect(framestart:framestart+saswindow-1)
    push!(sasposes,[Symbol("x$j") for j in sasstart+(i-1)*saswindow:sasstart+saswindow-1+(i-1)*saswindow])
end

fg = initfg();
rtkCov= Matrix(Diagonal([0.1;0.1].^2));
for i in 1:wlen
  addVariable!(fg, poses[i], Point2)
  pp = PriorPoint2(MvNormal(posData[i,:], rtkCov))
  addFactor!(fg, [poses[i];], pp, autoinit=true)
end

beacons = [Symbol("l$j") for j in 1:nsasfac]

for i=1:nsasfac
    addVariable!(fg, beacons[i], DynPoint2(ut=1_000_000+(i-1)*saswindow) )

    if i>1
        dpμ = [0.0;0;0;0];
        dpσ = Matrix(Diagonal([1.0;1;0.1;0.1].^2))
        vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
        addFactor!(fg, [beacons[i-1];beacons[i]], vp, autoinit=false)
    end

    waveformData = importdata_waveforms(sasdataframes[i],2, datadir=datadir);
    sas2d = prepareSAS2DFactor(saswindow, waveformData, rangemodel=:Correlator,
                               cfgd=cfgd, chirpFile=chirpFile)
    addFactor!(fg, [beacons[i];sasposes[i]], sas2d, autoinit=false)
    @show currfac = ls(fg,beacons[i])[1]

    tmpInit = approxConv(fg,currfac,beacons[i]);
    XXkde = manikde!(tmpInit,getManifolds(fg,beacons[i]));
    setValKDE!(fg,beacons[i],XXkde);
end

# visualization tools for debugging
# writeGraphPdf(fg,viewerapp="", engine="neato", filepath = "/tmp/test.pdf")

# getSolverParams(fg).drawtree = true
getSolverParams(fg).limititers=500
#getSolverParams(fg).showtree = true

## solve the factor graph
tree, smt, hist = solveTree!(fg)

# Some Debuggging
# @time L1 = approxConv(fg,:l1x1x2x3x4x5x6x7x8f1,:l1)
# @time L1 = approxConv(fg,:l2x9x10x11x12x13x14x15x16f1,:l2)

# plk = plotBeaconSolve(fg)
plk=[];
for sym in ls(fg) #plotting all syms labeled
    X1 = getKDEMean(getVertKDE(fg,sym))
    push!(plk, layer(x=[X1[1];],y=[X1[2];],label=String["$sym";], Geom.point,Geom.label), Theme(default_color=colorant"red",point_size = 1.5pt,highlight_width = 0pt))
end
# push!(plk,layer(x=L1[1,:],y=L1[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))

K1 = plotKDEContour(kde!(getVal(fg,:l1)[1:2,:]),xlbl="X (m)", ylbl="Y (m)",levels=5,layers=true);
push!(plk,K1...)
push!(plk,Gadfly.Theme(key_position = :none));

push!(plk,layer(x=posData[:,1],y=posData[:,2], Geom.point, Geom.path, Theme(default_color=colorant"blue")), layer(x=igt[:,1],y=igt[:,2], Geom.point, Geom.path, Theme(default_color=colorant"green" )))
push!(plk, Coord.cartesian(xmin=0, xmax=100, ymin=-50, ymax=20,fixed=true))
plkplot = Gadfly.plot(plk...); plkplot |> PDF("/tmp/test.pdf")

plk= [];

push!(plk,layer(x=posData[:,1],y=posData[:,2], Geom.path, Theme(default_color=colorant"green")), layer(x=igt[:,1],y=igt[:,2], Geom.point, Geom.path, Theme(default_color=colorant"red" )))

L1 = getVal(getVariable(fg, :l1))
push!(plk,layer(x=L1[1,:],y=L1[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))

# push!(plk, Coord.cartesian(xmin=40, xmax=100, ymin=-120, ymax=-50,fixed=true))
plkplot = Gadfly.plot(plk...); plkplot |> PDF("/tmp/test.pdf");

## Plot the SAS factor pairs
fsym = :l1x1x2x3x4x5x6x7x8f1
pl = plotSASPair(fg, fsym, show=true, filepath="/tmp/testDP_2pp_3vpvp_init.pdf");


text:
# setup model
# struct to hold values
struct Inputs
    parameters::Dict{Symbol,Float64}
    initial_conditions::Dict{Symbol,Float64}
    torque_parameters::Dict{Symbol,Dict{Symbol,Float64}}
    activation_parameters::Dict{Symbol,Vector{Float64}}
    matching_data::Dict{Symbol,Vector{Float64}}
end

function load_inputs(;
    parameters="data/parameters.csv",
    initial_conditions="data/initial_conditions.csv",
    torque_generators="data/torque_generator_parameters.csv",
    activations="data/activation_parameters.csv",
    matching="data/matchingData.csv")

    # load parameters
    input_p, headers_p = readdlm(parameters, ',', header=true)
    input_p = Vector{Float64}(input_p[2, :]) # remove units row
    p = Dict(Symbol(h) => v for (h, v) in zip(headers_p, input_p))

    # load initial conditions
    input_u, headers_u = readdlm(initial_conditions, ',', header=true)
    input_u = Vector{Float64}(input_u[2, :]) # remove units row
    u = Dict(Symbol(h) => v for (h, v) in zip(headers_u, input_u))

    # load torque generator parameters
    input_tq, headers_tq = readdlm(torque_generators, ',', header=true)
    tq = Dict(
        Symbol(x[1]) => Dict(
            Symbol(h) => Float64(v) for (h, v) in zip(headers_tq[2:end], x[2:end])
        ) for x in eachrow(input_tq)
    )

    # load activation profiles
    input_act, headers_act = readdlm(activations, ',', header=true)
    α = Dict(Symbol(x[1]) => Vector{Float64}(x[2:end]) for x in eachrow(input_act))

    # load matching data
    input_matching, headers_matching = readdlm(matching, ',', header=true)
    input_matching = Matrix{Float64}(input_matching[2:end, :]) # remove units row
    matching_data = Dict(Symbol(h) => ts for (h, ts) in zip(headers_matching, eachcol(input_matching)))

    # return in struct
    return Inputs(p, u, tq, α, matching_data)
end

# returns the parameters for a given speed from a csv file
function read_from_results(fname, speed)
    p, header = readdlm(fname, ',', Float64, header=true)
    row = p[p[:, 1].==speed, :]

    return row[3:end]
end

# returns u₀ and p created from results.csv
function load_from_results(inputs, fname, speed)
    p = read_from_results(fname, speed)

    # create new inputs
    _inputs = deepcopy(inputs)
    _inputs.initial_conditions[:vcmx] = speed
    # _inputs.activation_parameters[:, 2:end] = reshape(p, :, 6) |> permutedims # reshape to 6xn matrix, row major
    n = length(p)
    m = n ÷ 6
    [_inputs.activation_parameters[k] = p[i:i+m-1] for (k, i) in zip([:he, :ke, :ae, :hf, :kf, :af], collect(1:m:n))]

    p, u₀ = set_values(_inputs)
    return p, u₀
end

function set_values(inputs)
    @unpack parameters, initial_conditions, torque_parameters, activation_parameters = inputs
    parameters = NamedTuple(parameters)
    initial_conditions = NamedTuple(initial_conditions)

    # set parameters
    @inbounds (; footang, g, ina, inb, inc, ind, ine, inf, ing, inh, ini, k1, k2, k3, k4, k5, k6, k7, k8, l1, l10, l11, l2, l3, l4, l5, l6, l7, l8, l9, ma, mb, mc, md, me, mf, mg, mh, mi, mtpb, mtpk, pop1xi, pop2xi) = parameters
    mt = ma + mb + mc + md + me + mf + mg
    u8 = u9 = u10 = u11 = 0.0
    footang = deg2rad(footang)

    # set intial conditions
    @inbounds (; mtp_x, mtp_y, q3, mtp_angle, ankle_angle, knee_angle, hip_angle, vcmx, vcmy, u3, mtp_angular_velocity, ankle_angular_velocity, knee_angular_velocity, hip_angular_velocity) = initial_conditions

    # set torque parameters
    tq_p = Dict(k => [CCParameters(v), SECParameters(v)] for (k, v) in torque_parameters)

    # load activation profiles
    α_p = Dict(k => ActivationProfile(v) for (k, v) in activation_parameters)

    # convert initial conditions into generalised coordinates and speeds
    q3 = q3 |> deg2rad              # q3 = q3
    q1 = mtp_x + l2 * cos(q3 - q4 - q5 - q6 - q7)   # pop2x = q1 - l2*cos(q3-q4-q5-q6-q7)  (MTP contacts first)
    q2 = mtp_y + l2 * sin(q3 - q4 - q5 - q6 - q7)   # pop2y = q2 - l2*sin(q3-q4-q5-q6-q7)
    u3 = u3 |> deg2rad              # u3 = u3

    # initial angle driven values
    ra₀ = _ra(0.0)
    la₀ = _la(0.0)
    rk₀ = _rk(0.0)
    lk₀ = _lk(0.0)
    rh₀ = _rh(0.0)
    lh₀ = _lh(0.0)
    rap₀ = _rap(0.0)
    lap₀ = _lap(0.0)
    rkp₀ = _rkp(0.0)
    lkp₀ = _lkp(0.0)
    rhp₀ = _rhp(0.0)
    lhp₀ = _lhp(0.0)
    hato₀ = _hato(0.0)
    hatop₀ = _hatop(0.0)

    # u1 = 1.78872
    # u2 = -3.877485
    u1_fun(ra, la, rk, lk, rh, lh, hato, rap, lap, rkp, lkp, rhp, lhp, hatop) = vcmx + 0.5 * (2 * mhat * hato * sin(q3) * u3 + 2 * (lth * mff + lth * mrf + lth * msh + mth * (lth - ltho)) * sin(lh + q3) * (lhp + u3 + u5) + 2 * (lsh * mff + lsh * mrf + msh * (lsh - lsho)) * sin(lh + lk + q3) * (lhp + lkp + u3 + u5 + u7) + 2 * mff * (lff - lffo) * sin(la + lh + lk + lmtp + q3) * (lap + lhp + lkp + lmtpp + u3 + u5 + u7 + u9) + (2 * lrf * mff + mrf * (2 * lrf - lrfo)) * sin(la + lh + lk + q3) * (lap + lhp + lkp + u3 + u5 + u7 + u9) - 2 * mhat * hatop * cos(q3) - lrffo * mrf * sin(footang + la + lh + lk + q3) * (lap + lhp + lkp + u3 + u5 + u7 + u9) - lrffo * mrf * sin(footang + ra + rh + rk + q3) * (rap + rhp + rkp + u3 + u4 + u6 + u8) - 2 * (lth * mff + lth * mhat + lth * mrf + lth * msh + lth * mth + ltho * mth) * sin(rh + q3) * (rhp + u3 + u4) - 2 * (lsh * mff + lsh * mhat + lsh * mrf + lsh * msh + lsho * msh + 2 * lsh * mth) * sin(rh + rk + q3) * (rhp + rkp + u3 + u4 + u6) - 2 * (lff * mff + lff * mhat + lffo * mff + 2 * lff * mrf + 2 * lff * msh + 2 * lff * mth) * sin(ra + rh + rk + rmtp + q3) * (rap + rhp + rkp + rmtpp + u3 + u4 + u6 + u8) - (lrfo * mrf + 2 * lrf * mff + 2 * lrf * mhat + 2 * lrf * mrf + 4 * lrf * msh + 4 * lrf * mth) * sin(ra + rh + rk + q3) * (rap + rhp + rkp + u3 + u4 + u6 + u8)) / (mhat + 2 * mff + 2 * mrf + 2 * msh + 2 * mth)
    u2_fun(ra, la, rk, lk, rh, lh, hato, rap, lap, rkp, lkp, rhp, lhp, hatop) = vcmy - 0.5 * (2 * mhat * hatop * sin(q3) + 2 * mhat * hato * cos(q3) * u3 + 2 * (lth * mff + lth * mrf + lth * msh + mth * (lth - ltho)) * cos(lh + q3) * (lhp + u3 + u5) + 2 * (lsh * mff + lsh * mrf + msh * (lsh - lsho)) * cos(lh + lk + q3) * (lhp + lkp + u3 + u5 + u7) + 2 * mff * (lff - lffo) * cos(la + lh + lk + lmtp + q3) * (lap + lhp + lkp + lmtpp + u3 + u5 + u7 + u9) + (2 * lrf * mff + mrf * (2 * lrf - lrfo)) * cos(la + lh + lk + q3) * (lap + lhp + lkp + u3 + u5 + u7 + u9) - lrffo * mrf * cos(footang + la + lh + lk + q3) * (lap + lhp + lkp + u3 + u5 + u7 + u9) - lrffo * mrf * cos(footang + ra + rh + rk + q3) * (rap + rhp + rkp + u3 + u4 + u6 + u8) - 2 * (lth * mff + lth * mhat + lth * mrf + lth * msh + lth * mth + ltho * mth) * cos(rh + q3) * (rhp + u3 + u4) - 2 * (lsh * mff + lsh * mhat + lsh * mrf + lsh * msh + lsho * msh + 2 * lsh * mth) * cos(rh + rk + q3) * (rhp + rkp + u3 + u4 + u6) - 2 * (lff * mff + lff * mhat + lffo * mff + 2 * lff * mrf + 2 * lff * msh + 2 * lff * mth) * cos(ra + rh + rk + rmtp + q3) * (rap + rhp + rkp + rmtpp + u3 + u4 + u6 + u8) - (lrfo * mrf + 2 * lrf * mff + 2 * lrf * mhat + 2 * lrf * mrf + 4 * lrf * msh + 4 * lrf * mth) * cos(ra + rh + rk + q3) * (rap + rhp + rkp + u3 + u4 + u6 + u8)) / (mhat + 2 * mff + 2 * mrf + 2 * msh + 2 * mth)
    u1 = u1_fun(ra₀, la₀, rk₀, lk₀, rh₀, lh₀, hato₀, rap₀, lap₀, rkp₀, lkp₀, rhp₀, lhp₀, hatop₀)
    u2 = u2_fun(ra₀, la₀, rk₀, lk₀, rh₀, lh₀, hato₀, rap₀, lap₀, rkp₀, lkp₀, rhp₀, lhp₀, hatop₀)

    # parameters
    p = Params(footang, g, iff, ihat, irf, ish, ith, k1, k2, k3, k4, k5, k6, k7, k8, lff, lffo, lhat, lmtpxi, lrf, lrff, lrffo, lrfo, lsh, lsho, lth, ltho, ltoexi, mff, mhat, mrf, msh, mth, mtpb, mtpk, rmtpxi, rtoexi)

    # intial conditions
    u₀ = SVector(q1, q2, q3, u1, u2, u3)

    return p, u₀

end

function setup(;
    parameters="model/angle_driven/parameters.csv",
    initial_conditions="model/angle_driven/initial_conditions.csv",
    matching_data="data/matchingData.csv")

    # load parameters
    inputs = readdlm(parameters, ',', header=true)
    input_p, headers_p = convert(Vector{Float64}, inputs[1][:, 3]), Symbol.(inputs[1][:, 1])

    footang, g, iff, irf, ish, ith, ihat, k1, k2, k3, k4, k5, k6, k7, k8, lffo, lff, lrffo, lrfo, lrff, lrf, lsho, lsh, ltho, lth, lhat, mff, mrf, msh, mth, mhat, mtpb, mtpk, rtoexi, rmtpxi, ltoexi, lmtpxi = input_p
    u4 = u5 = u6 = u6 = u7 = u8 = u9 = 0.0
    footang = deg2rad(footang)

    # load initial conditions
    inputs = readdlm(initial_conditions, ',', header=true)
    input_u, header_u = convert(Vector{Float64}, inputs[1][:, 3]), Symbol.(inputs[1][:, 1])
    @inbounds Xmtp, Ymtp, q3, vcmx, vcmy, u3 = input_u

    # load matching data
    matching_data, headers_matching = readdlm(matching_data, ',', Float64, skipstart=1, header=true)
    matching_data = matching_data[1:2:end-2] # only use left side data

    # intial specified angles
    rmtp₀ = _rmtp(0.0)
    lmtp₀ = _lmtp(0.0)
    ra₀ = _ra(0.0)
    la₀ = _la(0.0)
    rk₀ = _rk(0.0)
    lk₀ = _lk(0.0)
    rh₀ = _rh(0.0)
    lh₀ = _lh(0.0)
    rmtpp₀ = _rmtpp(0.0)
    lmtpp₀ = _lmtpp(0.0)
    rap₀ = _rap(0.0)
    lap₀ = _lap(0.0)
    rkp₀ = _rkp(0.0)
    lkp₀ = _lkp(0.0)
    rhp₀ = _rhp(0.0)
    lhp₀ = _lhp(0.0)
    hato₀ = _hato(0.0)
    hatop₀ = _hatop(0.0)

    # convert initial conditions into generalised coordinates and speeds
    q3 = q3 |> deg2rad              # q3 = q3
    q1 = Xmtp - lff * cos(ra₀ + rh₀ + lmtp₀ - lk₀ - q3)   # POP2X = Q1 + LFF*COS(LA+LH+LMTP-LK-Q3)  (MTP contacts first)
    q2 = Ymtp + lff * sin(la₀ + lh₀ + lmtp₀ - lk₀ - q3)   # POP2Y = Q2 - LFF*SIN(LA+LH+LMTP-LK-Q3)
    u3 = u3 |> deg2rad              # u3 = u3

    # u1 = 1.78872
    # u2 = -3.877485
    u1_fun(rmtp, lmtp, ra, la, rk, lk, rh, lh, hato, rmtpp, lmtpp, rap, lap, rkp, lkp, rhp, lhp, hatop) = vcmx + 0.5 * (2 * mhat * hato * sin(q3) * u3 + 2 * (lth * mff + lth * mrf + lth * msh + mth * (lth - ltho)) * sin(rh - q3) * (rhp - u3 - u4) + (2 * lrf * mff + mrf * (2 * lrf - lrfo)) * sin(ra + rh - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) + 2 * (lsh * mff + lsh * mhat + lsh * mrf + lsh * msh + lsho * msh + 2 * lsh * mth) * sin(lh - lk - q3) * (lhp - lkp - u3 - u5 - u7) + 2 * (lff * mff + lff * mhat + lffo * mff + 2 * lff * mrf + 2 * lff * msh + 2 * lff * mth) * sin(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9) - 2 * mhat * hatop * cos(q3) - 2 * (lth * mff + lth * mhat + lth * mrf + lth * msh + lth * mth + ltho * mth) * sin(lh - q3) * (lhp - u3 - u5) - lrffo * mrf * sin(la + lh - footang - lk - q3) * (lap + lhp - lkp - u3 - u5 - u7 - u9) - lrffo * mrf * sin(ra + rh - footang - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) - 2 * (lsh * mff + lsh * mrf + msh * (lsh - lsho)) * sin(rh - rk - q3) * (rhp - rkp - u3 - u4 - u6) - 2 * mff * (lff - lffo) * sin(ra + rh + rmtp - rk - q3) * (rap + rhp + rmtpp - rkp - u3 - u4 - u6 - u8) - (lrfo * mrf + 2 * lrf * mff + 2 * lrf * mhat + 2 * lrf * mrf + 4 * lrf * msh + 4 * lrf * mth) * sin(la + lh - lk - q3) * (lap + lhp - lkp - u3 - u5 - u7 - u9)) / (mhat + 2 * mff + 2 * mrf + 2 * msh + 2 * mth)
    u2_fun(rmtp, lmtp, ra, la, rk, lk, rh, lh, hato, rmtpp, lmtpp, rap, lap, rkp, lkp, rhp, lhp, hatop) = vcmy + 0.5 * (2 * (lth * mff + lth * mrf + lth * msh + mth * (lth - ltho)) * cos(rh - q3) * (rhp - u3 - u4) + (2 * lrf * mff + mrf * (2 * lrf - lrfo)) * cos(ra + rh - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) + 2 * (lsh * mff + lsh * mhat + lsh * mrf + lsh * msh + lsho * msh + 2 * lsh * mth) * cos(lh - lk - q3) * (lhp - lkp - u3 - u5 - u7) + 2 * (lff * mff + lff * mhat + lffo * mff + 2 * lff * mrf + 2 * lff * msh + 2 * lff * mth) * cos(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9) - 2 * mhat * hatop * sin(q3) - 2 * mhat * hato * cos(q3) * u3 - 2 * (lth * mff + lth * mhat + lth * mrf + lth * msh + lth * mth + ltho * mth) * cos(lh - q3) * (lhp - u3 - u5) - lrffo * mrf * cos(la + lh - footang - lk - q3) * (lap + lhp - lkp - u3 - u5 - u7 - u9) - lrffo * mrf * cos(ra + rh - footang - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) - 2 * (lsh * mff + lsh * mrf + msh * (lsh - lsho)) * cos(rh - rk - q3) * (rhp - rkp - u3 - u4 - u6) - 2 * mff * (lff - lffo) * cos(ra + rh + rmtp - rk - q3) * (rap + rhp + rmtpp - rkp - u3 - u4 - u6 - u8) - (lrfo * mrf + 2 * lrf * mff + 2 * lrf * mhat + 2 * lrf * mrf + 4 * lrf * msh + 4 * lrf * mth) * cos(la + lh - lk - q3) * (lap + lhp - lkp - u3 - u5 - u7 - u9)) / (mhat + 2 * mff + 2 * mrf + 2 * msh + 2 * mth)
    u1 = u1_fun(rmtp₀, lmtp₀, ra₀, la₀, rk₀, lk₀, rh₀, lh₀, hato₀, rmtpp₀, lmtpp₀, rap₀, lap₀, rkp₀, lkp₀, rhp₀, lhp₀, hatop₀)
    u2 = u2_fun(rmtp₀, lmtp₀, ra₀, la₀, rk₀, lk₀, rh₀, lh₀, hato₀, rmtpp₀, lmtpp₀, rap₀, lap₀, rkp₀, lkp₀, rhp₀, lhp₀, hatop₀)

    # parameters
    p = Params(footang, g, iff, ihat, irf, ish, ith, k1, k2, k3, k4, k5, k6, k7, k8, lff, lffo, lhat, lmtpxi, lrf, lrff, lrffo, lrfo, lsh, lsho, lth, ltho, ltoexi, mff, mhat, mrf, msh, mth, mtpb, mtpk, rmtpxi, rtoexi)

    # intial conditions
    u₀ = SVector(q1, q2, q3, u1, u2, u3)

    return p, u₀, matching_data
end


text:
#!/usr/bin/env julia
using GraphIdx.Io.Snap: open_snap, parse_snap_edges
using HDF5: h5open, attrs

fname = get(ARGS, 1, "com-youtube.txt.gz")
out = get(ARGS, 2, "out.h5")
output_bin = false

@time n, m, head, tail = open_snap(parse_snap_edges, fname)

isfile(out) && rm(out)
h5open(out, "w") do io
    attrs(io)["n"] = n
    attrs(io)["m"] = m
    write(io, "head", head .-1)
    write(io, "tail", tail .-1)
end

output_bin && open("bin", "w") do io
    write(io, head .-1);
    write(io, tail .-1)
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
const regex = r"<x=(-?\d+), y=(-?\d+), z=(-?\d+)>";
const input = readlines(length(ARGS) > 0 ? ARGS[1] : "example1.txt");
pos = hcat(map(line -> parse.(Int,match(regex, line).captures), input)...);
const start = copy(pos)

function gravity(p,i,d)
    dv = 0
    for j in p[d,:]
        if p[d,i] > j
            dv -= 1
        elseif p[d,i] < j
            dv += 1
        end
    end
    return dv
end

function energy(positions, velocities, moon)
    potential = sum(abs.(positions[:,moon]))
    kinetic = sum(abs.(velocities[:,moon]))
    return (potential, kinetic, potential * kinetic)
end

vel = zeros(Int,3,4);
const period_length = zeros(Int,3)

const show_output = false

for i in 1:500000
    dv = [gravity(pos,body,dimension) for dimension in 1:3, body in 1:size(pos)[2]]
    global vel += dv
    global pos += vel
    if (show_output)
        println("After $(i) step$(i != 1 ? "s" : ""):")
        println("Positions:")
        display(pos')
        println()
        println("Velocities:")
        display(vel')
        println()
    end

    for j in 1:3
        if period_length[j] == 0 && pos[j,:] == start[j,:]
            period_length[j] = i + 1
        end
    end

    if i == 1000
        println("Day 12 Part 1: $(sum([energy(pos,vel,i)[3] for i in 1:4]))")
    end

    if i > 1000 && !(0 in period_length)
        break;
    end
end

println("Day 12 Part 2: $(foldl(lcm, period_length))")


text:
using FluxArchitectures
using Test

if Flux.CUDA.functional()
    @info "Testing GPU support"
else
    @warn "CUDA unavailable, not testing GPU support"
end

@testset "DARNN" begin
    include("DARNN.jl")
end

@testset "DSANet" begin
    include("DSANet.jl")
end

@testset "LSTnet" begin
    include("LSTnet.jl")
end

@testset "TPALSTM" begin
    include("TPALSTM.jl")
end

@testset "dataloader" begin
    include("dataloader.jl")
end

@testset "shared" begin
    include("shared.jl")
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
################################################################################
#
#   DIAMOND LATTICE
#
################################################################################

# Implementation
# - implementation 1 (only NN)
# - labels <: Any
function getUnitcellDiamond(
			unitcell_type  :: Type{U},
			implementation :: Val{1}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[0, 1, 1] .* 0.5,
			Float64[1, 0, 1] .* 0.5,
			Float64[1, 1, 0] .* 0.5,
		],
		# sites
		S[
			newSite(S, Float64[0.0, 0.0, 0.0], getDefaultLabelN(LS,1)),
			newSite(S, Float64[0.25, 0.25, 0.25], getDefaultLabelN(LS,2)),
		],
		# bonds
		B[
			newBond(B, 1, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 1, 2, getDefaultLabel(LB), (-1, 0, 0)),
			newBond(B, 1, 2, getDefaultLabel(LB), (0, -1, 0)),
			newBond(B, 1, 2, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 2, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 1, getDefaultLabel(LB), (1, 0, 0)),
			newBond(B, 2, 1, getDefaultLabel(LB), (0, 1, 0)),
			newBond(B, 2, 1, getDefaultLabel(LB), (0, 0, 1)),
		]
	)
end

# Implementation
# - implementation 2 (including NNN)
# - labels <: Any
function getUnitcellDiamond(
			unitcell_type  :: Type{U},
			implementation :: Val{2}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[0, 1, 1] .* 0.5,
			Float64[1, 0, 1] .* 0.5,
			Float64[1, 1, 0] .* 0.5,
		],
		# sites
		S[
			newSite(S, Float64[0.0, 0.0, 0.0], getDefaultLabelN(LS,1)),
			newSite(S, Float64[0.25, 0.25, 0.25], getDefaultLabelN(LS,2)),
		],
		# bonds
		B[
			newBond(B, 1, 2, getDefaultLabelN(LB, 1), (0, 0, 0)),
			newBond(B, 1, 2, getDefaultLabelN(LB, 1), (-1, 0, 0)),
			newBond(B, 1, 2, getDefaultLabelN(LB, 1), (0, -1, 0)),
			newBond(B, 1, 2, getDefaultLabelN(LB, 1), (0, 0, -1)),
			newBond(B, 2, 1, getDefaultLabelN(LB, 1), (0, 0, 0)),
			newBond(B, 2, 1, getDefaultLabelN(LB, 1), (1, 0, 0)),
			newBond(B, 2, 1, getDefaultLabelN(LB, 1), (0, 1, 0)),
			newBond(B, 2, 1, getDefaultLabelN(LB, 1), (0, 0, 1)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (1, 0, 0)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (-1, 0, 0)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (0, 1, 0)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (0, -1, 0)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (0, 0, 1)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (0, 0, -1)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (1, -1, 0)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (-1, 1, 0)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (0, 1, -1)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (0, -1, 1)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (-1, 0, 1)),
			newBond(B, 1, 1, getDefaultLabelN(LB, 2), (1, 0, -1)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (1, 0, 0)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (-1, 0, 0)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (0, 1, 0)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (0, -1, 0)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (0, 0, 1)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (0, 0, -1)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (1, -1, 0)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (-1, 1, 0)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (0, 1, -1)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (0, -1, 1)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (-1, 0, 1)),
			newBond(B, 2, 2, getDefaultLabelN(LB, 2), (1, 0, -1)),
		]
	)
end


text:
struct Lasso
  models::Vector{Tuple{Vector{Float64},Float64}}
end

function train(d::Data{Normalized}, ::Val{:lasso})
  M = length(d.y_T)

  models = map(1:M) do m
    X_m = d.X_T[m]
    y_m = d.y_T[m]
    n_m = length(y_m)

    if n_m < 10
      zeros(size(X_m, 2)), mean(y_m)
    else
      cv = GLMNet.glmnetcv(X_m, y_m)
      best = argmin(cv.meanloss)
      β = Vector(cv.path.betas[:, best])
      β0 = cv.path.a0[best]
      β, β0
    end
  end

  Lasso(models)
end

function getoutcomes(new_data::Data{Normalized}, lasso::Lasso)
  M = length(lasso.models)
  n = length(new_data.y)
  outcomes = Matrix{Float64}(undef, n, M)

  for m = 1:M
    β, β0 = lasso.models[m]
    outcomes[:, m] = new_data.X * β .+ β0
  end

  outcomes
end


text:
# Test that the settings adaptive_rho_max_adaptions is correctly implemented

using COSMO, SparseArrays, LinearAlgebra, Test

@testset "Max number of rho adaptions" begin
    q = [1; 1.];
    P = sparse([4. 1; 1 2]);
    A = [1. 1; 1 0; 0 1];
    l = [1.; 0; 0];
    u = [1; 0.7; 0.7];

    Aa = [-A; A]
    ba = [u; -l]
    constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);

    # Next, we define the settings object, the model and then assemble everything:
    settings = COSMO.Settings(verbose = false, adaptive_rho_interval = 25, adaptive_rho_max_adaptions = 2, rho = 1e-6, eps_abs = 1e-6, eps_rel = 1e-4);
    model = COSMO.Model();
    assemble!(model, P, q, constraint1, settings = settings);
    res = COSMO.optimize!(model);

    @test num_rho_adaptions(model) == 2

    # Solve again but restrict number of rho updates to 1
    settings = COSMO.Settings(verbose = false, adaptive_rho_interval = 25, adaptive_rho_max_adaptions = 1, rho = 1e-6, eps_abs = 1e-4, eps_rel = 1e-4);
    model = COSMO.Model();
    assemble!(model, P, q, constraint1, settings = settings);
    res = COSMO.optimize!(model);

    @test num_rho_adaptions(model) == 1
end


text:
using Base.Meta
export Arg, ∇, preprocess, @explicit_intercepts, @union_intercepts

"""
    @union_intercepts f type_tuple invoke_type_tuple [kwargs]

Interception strategy based on adding a method to `f` which accepts the union of each of
the types specified by `type_tuple`. If none of the arguments are `Node`s then the method
of `f` specified by `invoke_type_tuple` is invoked. If applicable, keyword arguments
should be provided as a `NamedTuple` and be added to the generated function's signature.
"""
macro union_intercepts(f::Symbol, type_tuple::Expr, invoke_type_tuple::Expr, kwargs::Expr=:(()))
    kwargs.head === :tuple || throw(ArgumentError("malformed keyword argument specification"))
    return esc(union_intercepts(f, type_tuple, invoke_type_tuple; parse_kwargs(kwargs)...))
end

"""
    union_intercepts(f::Symbol, type_tuple::Expr, invoke_type_tuple::Expr)

The work-horse for `@union_intercepts`.
"""
function union_intercepts(f::Symbol, type_tuple::Expr, invoke_type_tuple::Expr; kwargs...)
    call, arg_names = get_union_call(f, type_tuple; kwargs...)
    body = get_body(f, type_tuple, arg_names, invoke_type_tuple; kwargs...)
    return Expr(:macrocall, Symbol("@generated"), nothing, Expr(:function, call, body))
end

"""
    @explicit_intercepts(f::Symbol, type_tuple::Expr, is_node::Expr[, kwargs::Expr])
    @explicit_intercepts(f::Symbol, type_tuple::Expr)

Create a collection of methods which intecept the function calls to `f` in which at least
one argument is a `Node`. Types of arguments are specified by the type tuple expression
in `type_tuple`. If there are arguments which are not differentiable, they can be
specified by providing a boolean vector `is_node` which indicates those arguments that are
differentiable with `true` values and those which are not as `false`. Keyword arguments
to add to the function signature can be specified in `kwargs`, which must be a `NamedTuple`.
"""
macro explicit_intercepts(
    f::SymOrExpr,
    type_tuple::Expr,
    is_node::Expr=Expr(:vect, (true for _ in get_types(get_body(type_tuple)))...),
    kwargs::Expr=:(()),
)
    isnode = parse_is_node(is_node)
    ex = explicit_intercepts(f, type_tuple, isnode, __source__; parse_kwargs(kwargs)...)
    # The result contains all method definitions generated for f
    return esc(ex)
end

"""
    explicit_intercepts(f::Symbol, types::Expr, is_node::Vector)

Return a `:block` expression which evaluates to declare all of the combinations of methods
that could be required to catch if a `Node` is ever passed to the function specified in
`expr`.
"""
function explicit_intercepts(f::SymOrExpr, types::Expr, is_node::Vector{Bool}, linfo::LineNumberNode; kwargs...)
    function explicit_intercepts_(states::Vector{Bool})
        if length(states) == length(is_node)
            return any(states) ? boxed_method(f, types, states, linfo; kwargs...) : []
        else
            return vcat(
                explicit_intercepts_(vcat(states, false)),
                is_node[length(states) + 1] ? explicit_intercepts_(vcat(states, true)) : []
            )
        end
    end
    return Expr(:block, explicit_intercepts_(Vector{Bool}())...)
end

"""
    boxed_method(
        f::SymOrExpr,
        type_tuple::Expr,
        is_node::Vector{Bool},
        arg_names::Vector{Symbol};
        kwargs...
    )

Construct a method of the Function `f`, whose argument's types are specified by
`type_tuple`. Arguments which are potentially `Node`s should be indicated by `true` values
in `is_node`. Any provided keyword arguments will be added to the method.
"""
function boxed_method(
    f::SymOrExpr,
    type_tuple::Expr,
    is_node::Vector{Bool},
    arg_names::Vector{Symbol}=[gensym() for _ in is_node],
    linfo::LineNumberNode=LineNumberNode(0);
    kwargs...
)
    # Get the argument types and create the function call.
    types = get_types(get_body(type_tuple))
    noded_types = [node ? :(Node{<:$tp}) : tp for (node, tp) in zip(is_node, types)]
    call = replace_body(type_tuple, get_sig(f, arg_names, noded_types; kwargs...))

    # Construct body of call. If keyword arguments are absent, make a `Branch`, otherwise
    # call the internal function with positional arguments
    tuple_expr = Expr(:tuple, arg_names...)
    tape_expr = Expr(:call, :getfield, arg_names[findfirst(is_node)], quot(:tape))
    body = Expr(:call, :Branch, f, tuple_expr, tape_expr)

    # If keyword arguments are provided, add them to the generated call to `Branch`
    if !isempty(kwargs)
        params = Expr(:parameters)
        for (k, _) in kwargs
            push!(params.args, Expr(:kw, k, k))
        end
        insert!(body.args, 2, params)
    end

    # Combine call signature with the body to create a new function.
    return Expr(:function, call, Expr(:block, linfo, body))
end
boxed_method(f, t, n, l; kwargs...) = boxed_method(f, t, n, [gensym() for _ in n], l; kwargs...)

"""
    get_sig(f::SymOrExpr, arg_names::Vector{Symbol}, types::Vector; kwargs...)

Generate a function signature for `f` in which the arguments, whose names are `arg_names`,
specified by the `true` entires of `is_node` have type `Node`. The other arguments have
types specified by `types`. If keyword arguments are provided, they will be added to
the method signature.
"""
get_sig(f::SymOrExpr, arg_names::Vector{Symbol}, types::Vector; kwargs...) =
    add_kwargs!(Expr(:call, f, map((nm, tp)->:($nm::$tp), arg_names, types)...); kwargs...)

"""
    get_body(foo::Symbol, type_tuple::Expr, arg_names::Vector, invoke_type_tuple::Expr; kwargs...)

Get the body of the @generated function which is used to intercept the invocations
specified by type_tuple.
"""
function get_body(
    foo::Symbol,
    type_tuple::Expr,
    arg_names::Vector,
    invoke_type_tuple::Expr;
    kwargs...
)
    quot_arg_names = map(quot, arg_names)
    dots = Symbol("...")

    arg_tuple = any(isa_vararg.(get_types(get_body(type_tuple)))) ?
        Expr(:tuple, arg_names[1:end-1]..., Expr(dots, arg_names[end])) :
        Expr(:tuple, arg_names...)
    sym_arg_tuple = any(isa_vararg.(get_types(get_body(type_tuple)))) ?
        Expr(:tuple, quot_arg_names[1:end-1]..., quot(Expr(dots, arg_names[end]))) :
        Expr(:tuple, quot_arg_names...)

    args_dotted = Expr(dots, Expr(:vect, arg_names...))
    args_dotted_quot = Expr(dots, Expr(:vect, quot_arg_names...))

    branch = :(Nabla.branch_expr($(quot(foo)), is_node, x, x_syms, $(quot(arg_tuple))))
    add_kwargs!(branch; kwargs...)

    invoke = :(Nabla.invoke_expr($(quot(foo)), $(quot(invoke_type_tuple)), x_dots))
    add_kwargs!(invoke; kwargs...)

    return Expr(:block,
        Expr(Symbol("="), :x, Expr(:tuple, args_dotted)),
        Expr(Symbol("="), :x_syms, Expr(:tuple, args_dotted_quot)),
        Expr(Symbol("="), :x_dots, sym_arg_tuple),
        Expr(Symbol("="), :is_node, :([any((<:).(xj, Node)) for xj in x])),
        Expr(:return, Expr(:if, Expr(:call, :any, :is_node), branch, invoke))
    )
end

"""
    branch_expr(foo::Symbol, is_node::Vector{Bool}, x::Tuple, arg_tuple::Expr; kwargs...)

Generate an expression to call Branch.
"""
function branch_expr(
    foo::Symbol,
    is_node::Vector{Bool},
    x::Tuple,
    syms::NTuple{<:Any, Symbol},
    arg_tuple::Expr;
    kwargs...
)
    call = Expr(:call, :Branch, foo, arg_tuple, tape_expr(x, syms, is_node))
    add_kwargs!(call; kwargs...)
    return call
end

invoke_expr(f::Symbol, invoke_tuple::Expr, arg_syms; kwargs...) =
    Expr(:call, :invoke, f, invoke_tuple, arg_syms...; kwargs...)

"""
    tape_expr(x::Tuple, syms::NTuple{N, Symbol} where N, is_node::Vector{Bool})

Get an expression which will obtain the tape from a Node object in `x`.
"""
function tape_expr(x::Tuple, syms::NTuple{N, Symbol} where N, is_node::Vector{Bool})
    idx = findfirst(is_node)
    if idx == length(is_node) && isa(x[end], Tuple)
        node_idx = findfirst([varg <: Node for varg in x[end]])
        return Expr(:call, :getfield, Expr(:ref, syms[end], node_idx), quot(:tape))
    else
        return Expr(:call, :getfield, syms[idx], quot(:tape))
    end
end

"""
    preprocess(f, y, ȳ, xs...) = ()

Default implementation of preprocess returns an empty Tuple. Individual sensitivity
implementations should add methods specific to their use case. The output is passed
in to `∇` as the 3rd or 4th argument in the new-x̄ and update-x̄ cases respectively.

`preprocess` is invoked with `y` and `xs` still boxed.
The default implementation just calls `unbox` on them then calls `preprocess` on the unboxed
values.
If for preprocessing you need the boxed values you should overload
`preprocess(f, y::Node, ȳ, xs...)`.
If you need them unboxed, then overloading `preprocess(f, y, ȳ, xs...)` is fine.
"""
@inline preprocess(f, y, ȳ, xs...) = ()
@inline preprocess(f, y::Node, ȳ, xs...) = preprocess(f, unbox(y), ȳ, map(unbox, xs)...)


text:
@testset "Parallel.Johnson" begin
    g3 = PathGraph(5)
    d = Symmetric([0 1 2 3 4; 1 0 6 7 8; 2 6 0 11 12; 3 7 11 0 16; 4 8 12 16 0])
    for g in testgraphs(g3)
        z = @inferred(Parallel.johnson_shortest_paths(g, d))
        @test z.dists[3, :][:] == [7, 6, 0, 11, 27]
        @test z.parents[3, :][:] == [2, 3, 0, 3, 4]
        @test @inferred(enumerate_paths(z))[2][2] == []
        @test @inferred(enumerate_paths(z))[2][4] == enumerate_paths(z, 2)[4] == enumerate_paths(z, 2, 4) == [2, 3, 4]
    end

    g4 = PathDiGraph(4)
    for g in testdigraphs(g4)
        z = @inferred(Parallel.johnson_shortest_paths(g))
        @test length(enumerate_paths(z, 4, 3)) == 0
        @test length(enumerate_paths(z, 4, 1)) == 0
        @test length(enumerate_paths(z, 2, 3)) == 2
    end

    g5 = DiGraph([1 1 1 0 1; 0 1 0 1 1; 0 1 1 0 0; 1 0 1 1 0; 0 0 0 1 1])
    d = [0 3 8 0 -4; 0 0 0 1 7; 0 4 0 0 0; 2 0 -5 0 0; 0 0 0 6 0]
    for g in testdigraphs(g5)
        z = @inferred(Parallel.johnson_shortest_paths(g, d))
        @test z.dists == [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0]
    end
end


text:
function edit_distance{T}(a, b, submat::AbstractSubstitutionMatrix{T}, insertion_cost::T, deletion_cost::T)
    m = length(a)
    n = length(b)
    trace = Matrix{Trace}(m + 1, n + 1)
    D = Vector{T}(m + 1)
    D[1] = T(0)
    for i in 1:m
        D[i+1] = i * deletion_cost
        trace[i+1,1] = TRACE_INSERT
    end
    for j in 1:n
        d_diag = D[1]
        D[1] = j * insertion_cost
        trace[1,j+1] = TRACE_DELETE
        for i in 1:m
            del = D[i]   + deletion_cost
            ins = D[i+1] + insertion_cost
            mat = d_diag + submat[a[i],b[j]]
            # find the best score and its trace
            best = min(del, ins, mat)
            t = TRACE_NONE
            del == best && (t |= TRACE_DELETE)
            ins == best && (t |= TRACE_INSERT)
            mat == best && (t |= TRACE_MATCH)
            d_diag = D[i+1]
            D[i+1] = best
            trace[i+1,j+1] = t
        end
    end
    return D[m+1], trace
end

function edit_traceback(a, b, trace, endpos)
    # gap/character counts (reversed order)
    counts_a = [0, 0]
    counts_b = [0, 0]
    i, j = endpos
    while i ≥ 1 && j ≥ 2
        t = trace[i+1,j+1]
        if t & TRACE_MATCH > 0
            @match
        elseif t & TRACE_DELETE > 0
            @gap a
        elseif t & TRACE_INSERT > 0
            @gap b
        end
        # do not come here
        @assert false
    end
    while j ≥ 1 @gap a end
    while i ≥ 1 @gap b end
    reverse!(counts_a)
    reverse!(counts_b)
    return GappedSequence(a, counts_a), GappedSequence(b, counts_b)
end


text:
using GLVisualize, GLAbstraction

if !isdefined(:runtests)
	window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end

description = """
This example shows how you can define a function in the OpenGL shader language,
which will get sampled on the GPU on a per pixel basis.
"""

# create a glsl fragment shader
parametric_shader = """
uniform float arg1; // you can add arbitrary uniforms and supply them via the keyword args
float function(float x) {
    return arg1 * sin(1.0/tan(x));
}
"""
# _view the function on a 1700x800 pixel plane
paremetric = visualize(parametric_shader, :shader, arg1 = timesignal, dimensions = (1700, 800))
_view(paremetric, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
module NeuralOperators
    using Flux
    using FFTW
    using Tullio
    using CUDA
    using CUDAKernels
    using KernelAbstractions
    using Zygote
    using ChainRulesCore

    include("fourier.jl")
    include("model.jl")
end


text:
@testset "assembly.jl" begin
    comp1 = Item(Inventory(0, 50), Supplier(100, 10))
    comp2 = Item(Inventory(0, 20), Supplier(100, 10))
    ass = Assembly(100,10, comp1 => 2, comp2 =>1)
    @test state(ass) == []
    @test state_size(ass) == 0
    dest = []
    InventoryModels.pull!(ass, 10, dest)
    comp1(0)
    comp2(0)
    InventoryModels.dispatch!(comp1)
    InventoryModels.dispatch!(comp2)
    InventoryModels.dispatch!(ass)
    InventoryModels.reward!(comp1)
    InventoryModels.reward!(comp2)
    @test InventoryModels.reward!(ass) == -100 - 100
    @test dest == [10, ass.leadtime]
    @test state(ass) == []
    @test state(comp1) == [30]
    @test state(comp2) == [10]
    InventoryModels.pull!(ass, 12, dest)
    comp1(0)
    comp2(0)
    InventoryModels.dispatch!(comp1)
    InventoryModels.dispatch!(comp2)
    InventoryModels.dispatch!(ass)
    InventoryModels.reward!(comp1)
    InventoryModels.reward!(comp2)
    @test InventoryModels.reward!(ass) == -100 - 100
    @test dest == [10, ass.leadtime, 10, ass.leadtime]
    @test state(ass) == []
    @test state(comp1) == [10]
    @test state(comp2) == [0]

    comp1 = Item(Inventory(0, 50), Supplier(100, 10))
    comp2 = Item(Inventory(0, 20), Supplier(100, 10))
    ass = Assembly(100,10, comp1 => 2, comp2 =>1, leadtime = LeadTime(3, 20))
    @test state(ass) == [20,20,20]
    @test state_size(ass) == 3
    dest = []
    InventoryModels.pull!(ass, 10, dest)
    comp1(0)
    comp2(0)
    InventoryModels.dispatch!(comp1)
    InventoryModels.dispatch!(comp2)
    InventoryModels.dispatch!(ass)
    InventoryModels.reward!(comp1)
    InventoryModels.reward!(comp2)
    @test InventoryModels.reward!(ass) == -100 - 100
    @test dest == [20, ass.leadtime]
    @test state(ass) == [20,20,10]
    @test state(comp1) == [30]
    @test state(comp2) == [10]
    InventoryModels.pull!(ass, 12, dest)
    comp1(0)
    comp2(0)
    InventoryModels.dispatch!(comp1)
    InventoryModels.dispatch!(comp2)
    InventoryModels.dispatch!(ass)
    InventoryModels.reward!(comp1)
    InventoryModels.reward!(comp2)
    @test InventoryModels.reward!(ass) == -100 - 100
    @test dest == [20, ass.leadtime, 20, ass.leadtime]
    @test state(ass) == [20,10,10]
    @test state(comp1) == [10]
    @test state(comp2) == [0]
end


text:
#unit-boundstesting.jl
#testing that the endpoint bounds from unums conforms to the mathematical
#expectation
UT = Unum{3,5}

import Unums: inner_exact!, inner_ulp!, outer_exact!, outer_ulp!

left_neginf =   Ubound(neg_inf(UT), UT(-1))
left_negmmr_b = Ubound(neg_mmr(UT), UT(-1))
left_negmmr_u = neg_mmr(UT)
left_exact =    Ubound(UT(-2), UT(-1))
left_ulp =      Ubound(inner_ulp!(UT(-2)), UT(-1))
left_posinf =   inf(UT)

right_neginf =   neg_inf(UT)
right_ulp =      Ubound(UT(1), inner_ulp!(UT(2)))
right_exact =    Ubound(UT(1), UT(2))
right_posmmr_u = mmr(UT)
right_posmmr_b = Ubound(UT(1), mmr(UT))
right_posinf =   Ubound(UT(1), inf(UT))

#testing special ubound addition (NB: p. 113, TEoE).
#TOP TABLE
#top row, left to right.
@test left_neginf + left_neginf   == Ubound(neg_inf(UT), UT(-2))                #[-∞, -1] + [-∞, -1]  == [-∞, -2]
@test left_neginf + left_negmmr_b == Ubound(neg_inf(UT), UT(-2))                #[-∞, -1] + (-∞, -1]  == [-∞, -2]
@test left_neginf + left_negmmr_u == Ubound(neg_inf(UT), neg_mmr(UT))           #[-∞, -1] + (-∞, nmr) == [-∞, nmr)
@test left_neginf + left_exact    == Ubound(neg_inf(UT), UT(-2))                #[-∞, -1] + [-2, -1]  == [-∞, -2]
@test left_neginf + left_ulp      == Ubound(neg_inf(UT), UT(-2))                #[-∞, -1] + (-2, -1]  == [-∞, -2]
@test isnan(left_neginf + left_posinf)                                          #[-∞, -1] + ∞  == NaN
#second row (ubound), left to right
@test left_negmmr_b + left_neginf   == Ubound(neg_inf(UT), UT(-2))              #(-∞, -1] + [-∞, -1]  == [-∞, -2]
@test left_negmmr_b + left_negmmr_b == Ubound(neg_mmr(UT), UT(-2))              #(-∞, -1] + (-∞, -1]  == (-∞, -2]
@test left_negmmr_b + left_negmmr_u == neg_mmr(UT)                              #(-∞, -1] + (-∞, nmr) == (-∞, nmr)
@test left_negmmr_b + left_exact    == Ubound(neg_mmr(UT), UT(-2))              #(-∞, -1] + [-2, -1]  == (-∞, -2]
@test left_negmmr_b + left_ulp      == Ubound(neg_mmr(UT), UT(-2))              #(-∞, -1] + (-2, -1]  == (-∞, -2]
@test left_negmmr_b + left_posinf   == pos_inf(UT)                              #(-∞, -1] + ∞  == ∞
#second row (unum), left to right
@test left_negmmr_u + left_neginf   == Ubound(neg_inf(UT), neg_mmr(UT))         #(-∞, nmr) + [-∞, -1]  == [-∞, nmr)
@test left_negmmr_u + left_negmmr_b == neg_mmr(UT)                              #(-∞, nmr) + (-∞, -1]  == (-∞, nmr)
@test left_negmmr_u + left_negmmr_u == neg_mmr(UT)                              #(-∞, nmr) + (-∞, nmr) == (-∞, nmr)
@test left_negmmr_u + left_exact    == neg_mmr(UT)                              #(-∞, nmr) + [-2, -1]  == (-∞, nmr)
@test left_negmmr_u + left_ulp      == neg_mmr(UT)                              #(-∞, nmr) + (-2, -1]  == (-∞, nmr)
@test left_negmmr_u + left_posinf   == pos_inf(UT)                              #(-∞, nmr) + ∞  == ∞
#third row, left to right
@test left_exact + left_neginf   == Ubound(neg_inf(UT), UT(-2))                 #[-2, -1] + [-∞, -1]  == [-∞, -2]
@test left_exact + left_negmmr_b == Ubound(neg_mmr(UT), UT(-2))                 #[-2, -1] + (-∞, -1]  == (-∞, -2]
@test left_exact + left_negmmr_u == neg_mmr(UT)                                 #[-2, -1] + (-∞, nmr) == (-∞, nmr)
@test left_exact + left_exact    == Ubound(UT(-4), UT(-2))                      #[-2, -1] + [-2, -1]  == [-4, -2]
@test left_exact + left_ulp      == Ubound(inner_ulp!(UT(-4)), UT(-2))          #[-2, -1] + (-2, -1]  == (-4, -2]
@test left_exact + left_posinf   == pos_inf(UT)                                 #[-2, -1] + ∞  == ∞
#fourth row, left to right
@test left_ulp + left_neginf   == Ubound(neg_inf(UT), UT(-2))                   #(-2, -1] + [-∞, -1]  == [-∞, -2]
@test left_ulp + left_negmmr_b == Ubound(neg_mmr(UT), UT(-2))                   #(-2, -1] + (-∞, -1]  == (-∞, -2]
@test left_ulp + left_negmmr_u == neg_mmr(UT)                                   #(-2, -1] + (-∞, nmr) == (-∞, nmr)
@test left_ulp + left_exact    == Ubound(inner_ulp!(UT(-4)), UT(-2))            #(-2, -1] + [-2, -1]  == (-4, -2]
@test left_ulp + left_ulp      == Ubound(inner_ulp!(UT(-4)), UT(-2))            #(-2, -1] + (-2, -1]  == (-4, -2]
@test left_ulp + left_posinf   == pos_inf(UT)                                   #(-2, -1] + ∞  == ∞
#fifth row, left to right
@test isnan(left_posinf + left_neginf)                                          #∞ + [-∞, -1]  == NaN
@test left_posinf + left_negmmr_b == pos_inf(UT)                                #∞ + (-∞, -1]  == ∞
@test left_posinf + left_negmmr_u == pos_inf(UT)                                #∞ + (-∞, nmr) == ∞
@test left_posinf + left_exact    == pos_inf(UT)                                #∞ + [-2, -1]  == ∞
@test left_posinf + left_ulp      == pos_inf(UT)                                #∞ + (-2, -1]  == ∞
@test left_posinf + left_posinf   == pos_inf(UT)                                #∞ + ∞  == ∞

#BOTTOM TABLE
#top row, left to right.
@test right_neginf + right_neginf   == neg_inf(UT)                              #-∞ + -∞      == -∞
@test right_neginf + right_ulp      == neg_inf(UT)                              #-∞ + [1, 2)  == -∞
@test right_neginf + right_exact    == neg_inf(UT)                              #-∞ + [1, 2]  == -∞
@test right_neginf + right_posmmr_u == neg_inf(UT)                              #-∞ + (mr, ∞) == -∞
@test right_neginf + right_posmmr_b == neg_inf(UT)                              #-∞ + [1, ∞]  == -∞
@test isnan(right_neginf + right_posinf)                                        #-∞ + ∞       == NaN
#second row, left to right.
@test right_ulp + right_neginf   == neg_inf(UT)                                 #[1, 2) + -∞      == -∞
@test right_ulp + right_ulp      == Ubound(UT(2), inner_ulp!(UT(4)))            #[1, 2) + [1, 2)  == [2, 4)
@test right_ulp + right_exact    == Ubound(UT(2), inner_ulp!(UT(4)))            #[1, 2) + [1, 2]  == [2, 4)
@test right_ulp + right_posmmr_u == pos_mmr(UT)                                 #[1, 2) + (mr, ∞) == (mr, ∞)
@test right_ulp + right_posmmr_b == Ubound(UT(2), pos_mmr(UT))                  #[1, 2) + [1, ∞)  == [2, ∞)
@test right_ulp + right_posinf   == Ubound(UT(2), pos_inf(UT))                  #[1, 2) + [1 ,∞]  == [2, ∞]
#third row, left to right.
@test right_exact + right_neginf   == neg_inf(UT)                               #[1, 2] + -∞      == -∞
@test right_exact + right_ulp      == Ubound(UT(2), inner_ulp!(UT(4)))          #[1, 2] + [1, 2)  == [2, 4)
@test right_exact + right_exact    == Ubound(UT(2), UT(4))                      #[1, 2] + [1, 2]  == [2, 4]
@test right_exact + right_posmmr_u == pos_mmr(UT)                               #[1, 2] + (mr, ∞) == (mr, ∞)
@test right_exact + right_posmmr_b == Ubound(UT(2), pos_mmr(UT))                #[1, 2] + [1, ∞)  == [2, ∞)
@test right_exact + right_posinf   == Ubound(UT(2), pos_inf(UT))                #[1, 2] + [1 ,∞]  == [2, ∞]
#fourth row (unum), left to right.
@test right_posmmr_u + right_neginf   == neg_inf(UT)                            #(mr, ∞) + -∞      == -∞
@test right_posmmr_u + right_ulp      == pos_mmr(UT)                            #(mr, ∞) + [1, 2)  == (mr, ∞)
@test right_posmmr_u + right_exact    == pos_mmr(UT)                            #(mr, ∞) + [1, 2]  == (mr, ∞)
@test right_posmmr_u + right_posmmr_u == pos_mmr(UT)                            #(mr, ∞) + (mr, ∞) == (mr, ∞)
@test right_posmmr_u + right_posmmr_b == pos_mmr(UT)                            #(mr, ∞) + [1, ∞)  == (mr, ∞)
@test right_posmmr_u + right_posinf   == Ubound(pos_mmr(UT), pos_inf(UT))       #(mr, ∞) + [1 ,∞]  == (mr, ∞]
#fourth row (ubound), left to right.
@test right_posmmr_b + right_neginf   == neg_inf(UT)                            #[1, ∞) + -∞      == -∞
@test right_posmmr_b + right_ulp      == Ubound(UT(2), pos_mmr(UT))             #[1, ∞) + [1, 2)  == [2, ∞)
@test right_posmmr_b + right_exact    == Ubound(UT(2), pos_mmr(UT))             #[1, ∞) + [1, 2]  == [2, ∞)
@test right_posmmr_b + right_posmmr_u == pos_mmr(UT)                            #[1, ∞) + (mr, ∞) == (mr, ∞)
@test right_posmmr_b + right_posmmr_b == Ubound(UT(2), pos_mmr(UT))             #[1, ∞) + [1, ∞)  == [2, ∞)
@test right_posmmr_b + right_posinf   == Ubound(UT(2), pos_inf(UT))             #[1, ∞) + [1 ,∞]  == [2, ∞]
#fifth row, left to right.
@test isnan(right_posinf + right_neginf)                                        #[1, ∞] + -∞      == NaN
@test right_posinf + right_ulp      == Ubound(UT(2), pos_inf(UT))               #[1, ∞] + [1, 2)  == [2, ∞]
@test right_posinf + right_exact    == Ubound(UT(2), pos_inf(UT))               #[1, ∞] + [1, 2]  == [2, ∞]
@test right_posinf + right_posmmr_u == Ubound(pos_mmr(UT), pos_inf(UT))         #[1, ∞] + (mr, ∞) == (mr, ∞]
@test right_posinf + right_posmmr_b == Ubound(UT(2), pos_inf(UT))               #[1, ∞] + [1, ∞)  == [2, ∞]
@test right_posinf + right_posinf   == Ubound(UT(2), pos_inf(UT))               #[1, ∞] + [1, ∞]  == [2, ∞]

#testing special ubound multiplication (NB: p. 130, TEoE)
left_zero_exact = Ubound(UT(0), UT(1))
left_zero_ulp_b = Ubound(sss(UT), UT(1))
left_zero_ulp_u = sss(UT)
left_pos_exact  = Ubound(UT(1), UT(2))
left_pos_ulp    = Ubound(outer_ulp!(UT(1)), UT(2))

right_zero_exact = zero(UT)
#TOP TABLE
#top row, left to right.
@test left_zero_exact * left_zero_exact == Ubound(UT(0), UT(1))                 #[0, 1] * [0, 1]   == [0, 1]
@test left_zero_exact * left_zero_ulp_b == Ubound(UT(0), UT(1))                 #[0, 1] * (0, 1]   == [0, 1]
@test left_zero_exact * left_zero_ulp_u == Ubound(UT(0), sss(UT))               #[0, 1] * (0, ssn) == [0, ssn)
@test left_zero_exact * left_pos_exact  == Ubound(UT(0), UT(2))                 #[0, 1] * [1, 2]   == [0, 2]
@test left_zero_exact * left_pos_ulp    == Ubound(UT(0), UT(2))                 #[0, 1] * (1, 2]   == [0, 2]
@test isnan(left_zero_exact * left_posinf)                                      #[0, 1] * ∞        == NaN
#second row (ubound), left to right.
@test left_zero_ulp_b * left_zero_exact == Ubound(UT(0), UT(1))                 #(0, 1] * [0, 1]   == [0, 1]
@test left_zero_ulp_b * left_zero_ulp_b == Ubound(sss(UT), UT(1))               #(0, 1] * (0, 1]   == (0, 1]
@test left_zero_ulp_b * left_zero_ulp_u == sss(UT)                              #(0, 1] * (0, ssn) == (0, ssn)
@test left_zero_ulp_b * left_pos_exact  == Ubound(sss(UT), UT(2))               #(0, 1] * [1, 2]   == (0, 2]
@test left_zero_ulp_b * left_pos_ulp    == Ubound(sss(UT), UT(2))               #(0, 1] * (1, 2]   == (0, 2]
@test left_zero_ulp_b * left_posinf     == inf(UT)                              #(0, 1] * ∞        == ∞
#second row (unum), left to right.
ssn2 = Unum{3,5}(z64, z64, o16, 0x0007, 0x001e)
@test left_zero_ulp_u * left_zero_exact == Ubound(UT(0), sss(UT))               #(0, ssn) * [0, 1]   == [0, ssn)
@test left_zero_ulp_u * left_zero_ulp_b == sss(UT)                              #(0, ssn) * (0, 1]   == (0, ssn)
@test left_zero_ulp_u * left_zero_ulp_u == sss(UT)                              #(0, ssn) * (0, ssn) == (0, ssn)
@test left_zero_ulp_u * left_pos_exact  == ssn2                                 #(0, ssn) * [1, 2]   == (0, 2*ssn)
@test left_zero_ulp_u * left_pos_ulp    == ssn2                                 #(0, ssn) * (1, 2]   == (0, 2*ssn)
@test left_zero_ulp_u * left_posinf     == inf(UT)                              #(0, ssn) * ∞        == ∞
#third row, left to right.
@test left_pos_exact * left_zero_exact == Ubound(UT(0), UT(2))                  #[1, 2] * [0, 1]   == [0, 2]
@test left_pos_exact * left_zero_ulp_b == Ubound(sss(UT), UT(2))                #[1, 2] * (0, 1]   == (0, 2]
@test left_pos_exact * left_zero_ulp_u == ssn2                                  #[1, 2] * (0, ssn) == (0, 2 * ssn)
@test left_pos_exact * left_pos_exact  == Ubound(UT(1), UT(4))                  #[1, 2] * [1, 2]   == [1, 4]
@test left_pos_exact * left_pos_ulp    == Ubound(outer_ulp!(UT(1)), UT(4))      #[1, 2] * (1, 2]   == (1, 4]
@test left_pos_exact * left_posinf     == inf(UT)                               #[1, 2] * ∞        == ∞
#fourth row, left to right.
@test left_pos_ulp * left_zero_exact == Ubound(UT(0), UT(2))                    #(1, 2] * [0, 1]   == [0, 2]
@test left_pos_ulp * left_zero_ulp_b == Ubound(sss(UT), UT(2))                  #(1, 2] * (0, 1]   == (0, 2]
@test left_pos_ulp * left_zero_ulp_u == ssn2                                    #(1, 2] * (0, ssn) == (0, 2* ssn)
@test left_pos_ulp * left_pos_exact  == Ubound(outer_ulp!(UT(1)), UT(4))        #(1, 2] * [1, 2]   == (1, 4]
@test left_pos_ulp * left_pos_ulp    == Ubound(outer_ulp!(UT(1)), UT(4))        #(1, 2] * (1, 2]   == (1, 4]
@test left_pos_ulp * left_posinf     == inf(UT)                                 #(1, 2] * ∞        == ∞
#fifth row, left to right.
@test isequal(left_posinf * left_zero_exact, UT(NaN))                           # ∞ * [0, 1]   == NaN
@test left_posinf * left_zero_ulp_b == inf(UT)                                  # ∞ * (0, 1]   == ∞
@test left_posinf * left_zero_ulp_u == inf(UT)                                  # ∞ * (0, ssn) == ∞
@test left_posinf * left_pos_exact  == inf(UT)                                  # ∞ * [1, 2]   == ∞
@test left_posinf * left_pos_ulp    == inf(UT)                                  # ∞ * (1, 2]   == ∞
@test left_posinf * left_posinf     == inf(UT)                                  # ∞ * ∞        == ∞

#BOTTOM TABLE
#top row, left to right.
@test right_zero_exact * right_zero_exact == zero(UT)                           # 0 * 0        == 0
@test right_zero_exact * right_ulp        == zero(UT)                           # 0 * [1, 2)   == 0
@test right_zero_exact * right_exact      == zero(UT)                           # 0 * [1, 2]   == 0
@test right_zero_exact * right_posmmr_u   == zero(UT)                           # 0 * (mr, ∞)  == 0
@test right_zero_exact * right_posmmr_b   == zero(UT)                           # 0 * [1, ∞)   == 0
@test isequal(right_zero_exact * right_posinf, UT(NaN))                         # 0 * [1, ∞]   == NaN
#second row, left to right.
@test right_ulp * right_zero_exact == zero(UT)                                  #[1, 2) * 0        == 0
@test right_ulp * right_ulp        == Ubound(UT(1), inner_ulp!(UT(4)))          #[1, 2) * [1, 2)   == [1, 4)
@test right_ulp * right_exact      == Ubound(UT(1), inner_ulp!(UT(4)))          #[1, 2) * [1, 2]   == [1, 4)
@test right_ulp * right_posmmr_u   == mmr(UT)                                   #[1, 2) * (mr, ∞)  == (mr, ∞)
@test right_ulp * right_posmmr_b   == right_posmmr_b                            #[1, 2) * [1, ∞)   == [1, ∞)
@test right_ulp * right_posinf     == right_posinf                              #[1, 2) * [1, ∞]   == [1, ∞]
#third row, left to right.
@test right_exact * right_zero_exact == zero(UT)                                #[1, 2] * 0        == 0
@test right_exact * right_ulp        == Ubound(UT(1), inner_ulp!(UT(4)))        #[1, 2] * [1, 2)   == [1, 4)
@test right_exact * right_exact      == Ubound(UT(1), UT(4))                    #[1, 2] * [1, 2]   == [1, 4]
@test right_exact * right_posmmr_u   == mmr(UT)                                 #[1, 2] * (mr, ∞)  == (mr, ∞)
@test right_exact * right_posmmr_b   == right_posmmr_b                          #[1, 2] * [1, ∞)   == [1, ∞)
@test right_exact * right_posinf     == right_posinf                            #[1, 2] * [1, ∞]   == [1, ∞]
#fourth row (unum), left to right.
@test right_posmmr_u * right_zero_exact == zero(UT)                             #(mr, ∞) * 0        == 0
@test right_posmmr_u * right_ulp        == mmr(UT)                              #(mr, ∞) * [1, 2)   == (mr, ∞)
@test right_posmmr_u * right_exact      == mmr(UT)                              #(mr, ∞) * [1, 2]   == (mr, ∞)
@test right_posmmr_u * right_posmmr_u   == mmr(UT)                              #(mr, ∞) * (mr, ∞)  == (mr, ∞)
@test right_posmmr_u * right_posmmr_b   == mmr(UT)                              #(mr, ∞) * [1, ∞)   == (mr, ∞)
@test right_posmmr_u * right_posinf     == Ubound(mmr(UT), inf(UT))             #(mr, ∞) * [1, ∞]   == (mr, ∞]
#fourth row (ubound), left to right.
@test right_posmmr_b * right_zero_exact == zero(UT)                             #[1, ∞) * 0        == 0
@test right_posmmr_b * right_ulp        == right_posmmr_b                       #[1, ∞) * [1, 2)   == [1, ∞)
@test right_posmmr_b * right_exact      == right_posmmr_b                       #[1, ∞) * [1, 2]   == [1, ∞)
@test right_posmmr_b * right_posmmr_u   == mmr(UT)                              #[1, ∞) * (mr, ∞)  == (mr, ∞)
@test right_posmmr_b * right_posmmr_b   == right_posmmr_b                       #[1, ∞) * [1, ∞)   == [1, ∞)
@test right_posmmr_b * right_posinf     == Ubound(UT(1), inf(UT))               #[1, ∞) * [1, ∞]   == [1, ∞]
#fifth row, left to right.
@test isequal(right_posinf * right_zero_exact, UT(NaN))                         #[1, ∞] * 0        == NaN
@test right_posinf * right_ulp        == right_posinf                           #[1, ∞] * [1, 2)   == [1, ∞]
@test right_posinf * right_exact      == right_posinf                           #[1, ∞] * [1, 2]   == [1, ∞]
@test right_posinf * right_posmmr_u   == Ubound(mmr(UT), inf(UT))               #[1, ∞] * (mr, ∞)  == (mr, ∞]
@test right_posinf * right_posmmr_b   == right_posinf                           #[1, ∞] * [1, ∞)   == [1, ∞]
@test right_posinf * right_posinf     == right_posinf                           #[1, ∞] * [1, ∞]   == [1, ∞]

#testing special ubound division (NB: p. 138, TEoE)
#NB:  Why is dividing by zero always NaN except in the closed-bounds case?
#because without adjacent values, the signedness of the resulting infinity is not
#clear.

#special case:
@test isequal(left_zero_exact / left_zero_exact, UT(NaN))                       #[0, 1] / [0, 1]  == NaN
#if you want ops to throw infinities, do it by intercepting your function (eg pow(x, [even int]))
#CHART ONE
#top row.
@test isequal(left_zero_exact / right_zero_exact, UT(NaN))                      #[0, 1] / 0       == NaN
@test left_zero_exact / right_ulp            == Ubound(UT(0), UT(1))            #[0, 1] / [1, 2)  == [0, 1]
@test left_zero_exact / right_exact          == Ubound(UT(0), UT(1))            #[0, 1] / [1, 2]  == [0, 1]
@test left_zero_exact / right_posmmr_b       == Ubound(UT(0), UT(1))            #[0, 1] / [1, ∞)  == [0, 1]
@test left_zero_exact / right_posinf         == Ubound(UT(0), UT(1))            #[0, 1] / [1, ∞]  == [0, 1]
################################################################################
#second row
@test isequal(left_zero_ulp_b / right_zero_exact, UT(NaN))                      #(0, 1] / 0       == NaN
@test left_zero_ulp_b / right_ulp            == Ubound(sss(UT), UT(1))          #(0, 1] / [1, 2)  == (0, 1]
@test left_zero_ulp_b / right_exact          == Ubound(sss(UT), UT(1))          #(0, 1] / [1, 2]  == (0, 1]
@test left_zero_ulp_b / right_posmmr_b       == Ubound(sss(UT), UT(1))          #(0, 1] / [1, ∞)  == (0, 1]
@test left_zero_ulp_b / right_posinf         == Ubound(UT(0), UT(1))            #(0, 1] / [1, ∞]  == [0, 1]
################################################################################
#third row
@test isequal(left_pos_exact / right_zero_exact, UT(NaN))                       #[1, 2] / 0       == NaN
@test left_pos_exact / right_ulp         == Ubound(outer_ulp!(UT(0.5)), UT(2))  #[1, 2] / [1, 2)  == (0.5, 2]
@test left_pos_exact / right_exact       == Ubound(UT(0.5), UT(2))              #[1, 2] / [1, 2]  == [0.5, 2]
@test left_pos_exact / right_posmmr_b    == Ubound(sss(UT), UT(2))              #[1, 2] / [1, ∞)  == (0, 2]
@test left_pos_exact / right_posinf      == Ubound(UT(0), UT(2))                #[1, 2] / [1, ∞]  == [0, 2]
################################################################################
#fourth row
@test isequal(left_pos_ulp / right_zero_exact, UT(NaN))                         #(1, 2] / 0       == NaN
@test left_pos_ulp / right_ulp            == Ubound(outer_ulp!(UT(0.5)), UT(2)) #(1, 2] / [1, 2)  == (0.5, 2]
@test left_pos_ulp / right_exact          == Ubound(outer_ulp!(UT(0.5)), UT(2)) #(1, 2] / [1, 2]  == (0.5, 2]
@test left_pos_ulp / right_posmmr_b       == Ubound(sss(UT), UT(2))             #(1, 2] / [1, ∞)  == (0, 2]
@test left_pos_ulp / right_posinf         == Ubound(UT(0), UT(2))               #(1, 2] / [1, ∞]  == [0, 2]
################################################################################
#fifith row
@test isequal(left_posinf / right_zero_exact, UT(NaN))                          # ∞ / 0       == NaN
@test left_posinf / right_ulp            == inf(UT)                             # ∞ / [1, 2)  == ∞
@test left_posinf / right_exact          == inf(UT)                             # ∞ / [1, 2]  == ∞
@test left_posinf / right_posmmr_b       == inf(UT)                             # ∞ / [1, ∞)  == ∞
@test isequal(left_posinf / right_posinf, nan(UT))                              # ∞ / [1, ∞]  == NaN
#CHART TWO is covered by chart one.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
mutable struct ConstantMean{T<:Real,O} <: PriorMean{T}
    C::Vector{T}
    opt::O
end

"""
    ConstantMean(c::Real = 1.0; opt=ADAM(0.01))

## Arguments
- `c::Real` : Constant value

Construct a prior mean with constant `c`
Optionally set an optimiser `opt` (`ADAM(0.01)` by default)
"""
function ConstantMean(c::T=1.0; opt=ADAM(0.01)) where {T<:Real}
    return ConstantMean{T,typeof(opt)}([c], opt)
end

function Base.show(io::IO, ::MIME"text/plain", μ₀::ConstantMean)
    return print(io, "Constant Mean Prior (c = ", first(μ₀.C), ")")
end

(μ::ConstantMean{T})(::Real) where {T<:Real} = first(μ.C)
(μ::ConstantMean{T})(x::AbstractMatrix) where {T<:Real} = fill(first(μ.C), size(x, 1))

function update!(μ::ConstantMean{T}, grad::AbstractVector) where {T<:Real}
    return μ.C .+= Optimise.apply!(μ.opt, μ.C, grad)
end


text:
using PyCall, Interp
@pyimport scipy.interpolate as scii

##  pixel grid and Fourier conjugate grid in radians
immutable Grid_xandk
    pixel_size_arcmin::Float64
    n::Float64
    deltx::Float64
    deltk::Float64
    period::Float64
    nyq::Float64
    kside::Array{Float64,1}
    x::Array{Float64,2}
    y::Array{Float64,2}
    k1::Array{Float64,2}
    k2::Array{Float64,2}
    r::Array{Float64,2}
end
function Grid_xandk(pixel_size_arcmin, n)
    deltx = pixel_size_arcmin * pi / (180 * 60) #this is in radians
    x,y = meshgrid(1:n, 1:n)
    x .*= deltx
    x .-= minimum(x)
    y .*= deltx
    y .-= minimum(y)
    kside  = (2 * pi / deltx) .* (mod(1/2 .+ [0:(n-1)]./n, 1) .- 1/2)
    k1,k2  = meshgrid(kside, kside)
    period = deltx * n
    deltk  = 2 * pi / period
    nyq    = 2 * pi / (2 * deltx)
    r      = sqrt(k1.^2 .+ k2.^2)
    return Grid_xandk(pixel_size_arcmin, n, deltx, deltk, period, nyq, kside, x, y, k1, k2, r)
end

immutable SpectrumGrids # this contains all the signal + noise spectral information
    grd::Grid_xandk # everything is defined in reference to a grid
    ell::Array{Float64,1}
    ellLen::Array{Float64,1}
    #------------------------
    cMaskBool::BitArray{2} # this gives l_max for cmb
    pMaskBool::BitArray{2} # this gives l_max for phi
    cMaskInf::Array{Float64,2}
    #---------- noise
    nugget_at_each_pixel::Float64
    beamFWHM::Float64
    beam::Array{Float64,2}
    # radial profile
    CTell2d::Array{Float64,1}
    CTell2dLen::Array{Float64,1}
    CPell2d::Array{Float64,1}
    # matrices
    cTT::Array{Float64,2}
    cTTLen::Array{Float64,2}
    cPP::Array{Float64,2}
    cNT::Array{Float64,2}
end
function SpectrumGrids(
    ell,
    ellLen,
    CPell2d,
    CTell2d,
    CTell2dLen,
    pixel_size_arcmin,
    n,
    beamFWHM,
    nugget_at_each_pixel,
    maskupC,
    maskupP
)
    d = 2.0
    grd = Grid_xandk(pixel_size_arcmin, n)
    sig=(beamFWHM * (pi / (180 * 60))) / ( 2 * sqrt(2 * log(2)))
    beam = exp(-(sig^2) .* (grd.k1.^2 .+ grd.k2.^2) ./ 2)
    beamSQ =abs2(beam)
    cNT    = nugget_at_each_pixel .* (grd.deltx^d) ./ beamSQ

    cMaskBool = !((grd.r.<= maskupC) & (grd.r.>= 1.0)) # this is true when the masking is active Mat[cMaskBool] = 0 or Inf
    pMaskBool = !((grd.r.<= maskupP) & (grd.r.>= 1.0)) # this is true when the masking is active Mat[cMaskBool] = 0 or Inf
    cMaskInf = 1.0./((grd.r.<=maskupC) & (grd.r.>= 1.0))  # this equals Inf when the masking is active and 1.0 otherwise...

    index=ceil(grd.r)
    index[find(index.==0)]=1

    logCPP = linear_interp1(ell,log(CPell2d), index)
    logCPP[find(logCPP .== 0)] =-Inf
    logCPP[find(isnan(logCPP))]=-Inf
    cPP = exp(logCPP)
    logCTT = linear_interp1(ell,log(CTell2d), index)
    logCTT[find(logCTT .== 0)] =-Inf
    logCTT[find(isnan(logCTT))]=-Inf
    cTT = exp(logCTT)
    logCTTLen = linear_interp1(ellLen,log(CTell2dLen), index)
    logCTTLen[find(logCTTLen .== 0)] =-Inf
    logCTTLen[find(isnan(logCTTLen))]=-Inf
    logCTTLen[grd.r .> 8000]=-Inf
    cTTLen = exp(logCTTLen)

    inputs = {grd,
        ell,
        ellLen,
        cMaskBool,
        pMaskBool,
        cMaskInf,
        nugget_at_each_pixel,
        beamFWHM,
        beam,
        CTell2d,
        CTell2dLen,
        CPell2d,
        cTT,
        cTTLen,
        cPP,
        cNT }
    SpectrumGrids(inputs...)
end

# get spectra and an instance of SpectrumGrids
function setpar(pixel_size_arcmin, n, beamFWHM, nugget_at_each_pixel, maskupC, maskupP)
    Sfile = open("src/camb/test_scalCls.dat")
    lines = readlines(Sfile)
    close(Sfile)
    CLs = Array(Any,(length(lines),6))
    for i in 1:length(lines)
        CLs[i,:] = lines[i] |> strip |> chomp |> split |> transpose
    end
    CLs = map(parsefloat, CLs)
    ell=[1;CLs[:,1]]
    CTell2d=[0;CLs[:,2].*(2*pi)./(CLs[:,1].*(CLs[:,1] .+ 1))]
    CPell2d=[0;CLs[:,5]./(7.4311e12*CLs[:,1].^4)]

    LenSfile = open("src/camb/test_lensedCls.dat")
    Lenlines = readlines(LenSfile)
    close(LenSfile)
    LenCLs = Array(Any,(length(Lenlines),5))
    for i in 1:length(Lenlines)
        LenCLs[i,:] = Lenlines[i] |> strip |> chomp |> split |> transpose
    end
    LenCLs = map(parsefloat, LenCLs)  #[ell, CTT, CEE, CBB, CTE]
    ellLen = [1;LenCLs[:,1]]
    CTell2dLen =[0;LenCLs[:,2].*(2*pi)./(LenCLs[:,1].*(LenCLs[:,1].+1))]

   inputs = {ell,
        ellLen,
        CPell2d,
        CTell2d,
        CTell2dLen,
        pixel_size_arcmin,
        n,
        beamFWHM,
        nugget_at_each_pixel,
        maskupC,
        maskupP}
    SpectrumGrids(inputs...)
end

#  if you need to specify where src is...used in graphics
function setpar(pixel_size_arcmin, n, beamFWHM, nugget_at_each_pixel, maskupC, maskupP, pathtosrc::String)
    Sfile = open(joinpath(pathtosrc,"camb/test_scalCls.dat"))
    lines = readlines(Sfile)
    close(Sfile)
    CLs = Array(Any,(length(lines),6))
    for i in 1:length(lines)
        CLs[i,:] = lines[i] |> strip |> chomp |> split |> transpose
    end
    CLs = map(parsefloat, CLs)
    ell=[1;CLs[:,1]]
    CTell2d=[0;CLs[:,2].*(2*pi)./(CLs[:,1].*(CLs[:,1] .+ 1))]
    CPell2d=[0;CLs[:,5]./(7.4311e12*CLs[:,1].^4)]

    LenSfile = open(joinpath(pathtosrc,"camb/test_lensedCls.dat"))
    Lenlines = readlines(LenSfile)
    close(LenSfile)
    LenCLs = Array(Any,(length(Lenlines),5))
    for i in 1:length(Lenlines)
        LenCLs[i,:] = Lenlines[i] |> strip |> chomp |> split |> transpose
    end
    LenCLs = map(parsefloat, LenCLs)  #[ell, CTT, CEE, CBB, CTE]
    ellLen = [1;LenCLs[:,1]]
    CTell2dLen =[0;LenCLs[:,2].*(2*pi)./(LenCLs[:,1].*(LenCLs[:,1].+1))]

   inputs = {ell,
        ellLen,
        CPell2d,
        CTell2d,
        CTell2dLen,
        pixel_size_arcmin,
        n,
        beamFWHM,
        nugget_at_each_pixel,
        maskupC,
        maskupP}
    SpectrumGrids(inputs...)
end

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
# include the package
using LatticePhysics

# obtain the unitcell of the square octagon lattice (kitaev bonds)
unitcell = getUnitcellSquareOctagon(5)

# build a lattice for visual confirmation that this is the correct unitcell
lattice = getLatticeByBondDistance(unitcell, 7)
# plot the lattice for visual confirmation (including bond coloring)
plotLattice2D(lattice, openfile=true, colorcode_bonds_automation="COLOR")


text:
"""
    function ll_twogroup_newvariant_infboost(θ,model::KenyaCoVSD.CoVAreaModel,seroreversionrate;σ = 0.16,ω = 1/180,ct_min2 = 0.445)

Log-likelihood of the case and serology data in `model` assuming that new variant is more transmissible rather than immune evading.
"""
function ll_twogroup_newvariant_infboost(θ,model::KenyaCoVSD.CoVAreaModel,seroreversionrate;σ = 0.16,ω = 1/180,ct_min2 = 0.445)
    @unpack PCR_cases,sero_cases,baseline_sero_array,PCR_array,sero_sensitivity,sero_specificity,N,M_BB,prob,α,γ,relative_testing_rate = model
    @unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    #Set variance scalers
    clustering_factor_PCR = 0.5
    M_PCR = 30.
    T = eltype(R₀)

    #Set transmission parameters
	N₁ = P_eff*N
	N₂ = N - N₁
    p = convert.(T,[R₀,α,γ,ϵ,σ,N₁,N₂,ω,schooleffect,ct_min1,ct_min2])
    u0 = convert.(T,[N₁,E₀,0.,0.,0.,0.,
                     N₂,E₀,0.,0.,0.,0.,
                     0.0,0.0])

    #Create new variant introduction
    function new_variant_effect!(integrator)
        integrator.p[1] *= extra_transmissibility
        integrator.u[2] += influx_exposed_new_variant #Lower SES
        integrator.u[8] += influx_exposed_new_variant #Higher SES
    end

    janendpoint = (Date(2021,1,30) - Date(2020,2,20)).value
	aprilendpoint = (Date(2021,4,30) - Date(2020,2,20)).value
    variant_cb = PresetTimeCallback([janendpoint],new_variant_effect!)

    #Sero-waning
    sero_array = vcat(baseline_sero_array[1:30],[(1-seroreversionrate)^k for k in 1:length(baseline_sero_array[31:end])])

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    LL = T(0.)
    try
        sol = solve(prob, BS3();tspan = (0,aprilendpoint),callback = variant_cb,u0=u0, p=p, saveat = 1)

        ι₁ = diff(sol[:C₁])
        ι₂ = diff(sol[:C₂])
        ι_sero₁ = diff(sol[:C_sero₁])
        ι_sero₂ = diff(sol[:C_sero₂])
        PCR₁ = KenyaCoVSD.simple_conv(ι₁,PCR_array)
        PCR₂ = KenyaCoVSD.simple_conv(ι₂,PCR_array)
    	sero₁ = sero_sensitivity.*KenyaCoVSD.simple_conv(ι_sero₁,sero_array)
	    sero₂ = sero_sensitivity.*KenyaCoVSD.simple_conv(ι_sero₂,sero_array)

        #Calculate log-likelihood for PCR testing
        for t in 55:janendpoint
            #Convert from μ,α parameterisation to p,r parameterisation for using Neg. binomial model
            μ = relative_testing_rate[t]*(p_test₁*PCR₁[t] + p_test₂*PCR₂[t]) + 0.001
            σ² = μ + clustering_factor_PCR*μ^2
            p_negbin = 1 - (clustering_factor_PCR*μ^2/σ²)
            r_negbin  = 1/clustering_factor_PCR
            LL += logpdf(NegativeBinomial(r_negbin,p_negbin),PCR_cases[t,1])#likelihood contribution from PCR testing --- positive case detection

            if PCR_cases[t,2] >= 0 #Negative tests available
                p_PCR_pred = p_choose1*(χ₁*PCR₁[t]/((χ₁-1)*PCR₁[t] + sum(N₁))) + (1-p_choose1)*(χ₂*PCR₂[t]/((χ₂-1)*PCR₂[t] + sum(N₂))) + 0.001
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                LL += logpdf(BetaBinomial(PCR_cases[t,2],p_PCR_pred*M_PCR,(1-p_PCR_pred)*M_PCR),PCR_cases[t,1])#likelihood contribution from PCR testing --- proportion postive
            end
        end

        #Calculate log-likelihood for PCR testing after new variant introduction
        for t in (janendpoint+1):min(size(PCR_cases,1),size(PCR₁,1))
            #Convert from μ,α parameterisation to p,r parameterisation for using Neg. binomial model
            μ = p_test_boost*relative_testing_rate[t]*(p_test₁*PCR₁[t] + p_test₂*PCR₂[t]) + 0.001
            σ² = μ + clustering_factor_PCR*μ^2
            p_negbin = 1 - (clustering_factor_PCR*μ^2/σ²)
            r_negbin  = 1/clustering_factor_PCR
            LL += logpdf(NegativeBinomial(r_negbin,p_negbin),PCR_cases[t,1])#likelihood contribution from PCR testing --- positive case detection

            if PCR_cases[t,2] >= 0 #Negative tests available
                p_PCR_pred = p_choose1*(χ_boost*χ₁*PCR₁[t]/((χ_boost*χ₁-1)*PCR₁[t] + sum(N₁))) + (1-p_choose1)*(χ_boost*χ₂*PCR₂[t]/((χ_boost*χ₂-1)*PCR₂[t] + sum(N₂))) + 0.001
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                LL += logpdf(BetaBinomial(PCR_cases[t,2],p_PCR_pred*M_PCR,(1-p_PCR_pred)*M_PCR),PCR_cases[t,1])#likelihood contribution from PCR testing that day in that age group
            end
        end
        #Calculate log-likelihood contribution from serological testing
        for t in 1:min(size(sero_cases,1),size(sero₁,1))
            #Convert from the p_hat,M_BB parameterisation of the Betabinomial distribution to the α_BB,β_BB parameterisation

            #Overall serological testing
            if sero_cases[t,1,1] + sero_cases[t,2,1] > 0
                p_sero_pred  = (sero₁[t] + sero₂[t])./(N₁ + N₂)
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,1] + sero_cases[t,2,1],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,1])#Likelihood contribution from sero testing
            end
            #Lower SES serological testing
            if sero_cases[t,1,2] + sero_cases[t,2,2] > 0
                p_sero_pred  = sero₁[t]./N₁
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,2] + sero_cases[t,2,2],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,2])#Likelihood contribution from sero testing
            end
            #Higher SES serological testing
            if sero_cases[t,1,3] + sero_cases[t,2,3] > 0
                p_sero_pred  = sero₂[t]./N₂
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,3] + sero_cases[t,2,3],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,3])#Likelihood contribution from sero testing
            end

        end
    catch
        LL = T(-Inf)
    end
    return LL::T
end


text:
module Generators

using TOML

using ..Clang
using ..Clang.LibClang
using ..Clang:
    isBitField,
    isCursorDefinition,
    isMacroBuiltin,
    isMacroFunctionLike,
    isVariadic,
    getAlignOf,
    getArgType,
    getArgument,
    getCursorType,
    getCanonicalType,
    getCursorResultType,
    getElementType,
    getEnumDeclIntegerType,
    getFieldDeclBitWidth,
    getIncludedFile,
    getNamedType,
    getNumArguments,
    getNumElements,
    getOffsetOf,
    getOffsetOfField,
    getPointeeType,
    getSizeOf,
    getTranslationUnit,
    getTranslationUnitCursor,
    getTypedefDeclUnderlyingType,
    getTypeDeclaration,
    hasAttrs

using ..JLLEnvs
using ..JLLEnvs: get_system_dirs, triple2target

include("utils.jl")

include("jltypes.jl")
export AbstractJuliaType, AbstractJuliaSIT, AbstractJuliaSDT
export tojulia

include("definitions.jl")
export get_definition
export add_definition
export @add_def

include("types.jl")
export AbstractExprNodeType
export AbstractFunctionNodeType, AbstractTypedefNodeType, AbstractMacroNodeType
export AbstractStructNodeType, AbstractUnionNodeType, AbstractEnumNodeType
export ExprNode, ExprDAG
export get_nodes, get_exprs

include("macro.jl")
include("top_level.jl")
include("system_deps.jl")
include("nested.jl")
include("resolve_deps.jl")
include("preprocessing.jl")
include("documentation.jl")

include("translate.jl")
export translate

include("codegen.jl")
export emit!

include("print.jl")
export pretty_print

include("audit.jl")
export report_default_tag_types

include("mutability.jl")

include("passes.jl")
export AbstractPass
export Audit
export CatchDuplicatedAnonymousTags
export Codegen
export CodegenMacro
export CodegenPostprocessing
export CodegenPreprocessing
export CollectDependentSystemNode
export CollectNestedRecord
export CollectTopLevelNode
export CommonPrinter
export DeAnonymize
export EpiloguePrinter
export FindOpaques
export FunctionPrinter
export GeneralPrinter
export IndexDefinition
export LinkTypedefToAnonymousTagType
export ProloguePrinter
export RemoveCircularReference
export ResolveDependency
export StdPrinter
export TweakMutability
export TopologicalSort

include("context.jl")
export AbstractContext, Context
export parse_header!, parse_headers!
export create_context
export build!, BUILDSTAGE_ALL, BUILDSTAGE_NO_PRINTING, BUILDSTAGE_PRINTING_ONLY
export get_triple, get_default_args, detect_headers, find_dependent_headers
export get_identifier_node, get_tagtype_node

include("option.jl")
export load_options

function __init__()
   reset_definition()
end

end # module


text:
"Dual Soft Booleans (supports negation)"
struct DualSoftBool{S} <: AbstractSoftBool
  b0::S
  b1::S
end

# using ZenUtils
function DualSoftBool(a::A, b::B) where {A, B}
  a_, b_ = promote(a.logerr, b.logerr)
  DualSoftBool(SoftBool(a_), SoftBool(b_))
end

Base.:!(x::DualSoftBool{S}) where S = DualSoftBool{S}(x.b1, x.b0)
Bool(x::DualSoftBool) = Bool(x.b1)
err(x::DualSoftBool) = err(x.b1)
logerr(x::DualSoftBool) = logerr(x.b1)

dsofteq(x, y, k = globalkernel()) = DualSoftBool(SoftBool(x == y ? k(1.0) : 0.0), ssofteq(x, y, k)) # FIXME
dsoftgt(x, y, k = globalkernel()) = DualSoftBool(ssoftlt(x, y, k), ssoftgt(x, y, k))
dsoftlt(x, y, k = globalkernel()) = DualSoftBool(ssoftgt(x, y, k), ssoftlt(x, y, k))

Base.:&(x::DualSoftBool, y::DualSoftBool) = DualSoftBool(x.b0 | y.b0, x.b1 & y.b1)
Base.:|(x::DualSoftBool, y::DualSoftBool) = DualSoftBool(x.b0 & y.b0, x.b1 | y.b1)

dsofttrue(::Type{T} = Float64) where T = DualSoftBool(SoftBool(-inf(T)), SoftBool(zero(T)))
dsoftfalse(::Type{T} = Float64) where T = DualSoftBool(SoftBool(zero(T)), SoftBool(-inf(T)))

# Arithmetic #
Base.:*(x::DualSoftBool{SoftBool{T}}, y::T) where T <: Real = DualSoftBool(x.b0 * y, x.b1 * y)
Base.:*(x::T, y::DualSoftBool{SoftBool{T}}) where T <: Real = DualSoftBool(x * y.b0, x * y.b1)


text:
function empiricalpotential(x, q, t, a, b, c, pbc, rcutmax, eta, kappa, potential)

# compute full neighbor list
y, alist, neighlist, neighnum = fullneighborlist(x, a, b, c, pbc, rcutmax);
if !isempty(q)
    q = q[:,alist];
end

dim, N = size(x);
e = zeros(N);
f = zeros(dim,N);
v = zeros(6,N);
for i = 1:length(potential)

    style, bondtype, bondedatomtypes, potentialfunc, mu, rcut = getpotential(potential[i]);
    potentialfunc = getfield(Main, Symbol(potentialfunc));

    # single potentials
    if (style == "single") & (bondtype == "nonbonded")
        ilist = Array(1:N);
     elseif (style == "single") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        ilist = findatomtype(Array(1:N), t, typei);
    end
    if (style == "single")
        xi, qi, ai, ti = neighsingles(x, q, t, ilist);
        ei, fi = potentialfunc(xi, qi, ti, mu, eta, kappa);
        ea, fa = tallysingle(ei, fi, ai);
        va = tallysinglevirial(fi, xi, ai, -1.0)
    end

    # pair potentials
    if (style == "pair") & (bondtype == "nonbonded")
        ilist = Int64.(Array(1:N));
        pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
     elseif (style == "pair") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        typej = bondedatomtypes[2];

        ilist = findatomtype(Array(1:N), t, typei);
        pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, typej, neighlist, neighnum, rcut*rcut);
    end
    if (style == "pair")
        xij, ai, aj, ti, tj, qi, qj = neighpairs(y, q, pairlist, pairnum, t, ilist, alist);
        eij, fij = potentialfunc(xij, qi, qj, ti, tj, mu, eta, kappa);
        ea, fa = tallypair(eij, fij, ai, aj);
        va = tallypairvirial(fij, xij, ai, aj, -1.0/2)
    end

    # triplet potentials
    if (style == "triplet") & (bondtype == "nonbonded")
        ilist = Array(1:N);
        pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
     elseif (style == "triplet") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        typej = bondedatomtypes[2];
        typek = bondedatomtypes[3];

        ilist = findatomtype(Array(1:N), t, typei);
        pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, [typej; typek], neighlist, neighnum, rcut*rcut);
    end
    if (style == "triplet")
        tripletlist, tripletnum = neightripletlist(pairlist, pairnum, ilist);
        xij, xik, ai, aj, ak, ti, tj, tk, qi, qj, qk = neightriplets(y, q, tripletlist, tripletnum, atomtype, ilist, alist);
        eij, fij, fik = potentialfunc(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
        ea, fa = tallytriplet(eij, fij, fik, ai, aj, ak);
        va = tallytripletvirial(fij, fik, xij, xik, ai, aj, ak, -1.0/3)
    end

    # quadlet potentials
    if (style == "quadlet") & (bondtype == "nonbonded")
        ilist = Array(1:N);
        pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
     elseif (style == "quadlet") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        typej = bondedatomtypes[2];
        typek = bondedatomtypes[3];
        typel = bondedatomtypes[4];

        ilist = findatomtype(Array(1:N), t, typei);
        pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, [typej typek typel], neighlist, neighnum, rcut*rcut);
    end
    if (style == "quadlet")
        quadletlist, quadletnum = neighquadletlist(pairlist, pairnumsum, ilist);
        xij, xik, xil, ai, aj, ak, al, ti, tj, tk, tl, qi, qj, qk, ql =
                    neighquadlets(y, q, quadletlist, quadletnum, atomtype, ilist, alist);
        eij, fij, fik, fil = potentialfunc(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu, eta, kappa);
        ea, fa = tallyquadlet(eij, fij, fik, fil, ai, aj, ak, al);
        va = tallyquadletvirial(fij, fik, fil, xij, xik, xil, ai, aj, ak, al, -1.0/4)
    end

    e = e .+ ea[:];
    f = f .+ fa;
    v = v .+ va;
end

etot = sum(e);
vtot = sum(v,dims=2);

return etot, e, f, v, vtot

end

function empiricalpotential(x, q, t, alist, neighlist, neighnum, eta, kappa, potential)

# compute full neighbor list
if !isempty(q)
    q = q[:,alist];
end

dim = size(x,1);
N = length(t)

e = zeros(N);
f = zeros(dim,N);
v = zeros(6,N);
for i = 1:length(potential)

    style, bondtype, bondedatomtypes, potentialfunc, mu, rcut = getpotential(potential[i]);
    potentialfunc = getfield(Main, Symbol(potentialfunc));

    # single potentials
    if (style == "single") & (bondtype == "nonbonded")
        ilist = Array(1:N);
        elseif (style == "single") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        ilist = findatomtype(Array(1:N), t, typei);
    end
    if (style == "single")
        xi, qi, ai, ti = neighsingles(x, q, t, ilist);
        ei, fi = potentialfunc(xi, qi, ti, mu, eta, kappa);
        ea, fa = tallysingle(ei, fi, ai);
        va = tallysinglevirial(fi, xi, ai, -1.0)
    end

    # pair potentials
    if (style == "pair") & (bondtype == "nonbonded")
        ilist = Int64.(Array(1:N));
        pairlist, pairnum = neighpairlist(x, ilist, neighlist, neighnum, rcut*rcut);
        elseif (style == "pair") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        typej = bondedatomtypes[2];

        ilist = findatomtype(Array(1:N), t, typei);
        pairlist, pairnum = neighpairlisttypej(x, ilist, alist, t, typej, neighlist, neighnum, rcut*rcut);
    end
    if (style == "pair")
        xij, ai, aj, ti, tj, qi, qj = neighpairs(x, q, pairlist, pairnum, t, ilist, alist);
        eij, fij = potentialfunc(xij, qi, qj, ti, tj, mu, eta, kappa);
        ea, fa = tallypair(eij, fij, ai, aj);
        va = tallypairvirial(fij, xij, ai, aj, -1.0/2)
    end

    # triplet potentials
    if (style == "triplet") & (bondtype == "nonbonded")
        ilist = Array(1:N);
        pairlist, pairnum = neighpairlist(x, ilist, neighlist, neighnum, rcut*rcut);
        elseif (style == "triplet") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        typej = bondedatomtypes[2];
        typek = bondedatomtypes[3];

        ilist = findatomtype(Array(1:N), t, typei);
        pairlist, pairnum = neighpairlisttypej(x, ilist, alist, t, [typej; typek], neighlist, neighnum, rcut*rcut);
    end
    if (style == "triplet")
        tripletlist, tripletnum = neightripletlist(pairlist, pairnum, ilist);
        xij, xik, ai, aj, ak, ti, tj, tk, qi, qj, qk = neightriplets(x, q, tripletlist, tripletnum, atomtype, ilist, alist);
        eij, fij, fik = potentialfunc(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
        ea, fa = tallytriplet(eij, fij, fik, ai, aj, ak);
        va = tallytripletvirial(fij, fik, xij, xik, ai, aj, ak, -1.0/3)
    end

    # quadlet potentials
    if (style == "quadlet") & (bondtype == "nonbonded")
        ilist = Array(1:N);
        pairlist, pairnum = neighpairlist(x, ilist, neighlist, neighnum, rcut*rcut);
        elseif (style == "quadlet") & (bondtype == "bonded")
        typei = bondedatomtypes[1];
        typej = bondedatomtypes[2];
        typek = bondedatomtypes[3];
        typel = bondedatomtypes[4];

        ilist = findatomtype(Array(1:N), t, typei);
        pairlist, pairnum = neighpairlisttypej(x, ilist, alist, t, [typej typek typel], neighlist, neighnum, rcut*rcut);
    end
    if (style == "quadlet")
        quadletlist, quadletnum = neighquadletlist(pairlist, pairnumsum, ilist);
        xij, xik, xil, ai, aj, ak, al, ti, tj, tk, tl, qi, qj, qk, ql =
                    neighquadlets(x, q, quadletlist, quadletnum, atomtype, ilist, alist);
        eij, fij, fik, fil = potentialfunc(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu, eta, kappa);
        ea, fa = tallyquadlet(eij, fij, fik, fil, ai, aj, ak, al);
        va = tallyquadletvirial(fij, fik, fil, xij, xik, xil, ai, aj, ak, al, -1.0/4)
    end

    e = e .+ ea[:];
    f = f .+ fa;
    v = v .+ va;
end

etot = sum(e);
vtot = sum(v,dims=2);

return etot, e, f, v, vtot

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP