text:
# This file is a part of BoostFractor.jl, licensed under the MIT License (MIT).

# Note: Travis doesn't seem to like getting no output for 10 mins.
#       Re-added status updates for CHEERLEADER and DANCER for this reason.

using Distributed

@everywhere begin
    using Test
    using FFTW
end

@testset "Package BoostFractor" begin

@everywhere begin
    PATH = ".."
    #PATH = "."#"BoostFractor.jl"
    # General Interface
    include(PATH*"/src/boost_fractor.jl")

    # Differnt algorithms
    include(PATH*"/src/boost_fractor_recurprop.jl")
    include(PATH*"/src/boost_fractor_transformer.jl")

    # Some convenient tools
    include(PATH*"/src/beam_gaussian.jl")
    include(PATH*"/src/beam_coupling.jl")
end

### Test CHEERLEADER ###
@everywhere begin
    # Initialize Coordinate System
    dx = 0.01
    coords = SeedCoordinateSystem(X = -0.5:dx:0.5, Y = -0.5:dx:0.5)

    diskR = 0.15

    sbdry = SeedSetupBoundaries(coords, diskno=1)

end

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# beam_gaussian.jl testing (100%)
#
# Needs to be behind cheerleader init to define coords!
gauss_shape = gauss_profile(coords; z = 1e-9, omega0 = 0.097, f = 22e9)
@test gauss_shape[4000] ≈ 0.0954822450586725 - 4.408488908455253e-8im # relatively central element, 1 index skims through all elements of 2d array

dx = 0.0001 * floor(10000. * (coords.X[2] - coords.X[1]))
dy = 0.0001 * floor(10000. * (coords.Y[2] - coords.Y[1]))
match = unnormalized_match(gauss_shape, 3 .* gauss_shape, dx=dx, dy=dy)
@test match ≈ 0.04433886791643955 - 5.80356296618764e-26im
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Let us sweep over different disk phase depths for a disk at distance lambda/2 infront of the mirror
# Of course the same way you can create a parallelized frequency sweep, etc.
@everywhere Deps = 0:1.:3*pi
@everywhere zcat(args...) = cat(dims = 3, args...)
@time result = @sync @distributed (zcat) for i in 1:length(Deps)
    deps = Deps[i]
    dnu = pi # One could optimize this for maximal boost in 1D before; its left out here for simplicitly and just set to pi
    print("Running dEps = $deps (dNu = $dnu)... \n")

    # Set the distances (the wavelength in free space is 30mm, in the medium 10mm (free space devided by sqrt(epsilon))
    sbdry.distance = [0.0,(dnu)*3e-2/(2*pi), (deps)*1e-2/(2*pi),0.e-3];

    nmax = 100 # Maximum Number of Iterations
    #1D (if you calculate in 1D, use (vcat) above instead of (zcat) )
    #cheerleader(0,nmax,sbdry; prop=propagator1D,f=10e9,Xset=X,Yset=Y)
    #3D
    cheerleader(0,nmax,sbdry,coords; prop=propagator,f=10e9,diskR=diskR)
end;

boost_factor_cheerleader = sum(abs2.(result), dims=(1,2))[1,1,:]*dx*dx/(pi*diskR^2);

println("Finished testing for CHEERLEADER")

#--------------------------------------------------------------------------------------------
### Test DANCER ###
@everywhere begin
    # Initialize Coordinate System
    dx = 0.01
    coords = SeedCoordinateSystem(X = -0.5:dx:0.5, Y = -0.5:dx:0.5)

    diskR = 0.15

    # SetupBoundaries
    # No Tilt or surface roughness here.
    epsilon = 9
    eps = Array{Complex{Float64}}([NaN,1,epsilon,1])
    distance = [0.0, 0.15, 0.15/sqrt.(epsilon), 0.0]*1e-2

    sbdry = SeedSetupBoundaries(coords, diskno=1, distance=distance, epsilon=eps)

end

# Let us sweep over different disk phase depths for a disk at distance lambda/2 infront of the mirror
# Of course the same way you can create a parallelized frequency sweep, etc.
@everywhere Deps = 0:1.:3*pi
@everywhere zcat(args...) = cat(dims = 3, args...)
@time result_dancer = @sync @distributed (zcat) for i in 1:length(Deps)
    deps = Deps[i]
    dnu = pi # One could optimize this for maximal boost in 1D before; its left out here for simplicitly and just set to pi
    print("Running dEps = $deps (dNu = $dnu)... \n")

    # Set the distances (the wavelength in free space is 30mm, in the medium 10mm (free space devided by sqrt(epsilon))
    sbdry.distance = [0.0, (dnu)*3e-2/(2*pi), Deps[i]*1e-2/(2*pi),0];

    nmax = 200 # Maximum Number of Iterations
    #1D (if you calculate in 1D, use (vcat) above instead of (zcat) )
    dancer(0,nmax,sbdry,coords; prop=propagator1D,f=10e9,diskR=diskR)
    #3D
    #dancer(0,nmax,sbdry,coords; prop=propagator,f=10e9,diskR=diskR)
end;

boost_factor_dancer = sum(abs2.(result_dancer), dims=(1,2))[1,1,:]*dx*dx/(pi*diskR^2);

println("Finished testing for DANCER")

#--------------------------------------------------------------------------------------------
### Test TRANSFORMER ###
@everywhere begin
    # Initialize Coordinate System
    dx = 0.01
    coords = SeedCoordinateSystem(X = -0.5:dx:0.5, Y = -0.5:dx:0.5)

    diskR = 0.15

    # SetupBoundaries
    # No Tilt or surface roughness here.
    epsilon = 9
    eps = Array{Complex{Float64}}([NaN,1,epsilon,1])
    distance = [0.0, 0.15, 0.15/sqrt.(epsilon), 0.0]*1e-2

    sbdry = SeedSetupBoundaries(coords, diskno=1, distance=distance, epsilon=eps)

    # Initialize modes
    Mmax = 10
    Lmax = 0 # l-Modes are irrelevant for the azimuthally symmetric haloscope
    # For a 1D calculation:
    modes_1D = SeedModes(coords, ThreeDim=false, Mmax=Mmax, Lmax=Lmax, diskR=diskR)
    # For 3D:
    modes = SeedModes(coords, ThreeDim=true, Mmax=Mmax, Lmax=Lmax, diskR=diskR)

end

# Let us sweep over different disk phase depths for a disk at distance lambda/2 infront of the mirror
# Of course the same way you can create a parallelized frequency sweep, etc.
@everywhere Deps = 0:1.:3*pi
@time result_transformer = @sync @distributed (hcat) for i in 1:length(Deps)
    deps = Deps[i]
    dnu = pi # One could optimize this for maximal boost in 1D before; its left out here for simplicitly and just set to pi
    print("Running dEps = $deps (dNu = $dnu)... \n")

    # Set the distances (the wavelength in free space is 30mm, in the medium 10mm (free space devided by sqrt(epsilon))
    sbdry.distance = [0.0, (dnu)*3e-2/(2*pi), Deps[i]*1e-2/(2*pi),0];

    #1D (if you calculate in 1D, use (vcat) above instead of (hcat) )
    #transformer(sbdry; prop=propagator1D,f=10e9,Xset=[1e-7],Yset=[1e-7])
    #3D
    transformer(sbdry,coords,modes; prop=propagator,f=10e9,diskR=diskR)
end;

boost_factor_transformer = sum(abs2.(result_transformer), dims=1)[1,:]

mode_coeffs = field2modes(nothing, coords, modes,diskR=diskR)
field = modes2field(mode_coeffs, coords, modes)
@test field[5000] ≈ 0.03647557026951998 + 0.0im

propmat = calc_propagation_matrices(sbdry, coords, modes)
@test propmat[1][4,6] ≈ 0.00021462613718063045 + 0.0im

println("Finished testing for TRANSFORMER")
#--------------------------------------------------------------------------------------------
@test boost_factor_cheerleader ≈ [0.9541048836881387, 4.182537182036321, 18.70293960502623, 7.465918522414759, 10.764938064402624, 15.415614620919474, 1.2667394643264536, 1.9715312870420645, 27.119751483329694, 8.413847982210477]
#@test boost_factor_cheerleader ≈ [0.9541048836881387, 0.9428833993720742, 0.968920462905896, 1.0341129917881549, 1.144486522206282, 1.311281404330431, 1.5530311134049792, 1.8991893743757546, 2.3962627706758717, 3.1180836339325797, 4.182537182036321, 5.77631829548155, 8.180541344734435, 11.750104248034686, 16.683053458978392, 22.328276414700372, 26.584187537779602, 27.43539899358354, 25.287916750692773, 21.957711497592978, 18.70293960502623, 15.968880656976728, 13.809281087444926, 12.13342947587664, 10.836197159631562, 9.835302811952268, 9.062953656331851, 8.469762044594315, 8.021400572162888, 7.692579478167976, 7.465918522414759, 7.329994000109685, 7.2782449520717565, 7.308492142283055, 7.422848236181687, 7.6280184632230785, 7.936389069538961, 8.36766853387498, 8.951220499791912, 9.729781130188297, 10.764938064402624, 12.144089632575064, 13.986113592200057, 16.433651816395216, 19.591343308612494, 23.3074109612804, 26.695661251021864, 27.883116341877454, 25.476205860505605, 20.569262201812844, 15.415614620919474, 11.236389682727756, 8.18288180649727, 6.037430281779245, 4.5362561494977, 3.473161524397452, 2.713994325825552, 2.16586121911677, 1.7663720636616334, 1.4758050703937462, 1.2667394643264536, 1.1207427057453578, 1.025753566948445, 0.9741434315936932, 0.961917112768637, 0.9883369519635989, 1.0560617558520191, 1.1717821857313868, 1.3476155646302996, 1.6035770542564973, 1.9715312870420645, 2.501814156824845, 3.2745053970180558, 4.417660301811236, 6.13329514001611, 8.719975781044885, 12.530192745139441, 17.668994151816275, 23.226061744057777, 26.92556913789544, 27.119751483329694, 24.648270973331165, 21.3056595232044, 18.163653151444, 15.554966421842678, 13.502429127893384, 11.911407897285025, 10.675141498985058, 9.71876066215671, 8.981342488510444, 8.413847982210477, 7.984858529840069, 7.671899347016415, 7.458149018169155, 7.33352793792804]

@test boost_factor_dancer ≈ [1.003029819130254, 5.497534861256151, 19.03228878235599, 7.817723902385752, 12.610984927623164, 14.178962970482115, 1.1434594478152733, 2.4841683203606983, 31.00826150988098, 8.548916138820806]
#@test boost_factor_dancer ≈ [0.9465597604600772, 0.9353352136853974, 0.9613500833206581, 1.0265038271281812, 1.1368236550259119, 1.3035506496401894, 1.5452177859585288, 1.8912763453371952, 2.3882277762141007, 3.1098957209819647, 4.174151879779258, 5.767671034944011, 8.171539715302092, 11.740624274034273, 16.6729565465068, 22.31747736268758, 26.572764341320273, 27.423602923437432, 25.276015311731925, 21.94587203360396, 18.69124362921158, 15.957338828020001, 13.797876691456961, 12.122159028174904, 10.825042594928846, 9.824242192983457, 9.05197733385444, 8.458858565994458, 8.010564389351762, 7.68180970131266, 7.4552156931162274, 7.3193608505429895, 7.267685228467108, 7.298011723964006, 7.412455350301969, 7.6177193601016056, 7.926182866326155, 8.357549109013213, 8.941179575288242, 9.719807990020836, 10.755018242265454, 12.134204765439094, 13.97624329007452, 16.423779229700262, 19.58146822920849, 23.29757595456879, 26.685987721853778, 27.873800005446203, 25.4673997532372, 20.56095299917797, 15.407657470121954, 11.228632064574768, 8.175203196725251, 6.029722950269289, 4.52845052200948, 3.465206267228134, 2.7058182554234054, 2.157406992392685, 1.7575956843424423, 1.4666424302967727, 1.2571295318833156, 1.1106325891076447, 1.015090899155427, 0.9628837252091647, 0.9500194106597534, 0.9757574022771622, 1.0427511871062949, 1.1576808370114224, 1.3326197879771289, 1.5875240918437874, 1.954222605311172, 2.48302081927242, 3.2539303359739677, 4.394881308723679, 6.107709566775495, 8.690766655560651, 12.496389375253083, 17.629819645108952, 23.18174731282909, 26.878049544460755, 27.071896250013232, 24.602197593707817, 21.262213904229338, 18.123030553046704, 15.516867175830251, 13.466334807311917, 11.877078481350779, 10.642339693624438, 9.687101909202692, 8.950609350269554, 8.383897489809803, 7.955480021177518, 7.6429588384112215, 7.4295446543075006, 7.30513529006866]

@test boost_factor_transformer ≈ [0.9219455134201273, 4.047480657662905, 18.261152610706294, 7.2468392586506525, 10.476628877197278, 15.022228781938008, 1.1827220364110227, 1.917381440805424, 26.539784031486857, 8.127709799564798]
#@test boost_factor_transformer ≈ [0.9219455134201293, 0.9088200235288822, 0.932472788773596, 0.9944207593641872, 1.1003581017820354, 1.2611863926386724, 1.4949283227109704, 1.830218075631676, 2.312297724310878, 3.013071038825134, 4.04748065766292, 5.597811746319159, 7.939063571323361, 11.419642965451624, 16.237339442286526, 21.762153491791135, 25.94078520367847, 26.792095591972686, 24.702468292477356, 21.447386690530784, 18.261152610706343, 15.58295201463342, 13.466455659627735, 11.823057780821372, 10.551035900830932, 9.569460134669937, 8.811943274009819, 8.230495827772176, 7.791034891609431, 7.468848293066695, 7.246839258650652, 7.113709251557361, 7.063016802684704, 7.092586196139127, 7.2042465260699515, 7.404537357496414, 7.705836236739639, 8.127673560768761, 8.698972696379933, 9.461765129881503, 10.47662887719731, 11.829593496954038, 13.637801537929464, 16.041996446339745, 19.145637191639622, 22.79913203746825, 26.127752327000355, 27.285415126839293, 24.90699832909787, 20.08049151233705, 15.022228781938072, 10.924340523719358, 7.934461008467611, 5.835365830636519, 4.365617413305234, 3.326589761075334, 2.585469801238749, 2.0507559250463068, 1.6628485969509776, 1.3823877426172526, 1.1827220364110236, 1.0458874408021794, 0.9597096314284665, 0.9162076426336727, 0.9108154369913783, 0.942226086594848, 1.012766421097211, 1.1291151754297613, 1.3034670791757523, 1.5557165291069917, 1.9173814408054293, 2.4382970793440837, 3.1977944804460456, 4.322700464678026, 6.013034962036725, 8.564165028005242, 12.323314613788607, 17.38715683051501, 22.83894840807172, 26.422258677978554, 26.539784031486974, 24.058135524654077, 20.74801236360884, 17.65094597322297, 15.093879494161175, 13.086093286518972, 11.528435507644087, 10.323654342883119, 9.394058173303977, 8.67688545686285, 8.127709799564812, 7.714600528320907, 7.414617399406159, 7.2120011630079555, 7.096533830554453]

end # testset


text:
function compute_reference_frame(body, lat, lon, alt::Float64)
    landing_position = SurfacePosition(body, lat, lon, body_frame)
    q_long = (
        0.0,
        sin(-lon * 0.5 * pi / 180),
        0.0,
        cos(-lon * 0.5 * pi / 180))
    q_lat = (
        0.0,
        0.0,
        sin(lat * 0.5 * pi / 180),
        cos(lat * 0.5 * pi / 180))
    cr1 = CreateRelative(body_frame, position=landing_position, rotation=q_long)
    cr2 = CreateRelative(cr1,rotation=q_lat)
    return CreateRelative(cr2,position=(alt,0.0,0.0))
end

function restrict_to_circle(radius, origin, vector)
    if norm(vector+origin) < radius
        return vector
    end
    normal = vector/norm(vector)
    inner = dot(normal, origin)^2 - norm(origin)^2 + radius^2
    if inner < 0
        # not intersecting, this is a weird case
        # we have to truncate the base vector to the sphere
        # shouldn't happen in real life, but needs to be handled
        return norm(origin)*radius - origin
    elseif inner == 0
        # bas must be touching the outside of the sphere
        # truncate the entire vector
        return vector*0
    elseif inner > 0
        # two solutions
        st = -dot(normal,origin)
        sol1 = st+sqrt(inner)
        sol2 = st-sqrt(inner)

        #we want the closest positive solution
        if sol1 < 0 && sol2 > 0
            out = sol2
        elseif sol1 > 0 && sol2 < 0
            out = sol1
        elseif sol1 < 0 && sol2 < 0 # they're both behind... somehow
            return norm(origin)*radius - origin
        elseif sol1 < sol2
            out = sol1
        else
            out = sol2
        end
        if out < norm(vector)
            return out * normal
        else
            return vector
        end
    end
end
function restrict_to_cone(basev, vect, conedir, theta)
    a = basev
    b = vect
    v = conedir
    lambdasq = cosd(theta)^2

    quad_a = dot(b,b)*(1-lambdasq) - dot(b,v)^2
    quad_b = 2*(dot(a,b)- lambdasq * dot(a,b) - dot(a,v) * dot(b,v))
    quad_c = dot(a,a) - lambdasq * dot(a,a) - dot(a,v)^2

    determinant = quad_b^2 - 4*quad_a*quad_c
    if determinant >= 0 && a != 0
        #intersection - need to find which is right
        #will always be the positive solution
        sol1 = (-quad_b + sqrt(determinant))/(2*quad_a)
        sol2 = (-quad_b - sqrt(determinant))/(2*quad_a)
        if sol1 >= 0
            return b*clamp(sol1,0,1)
        elseif sol2 >= 0
            return b*clamp(sol2,0,1)
        else
            return b #dunno
        end
    else
        # no intersections
        return vect
    end
end

function trajectory_predictor(itime::Float64, bodypos::SVector{3,Float64}, mu::Float64,
                              pos::SVector{3,Float64}, vel::SVector{3,Float64})
    steps = 30
    dt = itime/steps
    println(dt)
    cpos = pos
    cvel = vel
    for i=1:steps
        bodyvec = cpos - bodypos
        accel = -mu / (norm(bodyvec)^3) * bodyvec
        cvel = cvel + accel*dt
        cpos = cpos + cvel*dt
    end
    return cvel, cpos
end

ext: jl


text:
using SimpleDirectMediaLayer
using Test

SDL2_pkg_dir = joinpath(@__DIR__, "..","..")

@testset "Image" begin
    @test SDL_Init(SDL_INIT_VIDEO) == 0
    @test IMG_Load(joinpath(SDL2_pkg_dir, "assets", "cat.bmp")) != C_NULL
    @test IMG_Load(joinpath(SDL2_pkg_dir, "assets", "cat.png")) != C_NULL
    @test IMG_Load(joinpath(SDL2_pkg_dir, "assets", "cat.jpg")) != C_NULL
    SDL_Quit()
end


text:
function aimpoint(target::TargetRect)
    return (target.a/2,target.b/2)
end


text:
using DocOpt

doc="""Training (Single) Gaussian Mixture model for global variance.

Usage:
    train_gv.jl [options] <tgt_dir> <dst_jld>
    train_gv.jl --version
    train_gv.jl -h | --help

Options:
    -h --help           show this message
    --add_delta         add delta feature
    --n_components=MIX  number of mixtures [default: 1]
    --n_iter=ITER       number of iterations [default: 200]
    --n_init=N          number of initialization [default: 2]
    --max=MAX           maximum number of data we use [default: 100]
    --min_covar=MIN     minimum covariance [default: 1.0e-7]
"""

using VoiceConversion
using HDF5, JLD
using PyCall

@pyimport sklearn.mixture as mixture

using Logging
@Logging.configure(level=DEBUG, output=STDOUT)

let
    args = docopt(doc, version=v"0.0.1")

    nmax = parse(Int, args["--max"])
    add_delta = args["--add_delta"]

    n_components  = parse(Int, args["--n_components"])
    n_iter = parse(Int, args["--n_iter"])
    n_init = parse(Int, args["--n_init"])
    min_covar = parse(Float64, args["--min_covar"])

    dataset = GVDataset(args["<tgt_dir>"], add_delta=add_delta, nmax=nmax)

    gmm = mixture.GMM(n_components=n_components,
                      covariance_type="full",
                      n_iter=n_iter,
                      n_init=n_init,
                      min_covar=min_covar
                      )

    @show gmm

    # pass transposed matrix because python is row-major language
    elapsed = @elapsed gmm[:fit](dataset.X')
    info("Elapsed time in training is $(elapsed) sec.")

    # save transposed parameters because julia is column-major language
    # convert means
    py_means = gmm[:means_]
    means = py_means'

    # convert covar tensor
    py_covars = gmm[:covars_] # shape: (n_components, order, order)
    order = size(py_covars, 2)
    covars = Array(eltype(py_covars), order, order, n_components)
    for m=1:n_components
        covars[:,:,m] = reshape(py_covars[m,:,:], order, order)
    end

    save(args["<dst_jld>"],
         "description", "Parameters of Gaussian Mixture Model for GV",
         "n_components", n_components,
         "weights", gmm[:weights_],
         "means", means,
         "covars", covars
         )
end

main()


text:
# Auxiliar functions for RHF computations
"""
    function RHF_core_guess(ints::IntegralHelper)

Compute an initial orbital coefficient matrix using the core guess.
"""
function RHF_core_guess(ints::IntegralHelper)
    output("Using Core Guess")
    S = ints["S"]
    Λ = S^(-1/2)
    F = ints["T"] + ints["V"]
    Ft = Λ*F*Λ'

    # Get orbital energies and transformed coefficients
    _, Ct = LinearAlgebra.eigen(Symmetric(Ft), sortby=x->x)

    # Reverse transformation to get MO coefficients
    C = Λ*Ct

    return C, Λ
end

"""
    function RHF_gwh_guess(ints::IntegralHelper)

Compute an initial orbital coefficient matrix using the GWH guess.
"""
function RHF_gwh_guess(ints::IntegralHelper)

    # Form GWH guess
    output("Using GWH Guess")
    molecule = ints.molecule
    S = ints["S"]
    d, U = LinearAlgebra.eigen(Symmetric(S), sortby = x->1/abs(x))
    Λ = Hermitian(S)^(-1/2) |> Array
    idxs = [abs(d[i]) > 1E-7 for i = eachindex(d)]

    H = ints["T"] + ints["V"]
    ndocc = molecule.Nα
    nvir = size(S,1) - ndocc
    F = similar(S)

    for i = 1:(ndocc+nvir)
        F[i,i] = H[i,i]
        for j = i+1:ndocc+nvir
            F[i,j] = 0.875*S[i,j]*(H[i,i] + H[j,j])
            F[j,i] = F[i,j]
        end
    end
    Ft = Λ'*F*Λ

    # Get orbital energies and transformed coefficients
    _, Ct = LinearAlgebra.eigen(Symmetric(Ft), sortby=x->x)

    # Reverse transformation to get MO coefficients
    C = Λ*Ct

    return C, Λ
end

"""
    RHFEnergy(D::Array{Float64}, H::Array{Float64}, F::Array{Float64})

Compute RHF energy given a density matrix `D`, Core Hamiltonian `H` and Fock matrix `F`.
"""
function RHFEnergy(D::Array{Float64}, H::Array{Float64}, F::Array{Float64})
    return sum(D .* (H .+ F))
end

"""
    build_fock!(F::Array{Float64}, H::Array{Float64}, D::Array{Float64}, ERI::Array{Float64,4})

Build a Fock matrix into `F` using the Core Hamiltonian `H`, density matrix `D` and two-electron repulsion integral `ERI`.
"""
function build_fock!(F::Array{Float64}, H::Array{Float64}, D::Array{Float64}, ints::IntegralHelper{Float64,Chonky,AtomicOrbitals})
    ERI = ints["ERI"]
    F .= H
    @tensor F[m,n] += 2*D[r,s]*ERI[m,n,r,s]
    @tensor F[m,n] -= D[r,s]*ERI[m,r,n,s]
end

"""
    build_fock!(F::Array{Float64,2}, H::Array{Float64,2}, D::Array{Float64,2}, ints::IntegralHelper{Float64,<:AbstractDFERI,AtomicOrbitals})

Build a Fock matrix into `F` using the Core Hamiltonian `H`, density matrix `D` and two-electron repulsion integral `ERI`
approximated by density fitting.
"""
function build_fock!(F::Array{Float64,2}, H::Array{Float64,2}, D::Array{Float64,2}, ints::IntegralHelper{Float64,<:AbstractDFERI,AtomicOrbitals})
    b = ints["ERI"]
    F .= H
    @tensoropt F[m,n] += 2*D[r,s]*b[Q,m,n]*b[Q,r,s]
    @tensoropt F[m,n] -= D[r,s]*b[Q,m,r]*b[Q,n,s]
end

function build_fock!(F::Array{Float64}, H::Array{Float64}, D::Array{Float64}, ints::IntegralHelper{Float64,SparseERI,AtomicOrbitals})
    D = D
    F .= H
    eri_vals = ints["ERI"].data
    idxs = ints["ERI"].indexes

    nbas = ints.orbitals.basisset.nbas
    Farrays = [zeros(nbas, nbas) for i = 1:Threads.nthreads()]

    Threads.@threads for z = eachindex(eri_vals)
    @inbounds @fastmath begin
        i,j,k,l = idxs[z] .+ 1
        ν = eri_vals[z]
        Ft = Farrays[Threads.threadid()]
        ij = Fermi.index2(i-1,j-1)
        kl = Fermi.index2(k-1,l-1)

        # Logical auxiliar: γpq (whether p and q are different) Xpq = δpq + 1
        γij = i !== j
        γkl = k !== l
        γab = ij !== kl

        Xik = i === k ? 2.0 : 1.0
        Xjk = j === k ? 2.0 : 1.0
        Xil = i === l ? 2.0 : 1.0
        Xjl = j === l ? 2.0 : 1.0

        if γij && γkl && γab
            #J
            Ft[i,j] += 4.0*D[k,l]*ν
            Ft[k,l] += 4.0*D[i,j]*ν

            # K
            Ft[i,k] -= Xik*D[j,l]*ν
            Ft[j,k] -= Xjk*D[i,l]*ν
            Ft[i,l] -= Xil*D[j,k]*ν
            Ft[j,l] -= Xjl*D[i,k]*ν

        elseif γkl && γab
            # J
            Ft[i,j] += 4.0*D[k,l]*ν
            Ft[k,l] += 2.0*D[i,j]*ν

            # K
            Ft[i,k] -= Xik*D[j,l]*ν
            Ft[i,l] -= Xil*D[j,k]*ν
        elseif γij && γab
            # J
            Ft[i,j] += 2.0*D[k,l]*ν
            Ft[k,l] += 4.0*D[i,j]*ν

            # K
            Ft[i,k] -= Xik*D[j,l]*ν
            Ft[j,k] -= Xjk*D[i,l]*ν

        elseif γij && γkl

            # Only possible if i = k and j = l
            # and i < j ⇒ i < l

            # J
            Ft[i,j] += 4.0*D[k,l]*ν

            # K
            Ft[i,k] -= D[j,l]*ν
            Ft[i,l] -= D[j,k]*ν
            Ft[j,l] -= D[i,k]*ν
        elseif γab
            # J
            Ft[i,j] += 2.0*D[k,l]*ν
            Ft[k,l] += 2.0*D[i,j]*ν
            # K
            Ft[i,k] -= Xik*D[j,l]*ν
        else
            Ft[i,j] += 2.0*D[k,l]*ν
            Ft[i,k] -= D[j,l]*ν
        end
    end
    end

    Fred = sum(Farrays)
    Threads.@threads for i = 1:nbas
        @inbounds @fastmath begin
            F[i,i] += Fred[i,i]
            for j = (i+1):nbas
                F[i,j] += Fred[i,j] + Fred[j,i]
                F[j,i] = F[i,j]
            end
        end
    end
    F
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from DCM to quaternion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angle_to_quat

"""
    angle_to_quat(θ₁::T1, θ₂::T2, θ₃::T3, rot_seq::Symbol = :ZYX) where {T1<:Number, T2<:Number, T3<:Number}
    angle_to_quat(eulerang::EulerAngles)

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq` to a quaternion.

The input values of the origin Euler angles can also be passed inside the
structure `Θ` (see [`EulerAngles`](@ref)).

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

!!! note
    The type of the new quaternion will be obtained by promiting `T1`, `T2`, and
    `T3`.

# Remarks

This function assigns `q = q1 * q2 * q3` in which `qi` is the quaternion related
with the *i*-th rotation, `i Є [1,2,3]`.

# Example

```jldoctest
julia> angle_to_quat(pi / 2, pi / 3, pi / 4, :ZYX)
Quaternion{Float64}:
  + 0.701057 - 0.092296⋅i + 0.560986⋅j + 0.560986⋅k
```
"""
function angle_to_quat(
    θ₁::T1,
    θ₂::T2,
    θ₃::T3,
    rot_seq::Symbol = :ZYX
) where {T1<:Number, T2<:Number, T3<:Number}
    T = promote_type(T1, T2, T3)

    # Compute the sines and cosines of half angle.
    s₁, c₁ = sincos(T(θ₁) / 2)
    s₂, c₂ = sincos(T(θ₂) / 2)
    s₃, c₃ = sincos(T(θ₃) / 2)

    if rot_seq == :ZYX
        q0 = c₁ * c₂ * c₃ + s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * c₂ * s₃),
            s * (s₁ * c₂ * c₃ - c₁ * s₂ * s₃)
        )
    elseif rot_seq == :XYX
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (s₁ * s₂ * c₃ - c₁ * s₂ * s₃)
        )
    elseif rot_seq == :XYZ
        q0 = c₁ * c₂ * c₃ - s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (s₁ * c₂ * c₃ + c₁ * s₂ * s₃),
            s * (c₁ * s₂ * c₃ - s₁ * c₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * s₂ * c₃)
        )
    elseif rot_seq == :XZX
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃)
        )
    elseif rot_seq == :XZY
        q0 = c₁ * c₂ * c₃ + s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (s₁ * c₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * c₂ * s₃)
        )
    elseif rot_seq == :YXY
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * s₃ - s₁ * s₂ * c₃)
        )
    elseif rot_seq == :YXZ
        q0 = c₁ * c₂ * c₃ + s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ + s₁ * c₂ * s₃),
            s * (s₁ * c₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ - s₁ * s₂ * c₃)
        )
    elseif rot_seq == :YZX
        q0 = c₁ * c₂ * c₃ - s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ + s₁ * s₂ * c₃),
            s * (s₁ * c₂ * c₃ + c₁ * s₂ * s₃),
            s * (c₁ * s₂ * c₃ - s₁ * c₂ * s₃)
        )
    elseif rot_seq == :YZY
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (s₁ * s₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃)
        )
    elseif rot_seq == :ZXY
        q0 = c₁ * c₂ * c₃ - s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ - s₁ * c₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * s₂ * c₃),
            s * (s₁ * c₂ * c₃ + c₁ * s₂ * s₃)
        )
    elseif rot_seq == :ZXZ
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (s₁ * s₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃)
        )
    elseif rot_seq == :ZYZ
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃)
        )
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

angle_to_quat(Θ::EulerAngles) = angle_to_quat(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)


text:
using Test, WorkerUtilities

@testset "WorkerUtilities" begin

    @testset "WorkerUtilities.@spawn" begin

        WorkerUtilities.init()
        threadid = fetch(WorkerUtilities.@spawn(Threads.threadid()))
        @test Threads.nthreads() == 1 ? (threadid == 1) : (threadid != 1)
        @test WorkerUtilities.@spawn(false, 1 + 1).storage === nothing

    end # @testset "WorkerUtilities.@spawn"

    @testset "Lockable" begin
        # Lockable{T, L<:AbstractLock}
        let # test the constructor `Lockable(value, lock)`
            lockable = Lockable(Dict("foo" => "hello"), ReentrantLock())
            @test lockable.value["foo"] == "hello"
            lock(lockable) do d
                @test d["foo"] == "hello"
            end
            lock(lockable) do d
                d["foo"] = "goodbye"
            end
            @test lockable.value["foo"] == "goodbye"
            lock(lockable) do d
                @test d["foo"] == "goodbye"
            end
        end
        let # test the constructor `Lockable(value)`
            lockable = Lockable(Dict("foo" => "hello"))
            @test lockable.value["foo"] == "hello"
            lock(lockable) do d
                @test d["foo"] == "hello"
            end
            lock(lockable) do d
                d["foo"] = "goodbye"
            end
            @test lockable.value["foo"] == "goodbye"
            lock(lockable) do d
                @test d["foo"] == "goodbye"
            end
        end
    end # @testset "Lockable"

end # @testset "WorkerUtilities"


text:
abstract type CSVWriter end

write_csv_row!(io::IO, vec::AbstractVector) = writedlm(io, reshape(vec, 1, :), ",")
write_csv_row!(io::IO, array::AbstractArray) = write_csv_row!(io, vec(array))

function write_inputs_csv(writer::CSVWriter, dir::String, robot_name::String, samples::AbstractVector{<:Sample})
    open(joinpath(dir, "$(robot_name)_inputs.csv"), "w") do io
        header = make_header(writer)
        write_csv_row!(io, header)
        for sample in samples
            sampledata = make_sample_data(writer, sample)
            write_csv_row!(io, vcat(sampledata...))
        end
    end
end

function write_expected_result_csvs(writer::CSVWriter, dir::String, robot_name::String, sample::Sample)
    Jlib = velocity_jacobian(writer, sample)

    state = MechanismState(writer.mechanism)
    set_configuration!(state, sample.q)
    set_velocity!(state, sample.v)
    M = mass_matrix(state)
    τid = inverse_dynamics(state, sample.v̇)
    c = dynamics_bias(state)

    Mlib = (Jlib' \ M) * Jlib'
    # τidlib = Jlib' \ τid
    # τlib = Jlib' \ sample.τ
    # clib = Jlib' \ c
    _, _, v̇lib, τlib = make_sample_data(writer, sample)

    open(joinpath(dir, "$(robot_name)_inverse_dynamics_expected.csv"), "w") do io
        write_csv_row!(io, τlib)
    end

    open(joinpath(dir, "$(robot_name)_mass_matrix_expected.csv"), "w") do io
        write_csv_row!(io, Mlib)
    end

    open(joinpath(dir, "$(robot_name)_dynamics_expected.csv"), "w") do io
        write_csv_row!(io, v̇lib)
    end
end

function make_header(writer::CSVWriter)
    header = String[]
    append!(header, configuration_coordinate_names(writer))
    append!(header, tangential_coordinate_names(writer, "v"))
    append!(header, tangential_coordinate_names(writer, "vd"))
    append!(header, tangential_coordinate_names(writer, "tau"))
    header
end

function configuration_coordinate_names(joint::Joint{<:Any, <:OneDegreeOfFreedomFixedAxis}, ::CSVWriter)
    ["q_$(string(joint))"]
end

function tangential_coordinate_names(joint::Joint{<:Any, <:OneDegreeOfFreedomFixedAxis}, prefix::String, ::CSVWriter)
    ["$(prefix)_$(string(joint))"]
end

function velocity_jacobian(writer::CSVWriter, sample::Sample)
    ForwardDiff.jacobian(sample.v) do v
        T = eltype(v)
        sample_dual = Sample(copyto!(similar(sample.q, T), sample.q), v, copyto!(similar(sample.v̇, T), sample.v̇), copyto!(similar(sample.τ, T), sample.τ))
        _, v_lib, _, _ = make_sample_data(writer, sample_dual)
        v_lib
    end
end

# RigidBodyDynamics.jl
struct RigidBodyDynamicsJlWriter <:CSVWriter
    mechanism::Mechanism{Float64}
end

RigidBodyDynamicsJlWriter(mechanism::Mechanism, robot_name::String) = RigidBodyDynamicsJlWriter(mechanism)

shortlibname(::RigidBodyDynamicsJlWriter) = "rbdjl"

function configuration_coordinate_names(writer::RigidBodyDynamicsJlWriter)
    configuration_coordinate_names.(tree_joints(writer.mechanism), Ref(writer)) |> flatten |> collect
end

function tangential_coordinate_names(writer::RigidBodyDynamicsJlWriter, prefix::String)
    tangential_coordinate_names.(tree_joints(writer.mechanism), Ref(prefix), Ref(writer)) |> flatten |> collect
end

function make_sample_data(writer::RigidBodyDynamicsJlWriter, sample::Sample)
    (sample.q, sample.v, sample.v̇, sample.τ)
end

function configuration_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, ::RigidBodyDynamicsJlWriter)
    jointprefix = "q_$(string(joint))"
    [
        jointprefix * "_qs",
        jointprefix * "_qx",
        jointprefix * "_qy",
        jointprefix * "_qz",
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
    ]
end

function tangential_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, prefix::String, ::RigidBodyDynamicsJlWriter)
    jointprefix = "$(prefix)_$(string(joint))"
    [
        jointprefix * "_wx",
        jointprefix * "_wy",
        jointprefix * "_wz",
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
    ]
end

# RBDL
struct RBDLWriter{S<:StateCache} <: CSVWriter
    mechanism::Mechanism
    statecache::S
    joints::Vector{Joint{Float64}}
    q_permutation::Vector{Int}
end

const RBDL_JOINT_NAMES = Dict{String, Vector{String}}()
RBDL_JOINT_NAMES["iiwa"] = [
    "lbr_iiwa_joint_1",
    "lbr_iiwa_joint_2",
    "lbr_iiwa_joint_3",
    "lbr_iiwa_joint_4",
    "lbr_iiwa_joint_5",
    "lbr_iiwa_joint_6",
    "lbr_iiwa_joint_7"
]
RBDL_JOINT_NAMES["hyq"] = [
    "base_link_0_to_world",
    "lf_haa_joint",
    "lf_hfe_joint",
    "lf_kfe_joint",
    "lh_haa_joint",
    "lh_hfe_joint",
    "lh_kfe_joint",
    "rf_haa_joint",
    "rf_hfe_joint",
    "rf_kfe_joint",
    "rh_haa_joint",
    "rh_hfe_joint",
    "rh_kfe_joint"
]
RBDL_JOINT_NAMES["atlas"] = [
    "pelvis_to_world",
    "back_bkz",
    "back_bky",
    "back_bkx",
    "l_arm_shz",
    "l_arm_shx",
    "l_arm_ely",
    "l_arm_elx",
    "l_arm_uwy",
    "l_arm_mwx",
    "l_arm_lwy",
    "neck_ay",
    "r_arm_shz",
    "r_arm_shx",
    "r_arm_ely",
    "r_arm_elx",
    "r_arm_uwy",
    "r_arm_mwx",
    "r_arm_lwy",
    "l_leg_hpz",
    "l_leg_hpx",
    "l_leg_hpy",
    "l_leg_kny",
    "l_leg_aky",
    "l_leg_akx",
    "r_leg_hpz",
    "r_leg_hpx",
    "r_leg_hpy",
    "r_leg_kny",
    "r_leg_aky",
    "r_leg_akx"
]

function RBDLWriter(mechanism::Mechanism, robot_name::String)
    # q order is permuted so that quaternion scalar components come last.
    # https://bitbucket.org/rbdl/rbdl/src/0879ee8c548a07b8a9e1296bcbb5b006e1a098c7/src/Model.cc#lines-374:375
    joint_names = RBDL_JOINT_NAMES[robot_name]
    joints_rbdl = Vector{Joint{Float64}}(undef, length(joints(mechanism)))
    for i in eachindex(joints_rbdl)
        joints_rbdl[i] = findjoint(mechanism, joint_names[i])
    end
    @assert isempty(setdiff(joints_rbdl, tree_joints(mechanism)))
    state = MechanismState(mechanism)
    q_permutation = Int[]
    quatjoints = Joint[]
    for joint in joints_rbdl
        qrange = configuration_range(state, joint)
        vrange = velocity_range(state, joint)
        if joint_type(joint) isa OneDegreeOfFreedomFixedAxis
            append!(q_permutation, qrange)
        elseif joint_type(joint) isa QuaternionFloating
            push!(quatjoints, joint)
            append!(q_permutation, qrange[5 : 7])
            append!(q_permutation, qrange[2 : 4])
        else
            error("Joint type not supported")
        end
    end
    for joint in quatjoints
        push!(q_permutation, configuration_range(state, joint)[1])
    end
    RBDLWriter(mechanism, StateCache(mechanism), joints_rbdl, q_permutation)
end

shortlibname(::RBDLWriter) = "rbdl"

function configuration_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, ::RBDLWriter)
    # NOTE: still transformed by q_permutation after (easier that way)
    jointprefix = "q_$(string(joint))"
    [
        jointprefix * "_qs",
        jointprefix * "_qx",
        jointprefix * "_qy",
        jointprefix * "_qz",
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
    ]
end

function tangential_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, prefix::String, ::RBDLWriter)
    jointprefix = "$(prefix)_$(string(joint))"
    [
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
        jointprefix * "_wx",
        jointprefix * "_wy",
        jointprefix * "_wz"
    ]
end

function configuration_coordinate_names(writer::RBDLWriter)
    names = configuration_coordinate_names.(tree_joints(writer.mechanism), Ref(writer)) |> flatten |> collect
    names[writer.q_permutation]
end

function tangential_coordinate_names(writer::RBDLWriter, prefix::String)
    tangential_coordinate_names.(writer.joints, Ref(prefix), Ref(writer)) |> flatten |> collect
end

function make_sample_data(writer::RBDLWriter, sample::Sample{T}) where T
    mechanism = writer.mechanism
    state = writer.statecache[T]
    rootframe = root_frame(mechanism)
    set_configuration!(state, sample.q)
    set_velocity!(state, sample.v)

    # q
    q_rbdl = sample.q[writer.q_permutation]

    # v, v̇, τ
    v_rbdl = fill!(similar(velocity(state)), NaN)
    v̇_rbdl = fill!(similar(velocity(state)), NaN)
    τ_rbdl = fill!(similar(velocity(state)), NaN)
    vstart = 1
    for joint in writer.joints
        jointid = JointID(joint)
        vrange = velocity_range(state, jointid)
        if joint_type(joint) isa OneDegreeOfFreedomFixedAxis
            v_rbdl[vstart : vstart] .= sample.v[vrange]
            v̇_rbdl[vstart : vstart] .= sample.v̇[vrange]
            τ_rbdl[vstart : vstart] .= sample.τ[vrange]
        elseif joint_type(joint) isa QuaternionFloating
            # RBDL rewrites
            #   [world] -- floating joint -- [floating body]
            # as
            #   [world] -- translation xyz -- [null body] -- spherical -- [floating body]
            # The velocity for an RBDL translation xyz joint is just the time derivative of the position.
            # The velocity for an RBDL spherical joint is the angular velocity in body frame
            # (https://bitbucket.org/rbdl/rbdl/src/0879ee8c548a07b8a9e1296bcbb5b006e1a098c7/include/rbdl/Joint.h#lines-135:144).

            v̇joint = sample.v̇[jointid]
            τjoint = sample.τ[jointid]
            jointtf = joint_transform(state, joint)
            jointtwist = RigidBodyDynamics.twist(state, joint)
            jointaccel = SpatialAcceleration(jointtwist.body, jointtwist.base, jointtwist.frame, SVector{3}(v̇joint[1 : 3]), SVector{3}(v̇joint[4 : 6]))
            jointwrench = Wrench(frame_after(joint), SVector{3}(τjoint[1 : 3]), SVector{3}(τjoint[4 : 6]))

            body_origin = Point3D(jointtwist.body, zero(T), zero(T), zero(T))
            body_origin_velocity = transform(point_velocity(jointtwist, body_origin), jointtf)
            body_origin_accel = transform(point_acceleration(jointtwist, jointaccel, body_origin), jointtf)

            # v
            v_rbdl[vstart .+ (0 : 2)] .= body_origin_velocity.v
            v_rbdl[vstart .+ (3 : 5)] .= angular(jointtwist)

            # v̇
            v̇_rbdl[vstart .+ (0 : 2)] .= body_origin_accel.v
            v̇_rbdl[vstart .+ (3 : 5)] .= angular(jointaccel)

            # τ
            τ_rbdl[vstart .+ (0 : 2)] .= linear(transform(jointwrench, jointtf))
            τ_rbdl[vstart .+ (3 : 5)] .= angular(jointwrench)
        else
            error("Joint type not yet supported")
        end
        vstart += num_velocities(joint)
    end

    (q_rbdl, v_rbdl, v̇_rbdl, τ_rbdl)
end

# Pinocchio
struct PinocchioWriter <: CSVWriter
    mechanism::Mechanism{Float64}
    q_permutation::Vector{Int}
    v_permutation::Vector{Int}
end

const PINOCCHIO_JOINT_NAMES = Dict{String, Vector{String}}()
PINOCCHIO_JOINT_NAMES["iiwa"] = [
    "lbr_iiwa_joint_1",
    "lbr_iiwa_joint_2",
    "lbr_iiwa_joint_3",
    "lbr_iiwa_joint_4",
    "lbr_iiwa_joint_5",
    "lbr_iiwa_joint_6",
    "lbr_iiwa_joint_7"
]
PINOCCHIO_JOINT_NAMES["hyq"] = [
    "base_link_0_to_world",
    "lf_haa_joint",
    "lf_hfe_joint",
    "lf_kfe_joint",
    "lh_haa_joint",
    "lh_hfe_joint",
    "lh_kfe_joint",
    "rf_haa_joint",
    "rf_hfe_joint",
    "rf_kfe_joint",
    "rh_haa_joint",
    "rh_hfe_joint",
    "rh_kfe_joint"
]
PINOCCHIO_JOINT_NAMES["atlas"] = [
    "pelvis_to_world",
    "back_bkz",
    "back_bky",
    "back_bkx",
    "l_arm_shz",
    "l_arm_shx",
    "l_arm_ely",
    "l_arm_elx",
    "l_arm_uwy",
    "l_arm_mwx",
    "l_arm_lwy",
    "neck_ay",
    "r_arm_shz",
    "r_arm_shx",
    "r_arm_ely",
    "r_arm_elx",
    "r_arm_uwy",
    "r_arm_mwx",
    "r_arm_lwy",
    "l_leg_hpz",
    "l_leg_hpx",
    "l_leg_hpy",
    "l_leg_kny",
    "l_leg_aky",
    "l_leg_akx",
    "r_leg_hpz",
    "r_leg_hpx",
    "r_leg_hpy",
    "r_leg_kny",
    "r_leg_aky",
    "r_leg_akx"
]

function PinocchioWriter(mechanism::Mechanism, robot_name::String)
    joint_names = PINOCCHIO_JOINT_NAMES[robot_name]
    joints_pinocchio = findjoint.(Ref(mechanism), joint_names)
    @assert isempty(setdiff(joints_pinocchio, tree_joints(mechanism)))
    state = MechanismState(mechanism)
    q_permutation = Int[]
    v_permutation = Int[]
    for joint in joints_pinocchio
        qrange = configuration_range(state, joint)
        vrange = velocity_range(state, joint)
        if joint_type(joint) isa OneDegreeOfFreedomFixedAxis
            append!(q_permutation, qrange)
            append!(v_permutation, vrange)
        elseif joint_type(joint) isa QuaternionFloating
            append!(q_permutation, qrange[[5, 6, 7, 2, 3, 4, 1]])
            append!(v_permutation, vrange[[4, 5, 6, 1, 2, 3]])
        else
            error("Joint type not supported")
        end
    end
    @assert length(q_permutation) == num_positions(mechanism)
    @assert length(v_permutation) == num_velocities(mechanism)
    PinocchioWriter(mechanism, q_permutation, v_permutation)
end

shortlibname(::PinocchioWriter) = "pinocchio"

function configuration_coordinate_names(writer::PinocchioWriter)
    names = configuration_coordinate_names.(tree_joints(writer.mechanism), Ref(writer)) |> flatten |> collect
    names[writer.q_permutation]
end

function tangential_coordinate_names(writer::PinocchioWriter, prefix::String)
    names = tangential_coordinate_names.(tree_joints(writer.mechanism), Ref(prefix), Ref(writer)) |> flatten |> collect
    names[writer.v_permutation]
end

function make_sample_data(writer::PinocchioWriter, sample::Sample)
    q_pinocchio = sample.q[writer.q_permutation]
    v_pinocchio = sample.v[writer.v_permutation]
    v̇_pinocchio = sample.v̇[writer.v_permutation]
    τ_pinocchio = sample.τ[writer.v_permutation]
    (q_pinocchio, v_pinocchio, v̇_pinocchio, τ_pinocchio)
end

function configuration_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, ::PinocchioWriter)
    # NOTE: still transformed by q_permutation after (easier that way)
    jointprefix = "q_$(string(joint))"
    [
        jointprefix * "_qs",
        jointprefix * "_qx",
        jointprefix * "_qy",
        jointprefix * "_qz",
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
    ]
end

function tangential_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, prefix::String, ::PinocchioWriter)
    # NOTE: still transformed by v_permutation after (easier that way)
    jointprefix = "$(prefix)_$(string(joint))"
    [
        jointprefix * "_wx",
        jointprefix * "_wy",
        jointprefix * "_wz",
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
    ]
end

# RobCoGen
struct RobCoGenWriter <:CSVWriter
    mechanism::Mechanism{Float64}
    non_floating_joints::Vector{Joint{Float64}}
    non_floating_q_indices::Vector{Int}

    function RobCoGenWriter(mechanism::Mechanism{Float64}, robot_name::String)
        non_floating_joints = Joint{Float64}[j for j in tree_joints(mechanism) if !isfloating(j)]
        state = MechanismState(mechanism)
        non_floating_q_indices = [configuration_range(state, j) for j in non_floating_joints] |> flatten |> collect
        new(mechanism, non_floating_joints, non_floating_q_indices)
    end
end

shortlibname(::RobCoGenWriter) = "robcogen"

function tangential_coordinate_names(joint::Joint{<:Any, <:QuaternionFloating}, prefix::String, ::RobCoGenWriter)
    jointprefix = "$(prefix)_$(string(joint))"
    [
        jointprefix * "_wx",
        jointprefix * "_wy",
        jointprefix * "_wz",
        jointprefix * "_x",
        jointprefix * "_y",
        jointprefix * "_z",
    ]
end

function configuration_coordinate_names(writer::RobCoGenWriter)
    configuration_coordinate_names.(writer.non_floating_joints, Ref(writer)) |> flatten |> collect
end

function tangential_coordinate_names(writer::RobCoGenWriter, prefix::String)
    tangential_coordinate_names.(tree_joints(writer.mechanism), Ref(prefix), Ref(writer)) |> flatten |> collect
end

function make_sample_data(writer::RobCoGenWriter, sample::Sample)
    (sample.q[writer.non_floating_q_indices], sample.v, sample.v̇, sample.τ)
end

max_stars_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_issues_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_forks_repo_name: CobbledSteel/rbd-benchmarks-riscv


text:
"""
    errorbars(x, y, error_both; kwargs...)
    errorbars(x, y, error_low, error_high; kwargs...)
    errorbars(x, y, error_low_high; kwargs...)

    errorbars(xy, error_both; kwargs...)
    errorbars(xy, error_low, error_high; kwargs...)
    errorbars(xy, error_low_high; kwargs...)

    errorbars(xy_error_both; kwargs...)
    errorbars(xy_error_low_high; kwargs...)

Plots errorbars at xy positions, extending by errors in the given `direction`.

If you want to plot intervals from low to high values instead of relative errors, use `rangebars`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Errorbars) do scene
    Theme(
        whiskerwidth = 0,
        color = theme(scene, :linecolor),
        linewidth = theme(scene, :linewidth),
        direction = :y,
        visible = theme(scene, :visible),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        inspectable = theme(scene, :inspectable),
        transparency = false
    )
end

"""
    rangebars(val, low, high; kwargs...)
    rangebars(val, low_high; kwargs...)
    rangebars(val_low_high; kwargs...)

Plots rangebars at `val` in one dimension, extending from `low` to `high` in the other dimension
given the chosen `direction`.

If you want to plot errors relative to a reference value, use `errorbars`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Rangebars) do scene
    Theme(
        whiskerwidth = 0,
        color = theme(scene, :linecolor),
        linewidth = theme(scene, :linewidth),
        direction = :y,
        visible = theme(scene, :visible),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        inspectable = theme(scene, :inspectable),
        transparency = false
    )
end

### conversions for errorbars

function Makie.convert_arguments(::Type{<:Errorbars}, x, y, error_both)
    xyerr = broadcast(x, y, error_both) do x, y, e
        Vec4f(x, y, e, e)
    end
    (xyerr,)
end

function Makie.convert_arguments(::Type{<:Errorbars}, x, y, error_low, error_high)
    xyerr = broadcast(Vec4f, x, y, error_low, error_high)
    (xyerr,)
end

function Makie.convert_arguments(::Type{<:Errorbars}, x, y, error_low_high::AbstractVector{<:VecTypes{2}})
    xyerr = broadcast(x, y, error_low_high) do x, y, (el, eh)
        Vec4f(x, y, el, eh)
    end
    (xyerr,)
end

function Makie.convert_arguments(::Type{<:Errorbars}, xy::AbstractVector{<:VecTypes{2}}, error_both)
    xyerr = broadcast(xy, error_both) do (x, y), e
        Vec4f(x, y, e, e)
    end
    (xyerr,)
end

function Makie.convert_arguments(::Type{<:Errorbars}, xy::AbstractVector{<:VecTypes{2}}, error_low, error_high)
    xyerr = broadcast(xy, error_low, error_high) do (x, y), el, eh
        Vec4f(x, y, el, eh)
    end
    (xyerr,)
end

function Makie.convert_arguments(::Type{<:Errorbars}, xy::AbstractVector{<:VecTypes{2}}, error_low_high::AbstractVector{<:VecTypes{2}})
    xyerr = broadcast(xy, error_low_high) do (x, y), (el, eh)
        Vec4f(x, y, el, eh)
    end
    (xyerr,)
end

function Makie.convert_arguments(::Type{<:Errorbars}, xy_error_both::AbstractVector{<:VecTypes{3}})
    xyerr = broadcast(xy_error_both) do (x, y, e)
        Vec4f(x, y, e, e)
    end
    (xyerr,)
end

### conversions for rangebars

function Makie.convert_arguments(::Type{<:Rangebars}, val, low, high)
    val_low_high = broadcast(Vec3f, val, low, high)
    (val_low_high,)
end

function Makie.convert_arguments(::Type{<:Rangebars}, val, low_high)
    val_low_high = broadcast(val, low_high) do val, (low, high)
        Vec3f(val, low, high)
    end
    (val_low_high,)
end

### the two plotting functions create linesegpairs in two different ways
### and then hit the same underlying implementation in `_plot_bars!`

function Makie.plot!(plot::Errorbars{T}) where T <: Tuple{AbstractVector{<:VecTypes{4}}}

    x_y_low_high = plot[1]

    is_in_y_direction = lift(plot.direction) do dir
        if dir == :y
            true
        elseif dir == :x
            false
        else
            error("Invalid direction $dir. Options are :x and :y.")
        end
    end

    linesegpairs = lift(x_y_low_high, is_in_y_direction) do x_y_low_high, in_y

        map(x_y_low_high) do (x, y, l, h)
            in_y ?
                (Point2f(x, y - l), Point2f(x, y + h)) :
                (Point2f(x - l, y), Point2f(x + h, y))
        end
    end

    _plot_bars!(plot, linesegpairs, is_in_y_direction)
end

function Makie.plot!(plot::Rangebars{T}) where T <: Tuple{AbstractVector{<:VecTypes{3}}}

    val_low_high = plot[1]

    is_in_y_direction = lift(plot.direction) do dir
        if dir == :y
            true
        elseif dir == :x
            false
        else
            error("Invalid direction $dir. Options are :x and :y.")
        end
    end

    linesegpairs = lift(val_low_high, is_in_y_direction) do vlh, in_y

        map(vlh) do (v, l, h)
            in_y ?
                (Point2f(v, l), Point2f(v, h)) :
                (Point2f(l, v), Point2f(h, v))
        end
    end

    _plot_bars!(plot, linesegpairs, is_in_y_direction)
end

function _plot_bars!(plot, linesegpairs, is_in_y_direction)

    f_if(condition, f, arg) = condition ? f(arg) : arg

    @extract plot (whiskerwidth, color, linewidth, visible, colormap, colorrange, inspectable, transparency)

    scene = parent_scene(plot)

    whiskers = lift(linesegpairs, scene.camera.projectionview,
        scene.camera.pixel_space, whiskerwidth) do pairs, _, _, whiskerwidth

        endpoints = [p for pair in pairs for p in pair]

        screenendpoints = scene_to_screen(endpoints, scene)

        screenendpoints_shifted_pairs = map(screenendpoints) do sep
            (sep .+ f_if(is_in_y_direction[], reverse, Point(0, -whiskerwidth/2)),
             sep .+ f_if(is_in_y_direction[], reverse, Point(0,  whiskerwidth/2)))
        end

        screen_to_scene([p for pair in screenendpoints_shifted_pairs for p in pair], scene)
    end

    whiskercolors = lift(Any, color) do color
        # we have twice as many linesegments for whiskers as we have errorbars, so we
        # need to duplicate colors if a vector of colors is given
        if color isa AbstractVector
            repeat(color, inner = 2)
        else
            color
        end
    end

    whiskerlinewidths = lift(Any, linewidth) do linewidth
        # same for linewidth
        if linewidth isa AbstractVector
            repeat(linewidth, inner = 2)
        else
            linewidth
        end
    end

    linesegments!(
        plot, linesegpairs, color = color, linewidth = linewidth, visible = visible,
        colormap = colormap, colorrange = colorrange, inspectable = inspectable,
        transparency = transparency
    )
    linesegments!(
        plot, whiskers, color = whiskercolors, linewidth = whiskerlinewidths,
        visible = visible, colormap = colormap, colorrange = colorrange,
        inspectable = inspectable, transparency = transparency
    )
    plot
end

function scene_to_screen(pts, scene)
    p4 = to_ndim.(Vec4f, to_ndim.(Vec3f, pts, 0.0), 1.0)
    p1m1 = Ref(scene.camera.projectionview[]) .* p4
    projected = Ref(inv(scene.camera.pixel_space[])) .* p1m1
    [Point2.(p[1:2]...) for p in projected]
end

function screen_to_scene(pts, scene)
    p4 = to_ndim.(Vec4f, to_ndim.(Vec3f, pts, 0.0), 1.0)
    p1m1 = Ref(scene.camera.pixel_space[]) .* p4
    projected = Ref(inv(scene.camera.projectionview[])) .* p1m1
    [Point2.(p[1:2]...) for p in projected]
end

function scene_to_screen(p::T, scene) where T <: Point
    p4 = to_ndim(Vec4f, to_ndim(Vec3f, p, 0.0), 1.0)
    p1m1 = scene.camera.projectionview[] * p4
    projected = inv(scene.camera.pixel_space[]) * p1m1
    T(projected[1:2]...)
end

function screen_to_scene(p::T, scene) where T <: Point
    p4 = to_ndim(Vec4f, to_ndim(Vec3f, p, 0.0), 1.0)
    p1m1 = scene.camera.pixel_space[] * p4
    projected = inv(scene.camera.projectionview[]) * p1m1
    T(projected[1:2]...)
end

# ignore whiskers when determining data limits
function data_limits(bars::Union{Errorbars, Rangebars})
    data_limits(bars.plots[1])
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b36fb690-3447-11eb-3629-edbddae4ebd1
input = readlines("./input.txt");

# ╔═╡ 3d1fca00-3449-11eb-3e0e-5784380ebdbc
expenseReport = map(x -> parse(Int32,x),input)

# ╔═╡ 4ab2f908-3449-11eb-299d-53d5a0a3044e
len = size(expenseReport)[1]

# ╔═╡ 473a9452-3449-11eb-1dd1-25d770220737
let p1 = nothing
	for i ∈ 1:len, j ∈ i+1:len
		if expenseReport[j]+ expenseReport[i] == 2020
			p1= expenseReport[i]*expenseReport[j]
			break
		end
	end
	p1
end

# ╔═╡ 6e67d06a-344b-11eb-15c6-ed3cfcfd8a83
let p2 = nothing
	for i ∈ 1:len, j ∈ i+1:len
		if (2020-(expenseReport[j]+ expenseReport[i])) in expenseReport
			p2= expenseReport[i]*expenseReport[j]*(2020 - (expenseReport[j] + expenseReport[i] ))
			break
		end
	end
	p2
end

# ╔═╡ Cell order:
# ╠═b36fb690-3447-11eb-3629-edbddae4ebd1
# ╠═3d1fca00-3449-11eb-3e0e-5784380ebdbc
# ╠═4ab2f908-3449-11eb-299d-53d5a0a3044e
# ╠═473a9452-3449-11eb-1dd1-25d770220737
# ╠═6e67d06a-344b-11eb-15c6-ed3cfcfd8a83


text:
"takes a json biomd metadata file and converts to DF"
function jsonfn_to_df(fn)
    json = read(fn, String);
    json = JSON3.read(json)
    haskey(json, :models) ? DataFrame(jsontable(json.models)) : missing
end

"curls the json files of all the "
function curl_biomd_metadata(meta_dir="$(datadir)/biomd_ode_meta")
    !ispath(meta_dir) && mkpath(meta_dir)
    offsets = 0:100:2200
    # "https://www.ebi.ac.uk/biomodels/search?query=&offset=0&numResults=10"
    sbml_query = "sbml"
    sbml_ode_query = "*%3A*%20AND%20modellingapproach%3A%22Ordinary%20differential%20equation%20model%22&domain=biomodels"
    urls = .*("https://www.ebi.ac.uk/biomodels/search?query=", sbml_ode_query,"&offset=", string.(offsets), "&numResults=100&format=json")
    # @sync for i in 1:length(urls)
    @sync for i in 1:10
        # Downloads.download(urls[i], "$(meta_dir)/sbml_$(i).json")
        @async run(`curl $(urls[i]) -o "$(meta_dir)/sbml_$(i).json"`)
    end
end

function biomd_metadata(meta_dir="$(datadir)/biomd_meta"; curl_meta=false)
    curl_meta && curl_biomd_metadata(meta_dir)
    fns = readdir(meta_dir; join=true)
    dfs = filter(!ismissing, jsonfn_to_df.(fns))
    vcat(dfs...)
end

"this downloads the xmls directly, without needing zip "
function curl_biomd_xmls(ids; verbose=false)
    base = "https://www.ebi.ac.uk/biomodels/model/download/"
    @sync Threads.@threads for id in ids
        verbose && @info("downloading $id")
        url = "$(base)$(id)?filename=$(id)_url.xml"
        fn = "$(datadir)/biomd/$(id).xml"
        # Downloads.download(url, fn)
        run(`curl $(url) -o "$(fn)"`)
    end
end

function biomodels(
    meta_dir="$(datadir)/biomd_meta",
    # zips_dir="$(datadir)/biomd_zips/",
    biomd_dir="$(datadir)/biomd/"; # the xml files are put here
    curl_meta=false,
    limit=nothing, # still gets all the metadata, just limits for curling zips
    verbose=true
    )
    verbose && @info("in biomodels()")

    mkpath.([meta_dir,
        # zips_dir,
        biomd_dir])

    df = biomd_metadata(meta_dir; curl_meta=curl_meta)
    verbose && display(df)
    CSV.write("$(datadir)/sbml_biomodels.csv", df)
    limit === nothing ? curl_biomd_xmls(df.id; verbose=verbose) : curl_biomd_xmls(df.id[1:limit]; verbose=verbose)
    df
end

###################################################
# I think below is all garbage
"uses the BioModels IDs and the download REST API"
function biomd_zip_urls(ids)
    base = "https://www.ebi.ac.uk/biomodels/search/download?models="
    N = 100 # api limits 100 at a time
    chunks = Iterators.partition(ids, N)
    qs = join.(chunks, ",")
    base .* qs
end

"""
takes the metadata dataframe from `biomd_metadata()`.

should probably do this async
"""
function curl_biomd_zips(urls, zips_dir="$(datadir)/biomd_zips/")
    for i in 1:length(urls) # @threads seems to not work
        run(`curl -X GET "$(urls[i])" -H "accept: application/zip" -o $(zips_dir)$i.zip`)
    end
    return urls
end

"needs unzip in shell path"
function unzip_biomd(zips_dir, unzip_dir)
    mkpath(unzip_dir)
    zips = readdir(zips_dir; join=true)
    for fn in zips
        run(`unzip $fn -d $(unzip_dir)`)
    end
    unzip_dir
end

# "to delete i guess"
# function grab_extract_biomd_zips()
#     urls = limit === nothing ? biomd_zip_urls(df.id) : biomd_zip_urls(df.id[1:limit])
#     curl_biomd_zips(urls, zips_dir)
#     unzip_biomd(zips_dir, unzip_dir)
# end

function ab()
    ids = df.id
    ms = []
    @sync for id in ids
        @async try
            m = readSBML("data/biomd/$id.xml")
            push!(ms, id => m)
        catch e
            @info e
        end
    end
    @show length(ms)
    df[first.(ms) .== df.id, :]
    filter(x -> x.id ∈ good, df)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using BenchmarkTools

f(n) = sum(abs2, 1:n)
g(n) = sum(abs2(x) for x in 1:n)

@show VERSION
@btime f(10^6)
@btime g(10^6)

# %%
using BenchmarkTools

f(n) = sum(abs2, 1:n)
g(n) = sum(abs2(x) for x in 1:n)

@show VERSION
@btime f(10^6)
@btime g(10^6)

# %%
using BenchmarkTools

f(n) = sum(abs2, 1:n)
g(n) = sum(abs2(x) for x in 1:n)

@show VERSION
@btime f(10^6)
@btime g(10^6)

# %%

max_stars_repo_path: 0005/regression of generator.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/regression of generator.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/regression of generator.jl

max_forks_repo_name: genkuroki/public


text:
# Autogenerated wrapper script for isoband_jll for armv6l-linux-musleabihf
export libisoband

JLLWrappers.@generate_wrapper_header("isoband")
JLLWrappers.@declare_library_product(libisoband, "libisoband.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libisoband,
        "lib/libisoband.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# io.jl

module IOTest

using Test
using QDates

@testset "string/show representation of QDate" begin
    qdt = QDates.QDate(2018, 11, 27)
    @test string(qdt) == sprint(show, qdt) == "旧2018年11月27日"
    qdtl = QDates.QDate(2017, 5, true, 1)
    @test string(qdtl) == sprint(show, qdtl) == "旧2017年閏05月01日"
end

end


text:
using Plots
pyplot()
using LaTeXStrings
using StatPlots     # needed for box plot support

n = 500
x = randn(n)  # N(0, 1)
x = exp.(x)   # Map x to lognormal
y = randn(n) + 2.0  # N(2, 1)
z = randn(n) + 4.0  # N(4, 1)
data = vcat(x, y, z)
l = [LaTeXString("\$X\$") LaTeXString("\$Y\$")  LaTeXString("\$Z\$") ]
xlabels = reshape(repmat(l, n), n*3, 1)

boxplot(xlabels, data, label="", ylims=(-2, 14))


text:
using GridInterpolations
using LocalFunctionApproximation
using StaticArrays
using POMDPs
using JLD2, FileIO
using Random
using ContinuumWorld
using CMSSPs

trials = 5
rng = MersenneTwister(2973)

params_fn = "grids-continuum-params1.toml"
params = continuum_parse_params(params_fn)

inhor_fn = "grids-continuum-params1-inhor.jld2"
outhor_fn = "grids-continuum-params1-outhor.jld2"
modal_horizon_policy = load_modal_horizon_policy_localapprox(inhor_fn, outhor_fn)

# Create CWorld actions
grid_cworld_actions = [Vec2(params.move_amt, 0.0), Vec2(-params.move_amt, 0.0),
                       Vec2(0.0, params.move_amt), Vec2(0.0, -params.move_amt), Vec2(0.0, 0.0)]

# Use one cworld for all of them
cworld = CWorld(xlim = (0.0, 1.0), ylim = (0.0, 1.0),
                reward_regions = [], rewards = [], terminal = [],
                stdev = params.move_std, actions = grid_cworld_actions, discount = 1.0)
cworlds = Dict(1=>cworld, 2=>cworld, 3=>cworld, 4=>cworld)
params.cworlds = cworlds

# Only need to find policy for one mode, since shared
modal_mdp = GridsContinuumMDPType(1, params, grid_cworld_actions, params.beta, params.horizon_limit)

for i = 1:trials

    start_horizon = convert(Int64, rand(rng, params.horizon_limit/2:params.horizon_limit))
    target_state = sample_continuum_state(rng)
    curr_state = Vec2(0.0, 0.0)
    tp_dist = TPDistribution([start_horizon],[1.0])

    @show target_state, start_horizon

    for t = 0:1:start_horizon-1
        a = get_best_intramodal_action(modal_horizon_policy, t, tp_dist, curr_state, target_state)
        if a == nothing
            @warn "ABORT!"
            break
        end

        (new_state, reward) = generate_sr(modal_mdp, curr_state, a.action, rng)
        curr_state = new_state
    end

    if isterminal(modal_mdp, curr_state - target_state)
        @info "Success!"
    else
        @info "Failure!"
    end
end


text:
function x()
	a = 12320657
	b = 9659666
    m = 20201227
    e = 1
    l = 0
    while e != a
        e *= 7
        e %= m
        l += 1
    end
    k = 1
    for i=1:l
        k *= b
        k %= m
    end
    println(k)
end


text:
"""
`Basics.estimate()`
    ... estimates various quantities, typically, due to semi-empirical formulas and parametrizations.

  + `("ionization potential: inner-shell", shell::Shell, Z::Int64)`
    ... to estimate the ionization potential = mean binding energy of an electron in the given shell; an energy value::Float64
        is returned. These ionization potentials are taken from the tabulation ....
"""
function Basics.estimate(sa::String, shell::Shell, Z::Int64)
    if      sa == "ionization potential: inner-shell"
        shellIndex = Dict( Subshell("1s_1/2") => 1, Subshell("2s_1/2") => 2, Subshell("2p_1/2") => 3, Subshell("2p_3/2") => 4,
                           Subshell("3s_1/2") => 5, Subshell("3p_1/2") => 6, Subshell("3p_3/2") => 7,
                           Subshell("3d_3/2") => 8, Subshell("3d_3/2") => 9  )
        idx = shellIndex[shell]
        return( JAC.store("binding energies: Williams (2000)", Z)[idx] )
    else
        error("Unsupported keystring = $sa")
    end
end


text:
  # Meus testes
  function teste1(n)
    a = 19
    println(a*f(n))
    println(f2(3*n,n))
  end

  teste1(10)


text:
#InspectDR: Base functionnality and types for Cairo layer
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#==Main types
===============================================================================#
#2D plot canvas; basic info:
mutable struct PCanvas2D <: PlotCanvas
	ctx::CairoContext
	bb::BoundingBox #Entire canvas
	graphbb::BoundingBox #Graph portion
	ext::PExtents2D #Extents of graph portion
	xf::Transform2D #Transform used to render data
end
PCanvas2D(ctx, bb, graphbb, ext) =
	PCanvas2D(ctx, bb, graphbb, ext, Transform2D(ext, graphbb))
PCanvas2D(ctx, bb, ginfo::Graph2DInfo) =
	PCanvas2D(ctx, bb, ginfo.graphbb, ginfo.ext, ginfo.xf)

#=
#TODO: use instead of PCanvas2D???
#TODO: Add axis scale info, etc???
#Canvas for 2D Graph:
mutable struct GCanvas2D
	ctx::CairoContext
	graphbb::BoundingBox
	ext::PExtents2D #Extents of graph portion
	xf::Transform2D #Transform used to render data
end
GCanvas2D(ctx, graphbb, ext) =
	GCanvas2D(ctx, graphbb, ext, Transform2D(ext, graphbb))
=#

#Used to buffer portions of plot for better GUI response times
mutable struct CairoBufferedPlot
	surf::Cairo.CairoSurface #Main surface where plot is drawn
	data::Cairo.CairoSurface #Cache of data image layer
end

#==Basic rendering
===============================================================================#

#Apply transform to Cairo (when easier to do so):
function setxfrm(ctx::CairoContext, xf::Transform2D)
	xorig = xf.x0*xf.xs
	yorig = xf.y0*xf.ys
	Cairo.translate(ctx, xorig, yorig)
	Cairo.scale(ctx, xf.xs, xf.ys)
end

#==Rendering Glyphs
===============================================================================#

function getglyphcolor(glyph::GlyphAttributes, line::LineAttributes)
	c = glyph.color
	if nothing == c
		c = line.color
	end
	return c
end

function getglyphfill(glyph::GlyphAttributes)
	fill = glyph.fillcolor
	if COLOR_TRANSPARENT == fill #Small optimization
		fill = nothing
	end
	return fill
end

function drawglyph(ctx::CairoContext, g::GlyphPolyline, pt::Point2D, size::DReal, fill)
	#TODO: could scale GlyphPolyline only once for efficiency (instead of every time we draw).
	x = g.x*size; y = g.y*(-1*size) #Compensate for y-reversal
	Cairo.move_to(ctx, pt.x, pt.y)
	Cairo.rel_move_to(ctx, x[1], y[1]) #rel: Probably less prone to rounding errors

	dxv = diff(x); dyv = diff(y)
	for (dx, dy) in zip(dxv, dyv)
		Cairo.rel_line_to(ctx, dx, dy)
	end
	if g.closepath
		Cairo.close_path(ctx)
		renderfill(ctx, fill)
	end
	Cairo.stroke(ctx)
	return
end
function drawglyph(ctx::CairoContext, g::GlyphLineSegments, pt::Point2D, size::DReal, fill)
	dxv = size * (g.x2 - g.x1)
	dyv = -size * (g.y2 - g.y1)
	x1v = size * g.x1
	y1v = -size * g.y1
	for (x1, y1, dx, dy) in zip(x1v, y1v, dxv, dyv)
		Cairo.move_to(ctx, pt.x, pt.y)
		Cairo.rel_move_to(ctx, x1, y1) #rel: Probably less prone to rounding errors
		Cairo.rel_line_to(ctx, dx, dy)
		Cairo.stroke(ctx)
	end
	return
end
function drawglyph(ctx::CairoContext, g::GlyphCircle, pt::Point2D, size::DReal, fill)
	cairo_circle(ctx, pt.x, pt.y, g.radius*size)
	renderfill(ctx, fill)
	Cairo.stroke(ctx)
end

#Correctly displays a glyph, given wfrm properties.
function drawglyph_safe(ctx::CairoContext, wfrm::DWaveform, pt::Point2D)
	_glyph = Glyph(wfrm.glyph.shape)
	if nothing == _glyph; return; end

	linecolor = getglyphcolor(wfrm.glyph, wfrm.line)
	linestyle = LineStyle(:solid, Float64(wfrm.line.width), linecolor)
	setlinestyle(ctx, linestyle)
	fill = getglyphfill(wfrm.glyph)
	gsize = Float64(wfrm.glyph.size)

	drawglyph(ctx, _glyph, pt, gsize, fill)
	return
end

#==Rendering base plot elements
===============================================================================#

#Render main plot annotation (titles, axis labels, ...)
#-------------------------------------------------------------------------------
function render(ctx::CairoContext, a::Annotation,
	bb::BoundingBox, databb::BoundingBox, graphbblist::Vector{BoundingBox}, lyt::PlotLayout)
	TIMESTAMP_OFFSET = 3 #WANTCONST

	#Title
#	xcenter = (bb.xmin+bb.xmax)/2 #Entire plot BB.
	xcenter = (databb.xmin+databb.xmax)/2 #Data-area BB only.
	pt = Point2D(xcenter, bb.ymin+lyt.voffset_title)
	render(ctx, a.title, pt, lyt.font_title, align=ALIGN_HCENTER|ALIGN_VCENTER)

	#X-axis label
	xcenter = (databb.xmin+databb.xmax)/2
	pt = Point2D(xcenter, bb.ymax-lyt.voffset_xaxislabel)
	render(ctx, a.xlabel, pt, lyt.font_axislabel, align=ALIGN_HCENTER|ALIGN_VCENTER)

	#Y-axis labels
	nstrips = min(length(a.ylabels), length(graphbblist))
	for i in 1:nstrips
		graphbb = graphbblist[i]
		ycenter = (graphbb.ymin+graphbb.ymax)/2
		pt = Point2D(bb.xmin+lyt.hoffset_yaxislabel, ycenter)
		render(ctx, a.ylabels[i], pt, lyt.font_axislabel, align=ALIGN_HCENTER|ALIGN_VCENTER, angle=-π/2)
	end

	#Time stamp
	if lyt.enable_timestamp
		pt = Point2D(bb.xmax-TIMESTAMP_OFFSET, bb.ymax-TIMESTAMP_OFFSET)
		render(ctx, a.timestamp, pt, lyt.font_time, align=ALIGN_RIGHT|ALIGN_BOTTOM)
	end
end

#Render frame around graph
#-------------------------------------------------------------------------------
function render_graphframe(canvas::PCanvas2D, aa::AreaAttributes)
	ctx = canvas.ctx #WANTCONST
Cairo.save(ctx)
	setlinestyle(ctx, LineStyle(aa.line))
	Cairo.rectangle(ctx, canvas.graphbb)
	Cairo.stroke(ctx)
Cairo.restore(ctx)
	return
end

#==Render grid
===============================================================================#
render_vlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, xlines::AbstractGridLines) = nothing
function render_vlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, xlines::GridLines)
	if xlines.displaymajor
		setlinestyle(ctx, lyt.line_gridmajor)
		for xline in xlines.major
			x = map2dev(xf, Point2D(xline, 0)).x
			drawline(ctx, Point2D(x, graphbb.ymin), Point2D(x, graphbb.ymax))
		end
	end
	if xlines.displayminor
		setlinestyle(ctx, lyt.line_gridminor)
		for xline in xlines.minor
			x = map2dev(xf, Point2D(xline, 0)).x
			drawline(ctx, Point2D(x, graphbb.ymin), Point2D(x, graphbb.ymax))
		end
	end
end
render_hlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, ylines::AbstractGridLines) = nothing
function render_hlines(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, ylines::GridLines)
	if ylines.displaymajor
		setlinestyle(ctx, lyt.line_gridmajor)
		for yline in ylines.major
			y = map2dev(xf, Point2D(0, yline)).y
			drawline(ctx, Point2D(graphbb.xmin, y), Point2D(graphbb.xmax, y))
		end
	end
	if ylines.displayminor
		setlinestyle(ctx, lyt.line_gridminor)
		for yline in ylines.minor
			y = map2dev(xf, Point2D(0, yline)).y
			drawline(ctx, Point2D(graphbb.xmin, y), Point2D(graphbb.xmax, y))
		end
	end
end

function render_grid(canvas::PCanvas2D, lyt::PlotLayout, grid::GridRect)
	ctx = canvas.ctx #WANTCONST
	Cairo.save(ctx) #-----
	render_vlines(ctx, canvas.graphbb, canvas.xf, lyt, grid.xlines)
	render_hlines(ctx, canvas.graphbb, canvas.xf, lyt, grid.ylines)
	Cairo.restore(ctx) #-----
end

#==Render ticks
===============================================================================#

#Default label:
function render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, ::AxisScale)
	tstr = ""
	if false #Old display system.
		#TODO: deprecate
		tstr = "$val"
		if length(tstr) > 7 #HACK!
			tstr = @sprintf("%0.1e", val)
		end
	else
		tstr = formatted(val, fmt.fmt, showexp=!fmt.splitexp)
	end
	render(ctx, tstr, pt, font, align=align)
end

render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, scale::LogScale{:e}) =
	render_power(ctx, "e", val, pt, font, align)
render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, scale::LogScale{2}) =
	render_power(ctx, "2", val, pt, font, align)
render_ticklabel(ctx::CairoContext, val::DReal, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, scale::LogScale{10}) =
	render_power(ctx, "10", val, pt, font, align)

function render_axisscalelabel(ctx::CairoContext, pt::Point2D, font::Font, align::CAlignment, fmt::TickLabelFormatting, ::AxisScale)
	tstr = formatted_exp(fmt.fmt)
	render(ctx, tstr, pt, font, align=align)
end

#Render ticks: Well-defined GridLines
#-------------------------------------------------------------------------------
function render_xticks(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, xlines::GridLines, xs::AxisScale, ticklabels::Bool)
	tframe = DReal(lyt.frame_data.line.width) #WANTCONST TODO: Fix LineAttributes to have concrete type
	fmt = TickLabelFormatting(lyt.format_xtick, xlines.rnginfo)
	yframe = graphbb.ymax
	ylabel = graphbb.ymax + lyt.voffset_xticklabel
	for xtick in xlines.major
		x = map2dev(xf, Point2D(xtick, 0)).x
		if ticklabels
			render_ticklabel(ctx, xtick, Point2D(x, ylabel), lyt.font_ticklabel, ALIGN_TOP|ALIGN_HCENTER, fmt, xs)
		end
		drawline(ctx, Point2D(x, yframe), Point2D(x, yframe-lyt.length_tickmajor))
	end
	for xtick in xlines.minor
		x = map2dev(xf, Point2D(xtick, 0)).x
		drawline(ctx, Point2D(x, yframe), Point2D(x, yframe-lyt.length_tickminor))
	end
	if fmt.splitexp && ticklabels
		xlabel = graphbb.xmax + tframe
		render_axisscalelabel(ctx, Point2D(xlabel, yframe), lyt.font_ticklabel, ALIGN_BOTTOM|ALIGN_LEFT, fmt, xs)
	end
end
function render_yticks(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, ylines::GridLines, ys::AxisScale)
	tframe = DReal(lyt.frame_data.line.width) #WANTCONST TODO: Fix LineAttributes to have concrete type
	fmt = TickLabelFormatting(lyt.format_ytick, ylines.rnginfo)
	xframe = graphbb.xmin
	xlabel = graphbb.xmin - lyt.hoffset_yticklabel
	for ytick in ylines.major
		y = map2dev(xf, Point2D(0, ytick)).y
		render_ticklabel(ctx, ytick, Point2D(xlabel, y), lyt.font_ticklabel, ALIGN_RIGHT|ALIGN_VCENTER, fmt, ys)
		drawline(ctx, Point2D(xframe, y), Point2D(xframe+lyt.length_tickmajor, y))
	end
	for ytick in ylines.minor
		y = map2dev(xf, Point2D(0, ytick)).y
		drawline(ctx, Point2D(xframe, y), Point2D(xframe+lyt.length_tickminor, y))
	end
	if fmt.splitexp
		ylabel = graphbb.ymin - tframe
		render_axisscalelabel(ctx, Point2D(xframe, ylabel), lyt.font_ticklabel, ALIGN_BOTTOM|ALIGN_LEFT, fmt, ys)
	end
end

#Render ticks: UndefinedGridLines
#-------------------------------------------------------------------------------
function render_xticks(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, xlines::UndefinedGridLines, xs::AxisScale, ticklabels::Bool)
	fmt = TickLabelFormatting(NoRangeDisplayInfo())
	yframe = graphbb.ymax
	ylabel = graphbb.ymax + lyt.voffset_xticklabel
	for (x, xlabel) in [(graphbb.xmin, xlines.minline), (graphbb.xmax, xlines.maxline)]
		if ticklabels
			render_ticklabel(ctx, xlabel, Point2D(x, ylabel), lyt.font_ticklabel, ALIGN_TOP|ALIGN_HCENTER, fmt, LinScale())
		end
		drawline(ctx, Point2D(x, yframe), Point2D(x, yframe-lyt.length_tickmajor))
	end
end
function render_yticks(ctx::CairoContext, graphbb::BoundingBox, xf::Transform2D, lyt::PlotLayout, ylines::UndefinedGridLines, ys::AxisScale)
	fmt = TickLabelFormatting(NoRangeDisplayInfo())
	xframe = graphbb.xmin
	xlabel = graphbb.xmin - lyt.hoffset_yticklabel
	for (y, ylabel) in [(graphbb.ymax, ylines.minline), (graphbb.ymin, ylines.maxline)]
		render_ticklabel(ctx, ylabel, Point2D(xlabel, y), lyt.font_ticklabel, ALIGN_RIGHT|ALIGN_VCENTER, fmt, LinScale())
		drawline(ctx, Point2D(xframe, y), Point2D(xframe+lyt.length_tickmajor, y))
	end
end

function render_ticks(canvas::PCanvas2D, lyt::PlotLayout, grid::GridRect, xs::AxisScale, ys::AxisScale, xticklabels::Bool)
	render_xticks(canvas.ctx, canvas.graphbb, canvas.xf, lyt, grid.xlines, xs, xticklabels)
	render_yticks(canvas.ctx, canvas.graphbb, canvas.xf, lyt, grid.ylines, ys)
end

#==High-level rendering
===============================================================================#

#Render axis labels, ticks, ...
#-------------------------------------------------------------------------------
function render_axes(canvas::PCanvas2D, lyt::PlotLayout, grid::GridRect, xs::AxisScale, ys::AxisScale, xticklabels::Bool)
	render_graphframe(canvas, lyt.frame_data)
	render_ticks(canvas, lyt, grid, xs, ys, xticklabels)
end

#Render NaNs on plot:
#-------------------------------------------------------------------------------
function rendernans(canvas::PCanvas2D, wfrm::IWaveform)
	NAN_LINE = LineStyle(
		:solid, Float64(4), ARGB32(1, 0, 0, 0.5)
	) #WANTCONST
	ctx = canvas.ctx #WANTCONST
	graphbb = canvas.graphbb #WANTCONST
	x = wfrm.ds.x #WANTCONST
	y = wfrm.ds.y #WANTCONST

	if length(x) != length(y)
		error("x & y - vector length mismatch.")
	end

	setlinestyle(ctx, NAN_LINE)
	hasNaNNaN = false

	for i in 1:length(x)
		pt = Point2D(x[i], y[i])
		xnan = isnan(pt.x); ynan = isnan(pt.y)

		if xnan && ynan
			hasNaNNaN = true #nothing practical to display
		elseif xnan
			pt = map2dev(canvas.xf, pt)
			Cairo.move_to(ctx, graphbb.xmin, pt.y)
			Cairo.line_to(ctx, graphbb.xmax, pt.y)
			Cairo.stroke(ctx)
		elseif ynan
			pt = map2dev(canvas.xf, pt)
			Cairo.move_to(ctx, pt.x, graphbb.ymin)
			Cairo.line_to(ctx, pt.x, graphbb.ymax)
			Cairo.stroke(ctx)
		end
	end

	#TODO: warn hasNaNNaN?
	return hasNaNNaN
end
function rendernans(canvas::PCanvas2D, wfrmlist::Vector{IWaveform}, strip::Int)
	for wfrm in wfrmlist
		if 0 == wfrm.strip || wfrm.strip == strip
			rendernans(canvas, wfrm)
		end
	end
	return
end

#Render an actual waveform
#-------------------------------------------------------------------------------
function render(canvas::PCanvas2D, wfrm::DWaveform)
	ctx = canvas.ctx #WANTCONST
	ds = wfrm.ds #WANTCONST

	if length(ds) < 1; return; end

if hasline(wfrm)
	setlinestyle(ctx, LineStyle(wfrm.line))
	newsegment = true
	for i in 1:length(ds)
		pt = map2dev(canvas.xf, ds[i])
		xnan = isnan(pt.x); ynan = isnan(pt.y)

		#TODO: log NaNs?  Display NaNs?
		if xnan && ynan
			Cairo.stroke(ctx) #Close last line
			newsegment = true
		elseif xnan || ynan
			#Simply skip this point
		elseif newsegment
			Cairo.move_to(ctx, pt.x, pt.y)
			newsegment = false
		else
			Cairo.line_to(ctx, pt.x, pt.y)
		end
	end
	Cairo.stroke(ctx)
end

	_glyph = Glyph(wfrm.glyph.shape)
	if nothing == _glyph; return; end

	#TODO: do not draw when outside graph extents.

	#Draw symbols:
	linecolor = getglyphcolor(wfrm.glyph, wfrm.line)
	linestyle = LineStyle(:solid, Float64(wfrm.line.width), linecolor)
	setlinestyle(ctx, linestyle)
	fill = getglyphfill(wfrm.glyph)
	gsize = Float64(wfrm.glyph.size)
	for i in 1:length(ds)
		pt = map2dev(canvas.xf, ds[i])
		drawglyph(ctx, _glyph, pt, gsize, fill)
	end

	return
end
function render(canvas::PCanvas2D, wfrmlist::Vector{DWaveform}, strip::Int)
	for wfrm in wfrmlist
		if 0 == wfrm.strip || wfrm.strip == strip
			render(canvas, wfrm)
		end
	end
	return
end

#Last line


text:
#=
    This code was generated by heta-compiler *

=#

__platform__ = (function()

### create default constants
mm_constants_num_ = NamedTuple{(
  :Vmax,:Km,
)}(Float64[
  0.1,2.5,
])

### create default observables
mm_records_output_ = NamedTuple{(
  :default_comp,:S,:P,:r1,
)}(Bool[
  false,true,true,false,
])

### create default events
mm_events_active_ = NamedTuple{(

)}(Bool[

])

### initialization of ODE variables and Records
function mm_init_func_(cons)
    #(Vmax,Km,) = cons

    # Heta initialize
    t = 0.0 # initial time
    P = 0.0
    S = 10.0
    default_comp = 1.0
    r1 = cons[1] * S / (cons[2] + S) * default_comp

    # save results

    return (
        Float64[
            S * default_comp,
            P * default_comp,
        ],
        Float64[
            default_comp,
        ]
    )
end

### calculate RHS of ODE
function mm_ode_func_(du, u, p, t)
    cons = p.constants
    (default_comp,) = p.static
    (S_,P_,) = u

    # Heta rules
    P = P_ / default_comp
    S = S_ / default_comp
    r1 = cons[1] * S / (cons[2] + S) * default_comp

    #p.static .= [default_comp,]
    du .= [
      -r1,  # dS_/dt
      r1,  # dP_/dt
    ]
end

### output function
function mm_saving_generator_(outputIds::Vector{Symbol})
    function saving_(u, t, integrator)
        cons = integrator.p.constants
        (default_comp,) = integrator.p.static
        (S_,P_,) = u

        # Heta rules
        P = P_ / default_comp
        S = S_ / default_comp
        r1 = cons[1] * S / (cons[2] + S) * default_comp

        # force amount

        d = Base.@locals
        return [d[id] for id in outputIds]
    end
end

### TIME EVENTS ###

### D EVENTS ###

### STOP EVENTS ###

### event assignments

### MODELS ###

mm_model_ = (
  mm_init_func_,
  mm_ode_func_,
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  mm_saving_generator_,
  mm_constants_num_,
  mm_events_active_,
  mm_records_output_
)

### OUTPUT ###

return (
  (
    mm = mm_model_,
  ),
  (),
  "*"
)

end
)()

max_stars_repo_name: hetalang/heta-compiler

max_issues_repo_name: hetalang/heta-compiler

max_forks_repo_name: hetalang/heta-compiler


text:
using OffsetArrays: Origin
@testset "Origin" begin
    get_origin(A::AbstractArray) = first.(axes(A))

    @test Origin(0) != Origin((0, ))
    @test Origin(CartesianIndex(1, 2)) === Origin((1, 2)) === Origin(1, 2)

    # 1d
    v = [1, 2]
    @test get_origin(OffsetArray(v, Origin(2))) == (2, )
    ov = OffsetArray(v, -3)
    @test get_origin(OffsetArray(ov, Origin(2))) == (2, )
    @test get_origin(OffsetVector(ov, Origin(2))) == (2, )
    @test get_origin(OffsetArray(ov, Origin((2, )))) == (2, )

    # 2d
    a = [1 2;3 4]
    @test get_origin(OffsetArray(a, Origin(0))) == (0, 0)
    oa = OffsetArray(a, -3, -3)
    @test get_origin(OffsetArray(oa, Origin(0))) == (0, 0)
    @test get_origin(OffsetMatrix(oa, Origin(0))) == (0, 0)
    @test get_origin(OffsetArray(oa, Origin(1, 2))) == (1, 2)

    # 3d
    a = ones(3, 3, 3)
    @test get_origin(OffsetArray(a, Origin(0))) == (0, 0, 0)
    oa = OffsetArray(a, -3, -3, -3)
    @test get_origin(OffsetArray(oa, Origin(0))) == (0, 0, 0)
    @test get_origin(OffsetArray(oa, Origin(1, 2, 3))) == (1, 2, 3)

    # Scalar broadcasting
    let
        a = [ [1,2,3], [4,5,6] ]
        oa = OffsetVector.(a, Origin(0))
        @test get_origin.(oa) == [ (0,), (0,) ]

        a = [ [1 2; 3 4], [5 6 7; 8 9 10] ]
        oa = OffsetArray.(a, Origin(0, -1))
        @test get_origin.(oa) == [ (0,-1), (0,-1) ]
    end
end


text:
function parseinput(lines)
    sep1, sep2 = findall(==(""), lines)
    rules = Dict{String,Vector{UnitRange{Int}}}()
    for i in 1:sep1 - 1
        key, values = split(lines[i], ": ")
        regex = r"(\d+)-(\d+)"
        rule = [parse(Int, m.captures[1]):parse(Int, m.captures[2]) for m in eachmatch(regex, values)]
        rules[key] = rule
    end
    parseticket = x -> parse.(Int, [m.match for m in eachmatch(r"\d+", x)])
    ticket = parseticket(lines[sep1 + 2])
    tickets = Vector{Int}[]
    for i in sep2 + 2:length(lines)
        push!(tickets, parseticket(lines[i]))
    end
    return rules, ticket, tickets
end

function invalid(ticket::Vector{Int}, rules_values)
    mask = [sum(sum(t .∈ r) for r in rules_values) == 0 for t in ticket]
    return  sum(ticket[mask])
end

const rules, ticket, tickets = parseinput(readlines("day16.txt"))

println("Part 1: $(sum(invalid(t, values(rules)) for t in tickets))")

function isvalidticket(ticket::Vector{Int}, rules_values)
    return all(sum(sum(t .∈ r) for r in rules_values) >= 1 for t in ticket)
end

function part2(tickets::Vector{Vector{Int}}, rules::Dict{String,Vector{UnitRange{Int}}})
    valid_tickets = filter(x -> isvalidticket(x, values(rules)), tickets)
    # create a 2D array where each row is a ticket and each colum a field
    T = [t[i] for t in valid_tickets, i = 1:length(valid_tickets[1])]
    # We want to check each column of the array, to see how many rules can apply. We
    # will find one rule that has only one possibility, we will remove it from the
    # others possible values and so on, until we know which field corresponds to which
    # rule
    possibilities = Dict{String,Vector{Int}}()
    num_tickets = length(valid_tickets)
    for (name, rule) in rules
        # this is ugly because it collects everything, but I had trouble with
        # broadcasting when I kept the 2 unit ranges.
        ranges = union(rule...)
        # Ref is necessary to avoid the broadcasting, we just want to check is each
        # element of the column (i.e. each field) is possible with this rule
        # We sum to know if the rule is valid for all tickets and findall to know which
        # columns are applicable for this rule
        possibles = findall([sum(col .∈ Ref(ranges)) == num_tickets for col in eachcol(T)])
        possibilities[name] = possibles
    end
    seen = Set{Int}()
    fields = Dict{String,Int}()
    sorted = sort(collect(possibilities), by=x -> length(x[2]))
    for (name, values) in sorted
        field = setdiff(values, seen)[1]
        seen = seen ∪ values
        fields[name] = field
    end
    return prod(ticket[fields[i]] for i in keys(fields) if startswith(i, "departure"))
end

println("Part 2: $(part2(tickets, rules))")


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

"""
    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)

The classical inviscid Taylor-Green vortex.
"""
function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
  A  = 1.0 # magnitude of speed
  Ms = 0.1 # maximum Mach number

  rho = 1.0
  v1  =  A * sin(x[1]) * cos(x[2]) * cos(x[3])
  v2  = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
  v3  = 0.0
  p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
  p   = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3]))

  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
end
initial_condition = initial_condition_taylor_green_vortex

volume_flux = flux_ranocha
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-1.0, -1.0, -1.0) .* pi
coordinates_max = ( 1.0,  1.0,  1.0) .* pi
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                n_cells_max=100_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.4)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using Jtac

chain = @stack Morris stack_input=true begin
  Conv(128, relu)
  Dense(200, relu)
end

model = NeuralModel(Morris, chain)
player = MCTSPlayer(model, power = 100)

loss = Loss(value = 10, regularization = 1e-5)

opponents = [MCTSPlayer(power = p) for p in [10, 25, 50, 100, 250, 500]]

with_contest( train_self!
            , player
            , loss = loss
            , playings = 500
            , batchsize = 100
            , branching = 0.1
            , epochs = 10
            , opponents = opponents
            , interval = 3
            , length = 1000
            , cache = 2000 )

max_stars_repo_name: roSievers/jtac


text:
# LSS solve
using SparseArrays
using BlockArrays
function lsssolve(R, b)
    d_u, n = size(b)
    ndu = n*d_u
    eye = sparse(Matrix([zeros(ndu, d_u) 1.0I(ndu)]))
    D = (BlockArray{Float64}(zeros(ndu,ndu),
               d_u*ones(Int64,n),
               d_u*ones(Int64,n)))
    [D[Block(i,i)] = R[:,:,i+1] for i =1:n-1]
    D = sparse([Array(D) zeros(ndu, d_u)])
    B = D - eye
    BB = B*transpose(B)
	a = -transpose(B)*(BB\b[:])
    a = reshape(a, d_u, n+1)[:,1:end-1]
	return a
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

run tests using HSL sparse linear system solver caches
requires that HSL.jl be installed and built successfully; Requires.jl handles this optional dependency
=#

using Test
import HSL
import Hypatia
import Hypatia.Solvers
blas_reals = [
    Float64,
    Float32,
    ]
include(joinpath(@__DIR__, "nativeinstances.jl"))
include(joinpath(@__DIR__, "nativesets.jl"))

options = (verbose = false,)

@testset "HSL cache tests" begin

@testset "cache setup: $cache_type" for cache_type in [Hypatia.HSLSymCache] # TODO wrap and test a HSLNonSymCache
    cache = cache_type()
    @test !cache.analyzed
    @test cache_type{Float32}().analyzed == false
    @test_throws Exception cache_type{BigFloat}()
end

@testset "SymIndefSparse tests: $inst_name, $T" for inst_name in inst_cones_many,
    T in blas_reals
    inst_function = eval(Symbol(inst_name))
    inst_function(T, solver = Solvers.Solver{T}(syssolver =
        Solvers.SymIndefSparseSystemSolver{T}(
        fact_cache = Hypatia.HSLSymCache{T}()); options...))
end

end
;


text:
using CSV

# TODO
#function series_eq(x, y) end
#function read_csv(fh, **kwargs) end
#function sub_chr(s, chrom) end
#function get_present_chrs(fh, num) end
#function which_compression(fh) end
#function read_cts(fh, match_snps) end
#function ldscore_fromlist(flist, num=nothing) end
#function l2_parser(fh, compression) end
#function annot_parser(
#    fh, compression, frqfile_full=nothing, compression_frq=nothing,
#) end
#function frq_parser(fh, compression) end
#function ldscore(fh, num=nothing) end
#function M(fh; num=nothing, N=2, common=False) end
#function M_fromlist(flist; num=nothing, N=2, common=False) end
#function annot(fh_list; num=nothing, frqfile=nothing) end
#function __ID_List_Factory__(
#    colnames, keepcol, fname_end; header=nothing, usecols=nothing,
#) end
#    struct IDContainer end
#        function __init__(container::IDContainer, fname) end
#        function __read__(container::IDContainer, fname) end
#        function loj(container::IDContainer, externalDf) end

function get_compression(fh)
    if endswith(fh, "gz") return "gzip"
    elseif endswith(fh, "bz2") return "bz2"
    else return nothing end
end

function parse_sumstats(fh; alleles=false, dropna=true)
    @info "Reading from:" fh
    x = CSV.File(fh)

    if dropna x = x.dropna(how="any") end # TODO double check Julian equivalent

    @info x
    return x
end