text:
using KeplerTools
using LinearAlgebra
using Test
using Statistics

const KT = KeplerTools

@testset "Orbits and Bodies" begin
    sun = Star("Kerbol", 261600000., 1.1723328e18);
    @test length(sun.satellite_bodies) == 0
    kerbin =    CelestialBody("Kerbin", 600000.,  3.5316000e12,       21549.425,      0.,     Orbit(13599840256., 0., 0., 0., 0., 3.14, 0., sun));
    @test length(sun.satellite_bodies) == 1
    moho =      CelestialBody("Moho",   250000.,  168.609378655e9,    1210000.,       190.,   Orbit(5263138304.,0.200000003, deg2rad(7.), deg2rad(70.), deg2rad(15.), 3.14, 0., sun));
    @test length(sun.satellite_bodies) == 2
    eve =       CelestialBody("Eve",    700000.,  8171.730229211e9,   80500.,         0.,     Orbit(9832684544.,0.01, deg2rad(2.099999905), deg2rad(15.000000000), deg2rad(0.), 3.14, 0., sun));
    @test length(sun.satellite_bodies) == 3
    duna =      CelestialBody("Duna",   320000.,  301.363211975e9,    65517.859375,   90.,    Orbit(20726155264., 0.051, deg2rad(0.06), deg2rad(135.5), 0., 3.14, 0., sun));
    @test length(sun.satellite_bodies) == 4
end

@testset "Kepler's Equation Solver" begin
    # elliptical case
    es = collect(range(0,stop=0.999,length=1000))
    Ms = collect(range(-π,stop=π,length=1000))
    errs = fill(NaN, length(es), length(Ms))
    for (i,e) in enumerate(es) for (j,M) in enumerate(Ms)
        θ = KT.mean_to_true(M, e)
        errs[i,j] = abs(KT.angledist(KT.true_to_mean(θ, e), M))
    end end
    @test maximum(errs) < 1e-9
    @test mean(errs) < 1e-12
    # hyperbolic case
    es = collect(range(1.01,stop=100,length=1000))
    Ms = collect(range(-100*π,stop=100*π,length=1000))
    errs = fill(NaN, length(es), length(Ms))
    for (i,e) in enumerate(es) for (j,M) in enumerate(Ms)
        θ = KT.mean_to_true(M, e)
        errs[i,j] = abs(KT.angledist(KT.true_to_mean(θ, e), M))
    end end
    @test maximum(errs) < 1e-9
    @test mean(errs) < 1e-12
end

@testset "Anomalies" begin
    # elliptical orbits from Kerbol system
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    # mean anomaly
    for bd in sun.satellite_bodies
        t = rand()*10*bd.orbit.period  - 5*bd.orbit.period + bd.orbit.epoch
        M = time_to_mean(t, bd.orbit)
        @test isapprox(t, mean_to_time(M, bd.orbit, t-1), rtol=1e-12)
    end
    # true anomaly
    for bd in sun.satellite_bodies
        t = rand()*10*bd.orbit.period  - 5*bd.orbit.period + bd.orbit.epoch
        θ = time_to_true(t, bd.orbit)
        @test isapprox(t, true_to_time(θ, bd.orbit, t-1), rtol=1e-12)
    end

    # random elliptical orbits
    for i in 1:50
        orb = Orbit(moho.orbit.a + 3*duna.orbit.a*rand(),
                    rand(),
                    π*rand(),
                    2π*rand(),
                    2π*rand(),
                    2π*rand(),
                    1e6*rand(),
                    sun
                    )
        t = rand()*10*orb.period  - 5*orb.period + orb.epoch
        M = time_to_mean(t, orb)
        @test isapprox(t, mean_to_time(M, orb, t-1), rtol=1e-12)
        θ = time_to_true(t, orb)
        @test isapprox(t, true_to_time(θ, orb, t-1), rtol=1e-12)
    end

    # random hyperbolic orbits
    for i in 1:50
        e = 1+4*rand()
        a = -(moho.orbit.a + 3*duna.orbit.a*rand())/(e-1)
        orb = Orbit(a,
                    e,
                    π*rand(),
                    2π*rand(),
                    2π*rand(),
                    2π*rand(),
                    1e6*rand(),
                    sun
                    )
        t = rand()*10*orb.period  - 5*orb.period + orb.epoch
        M = time_to_mean(t, orb)
        @test isapprox(t, mean_to_time(M, orb, t-1), rtol=1e-12)
        θ = time_to_true(t, orb)
        @test isapprox(t, true_to_time(θ, orb, t-1), rtol=1e-12)
    end

end

@testset "Bases and Rotations" begin
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    #TODO
end

@testset "State vectors" begin
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    # recover orbits from state vectors
    for bd in sun.satellite_bodies
        t = (rand()-2)*10*bd.orbit.period + bd.orbit.epoch
        pos, vel = time_state_vector(t, bd.orbit)
        orb = Orbit(t, pos, vel, bd.primary, bd.orbit.epoch)
        @test orb ≈ bd.orbit
        pos2, vel2 = time_state_vector(t, orb)
        @test isapprox(pos, pos2, atol=norm(pos)*1e-6)
        @test isapprox(vel, vel2, atol=norm(vel)*1e-6)
        orb2 = Orbit(t, pos2, vel2, bd.primary)
        @test orb2 ≈ bd.orbit
    end
    # generate random elliptical orbits for testing, and try to recover them from state vectors
    for i in 1:50
        orb = Orbit(moho.orbit.a + 3*duna.orbit.a*rand(),
                    rand(),
                    π*rand(),
                    2π*rand(),
                    2π*rand(),
                    2π*rand(),
                    1e6*rand(),
                    sun
                    )
        time = orb.epoch + 1e6*rand()
        pos, vel = time_state_vector(time, orb)
        @test Orbit(time, pos, vel, sun) ≈ orb
    end

    # generate random hyperbolic orbits for testing, and try to recover them from state vectors
    for i in 1:50
        e = 1+4*rand()
        a = -(moho.orbit.a + 3*duna.orbit.a*rand())/(e-1)
        orb = Orbit(a,
                    e,
                    π*rand(),
                    2π*rand(),
                    2π*rand(),
                    2π*rand(),
                    1e6*rand(),
                    sun
                    )
        time = orb.epoch + 1e6*rand()
        pos, vel = time_state_vector(time, orb)
        @test Orbit(time, pos, vel, sun) ≈ orb
    end
end

@testset "Plane angles" begin
    @test angle_in_plane([1,0,0], [1,0,0]) == 0.
    @test angle_in_plane([0,1,0], [1,0,0]) == π/2
    @test angle_in_plane([1,0,0], [0,1,0]) == 3π/2
    @test angle_in_plane([1,1,0], [1,0,0]) == π/4
    @test angle_in_plane([1,0,0], [1,1,0]) == 7π/4

    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    # phase angle between Kerbin and Duna for a typical transfer, from from alexmoon's app: https://alexmoon.github.io/ksp/
    @test isapprox(angle_in_plane(duna.orbit, kerbin.orbit, 5091552), deg2rad(36.49); rtol=0.001)

end

@testset "Lambert Solver" begin
    # test from Problems 5.3 and 5.4, http://www.braeunig.us/
    #  au = 1.495978707e11     # number of meters per astronomical unit
    μ = 3.964016e-14        # the real sun AU^3/m^2
    Δt = 207*24*60*60
    r̄ₛ = [0.473265, -0.899215, 0.] # .* au
    r̄ₑ = [0.066842, 1.561256, 0.030948] # .* au
    v̄ = KT.p_lambert_velocity(r̄ₛ, r̄ₑ, Δt, μ)
    @test isapprox(v̄, [0.000000193828, 0.000000101824, 0.00000000861759]; rtol=1e-3)

    # transfer from Kerbin to Duna, from alexmoon's app: https://alexmoon.github.io/ksp/
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    torb = KT.p_lambert(kerbin.orbit, duna.orbit, 5091522., 5588238. +5091522.)
    @test isapprox(torb.a, 1.68e10, rtol=0.01)
    @test isapprox(torb.e, 0.194, rtol=0.01)

    # generate random elliptical orbits for testing, and recover them with the Lambert solver
    μ = sun.μ
    for i in 1:50
        orb = Orbit(moho.orbit.a + 3*duna.orbit.a*rand(),
                    rand(),
                    π*rand(),
                    2π*rand(),
                    2π*rand(),
                    2π*rand(),
                    1e6*rand(),
                    sun
                    )
        starttime = orb.epoch + 1e6*rand()
        endtime = starttime + orb.period*rand()
        (spos, svel), (epos, evel) = time_state_vector(starttime, orb), time_state_vector(endtime, orb)
        orb2 = p_lambert(orb, orb, starttime, endtime)
        @test orb2 ≈ orb
        @test isapprox(time_orbital_position(starttime, orb2), spos)
        @test isapprox(time_orbital_velocity(starttime, orb2), svel)
        @test isapprox(time_orbital_position(endtime, orb2), epos)
        @test isapprox(time_orbital_velocity(endtime, orb2), evel)
        d = orb.i<π/2 ? 1 : -1
        @test isapprox(KT.p_lambert_velocity(spos, epos, endtime-starttime, μ; dir=d), svel, rtol=1e-6)
    end

    # generate random hyperbolic orbits for testing, and recover them with the Lambert solver
    μ = sun.μ
    for i in 1:50
        e = 1+4*rand()
        a = -(moho.orbit.a + 3*duna.orbit.a*rand())/(e-1)
        orb = Orbit(a,
                    e,
                    π*rand(),
                    2π*rand(),
                    2π*rand(),
                    2π*rand(),
                    1e6*rand(),
                    sun
                    )
        starttime = orb.epoch + 1e6*rand()
        endtime = starttime + orb.period*rand()
        (spos, svel), (epos, evel) = time_state_vector(starttime, orb), time_state_vector(endtime, orb)
        orb2 = p_lambert(orb, orb, starttime, endtime)
        @test orb2 ≈ orb
        @test isapprox(time_orbital_position(starttime, orb2), spos)
        @test isapprox(time_orbital_velocity(starttime, orb2), svel)
        @test isapprox(time_orbital_position(endtime, orb2), epos)
        @test isapprox(time_orbital_velocity(endtime, orb2), evel)
        d = orb.i<π/2 ? 1 : -1
        @test isapprox(KT.p_lambert_velocity(spos, epos, endtime-starttime, μ; dir=d), svel, rtol=1e-6)
    end
end

@testset "Departure and arrival trajectory calculation" begin
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    # define transfer trajectory
    starttime = 5091522.
    endtime = starttime + 5588238.
    torb = p_lambert(kerbin.orbit, duna.orbit, starttime, endtime)
    v̄rel = time_orbital_velocity(starttime, torb) - time_orbital_velocity(starttime, kerbin.orbit)

    # define parking orbit
    pkorb = Orbit(kerbin.eqradius+100000, kerbin)

    # get fast departure orbit (valid for circular parking orbits)
    dorb, Δv̄ = fast_departure_orbit(pkorb, v̄rel, starttime)
    @test isapprox(time_orbital_velocity(KT.ejection_time(dorb), dorb), v̄rel; rtol=1e-6)

    # recover (something similar to) the fast departure orbit from the optimizer
    θₒ = time_to_true(dorb.epoch, pkorb)
    d2orb, eccobj, Δv̄2 = departure_orbit(pkorb, v̄rel, starttime)
    @test isapprox(d2orb, dorb; atol=0.01, rtol=0.01)
end

@testset "Instantaneous Burns" begin
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    #TODO
end

@testset "Patched Conic Transfer Calculations" begin
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    # Transfer from Kerbin to Duna
    startorb = Orbit(kerbin.eqradius+100000, kerbin)
    endorb   = Orbit(duna.eqradius  +100000, duna)
    starttime = 5091522.
    endtime = starttime + 5588238.;
    tfer = Transfer(startorb, endorb, starttime, endtime)
    ftfer = fastTransfer(startorb, endorb, starttime, endtime)
    @test isapprox(tfer.Δv, ftfer.Δv, rtol=0.01)

    # Porkchop plot containing time ranges around above transfer
    startrange = [0.,               852*6*3600.]
    flightrange =[151 *6 *3600.,    453*6. *3600.]
    pc = Porkchop(startorb, endorb, startrange..., flightrange...)
    @test isapprox(minimum(pc.Δv), tfer.Δv; rtol=0.01)
    @test Transfer(startorb, endorb, KT.best_start_end_times(pc)...).Δv == minimum(pc.Δv)

    fpc = fastPorkchop(startorb, endorb, startrange..., flightrange...)
    @test isapprox(minimum(fpc.Δv), tfer.Δv; rtol=0.01)
    @test fastTransfer(startorb, endorb, KT.best_start_end_times(fpc)...).Δv == minimum(fpc.Δv)

    # optimize patch times and positions
    @test sum(norm, KT.patch_position_errors(tfer)) > 0
    tfer = match_patch_positions(tfer)

    @test sum(norm, KT.patch_position_errors(tfer)) < 1

end

@testset "SoI patch propagation" begin
    include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
    # TODO
end


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using BenchmarkTools
include("./basiccomputer.jl")

lines = readlines("day08input.txt")

function part1(lines)
    state = initializecomputer(lines)
    visited = BitSet()
    repeat_visit(state::ComputerState) = state.pc in visited
    execute(state, visited, repeat_visit)
    return state.accumulator
end

function swap_jmp_and_nop!(instr::Instruction)
    if instr.op == jmp
        instr.op = nop
    elseif instr.op == nop
        instr.op = jmp
    else
        throw(ErrorException("No jmp or nop found"))
    end
end

function part2(lines::AbstractArray{String})
    visited = BitSet()
    n = length(lines)
    n_next = n+1
    instructions = map(parseline, lines)
    state = ComputerState(0, 1, instructions)
    function repeat_visit_or_out_of_range(state::ComputerState)
        state.pc in visited || state.pc > n
    end
    for i in 1:n
        if instructions[i].op == acc
            continue
        end
        original_instr = instructions[i]
        swap_jmp_and_nop!(instructions[i])

        execute(state, visited, repeat_visit_or_out_of_range)

        # check for completion
        if state.pc == n_next
            return state.accumulator
        end

        # reset
        empty!(visited)
        swap_jmp_and_nop!(instructions[i])
        state.accumulator = 0
        state.pc = 1
    end
end

println("Part 1:")
answer1 = @btime part1(lines)
println(answer1)
println("Part 2:")
answer2 = @btime part2(lines)
println(answer2)

ext: jl


text:
#=  ADMM Algorithmic Regularization Paths for Sparse Statistical Machine Learning
#   Yue Hu, Eric C. Chi and Genevera I. Allen
#   arXiv:1504.06637v1
# 3.2 Reduced-Rank Multi-Task Learning
  Goal:  fit for B in Y = X*B+eps
  n # num iid samples
  p # num covariates
  q # num outcomes
# Algorithm 4 Algorithmic Regularization Path for Reduced-Rank Regression
  1. Initialize: Z0 = 0, U0 = 0, g0 = e, and step size t > 0.
  2. Precompute: H = (XTX=n+I)􀀀1 and HXTY.
  3. While kZkk 6= 0
     gk = gk􀀀1 +t (or gk = gk􀀀1t).
     Bk = HXTY+H(Zk􀀀1􀀀Uk􀀀1).
     Zk = SVTgk (Bk+Uk􀀀1). (Record Zk at each iteration.)
  Uk = Uk􀀀1+Bk􀀀Uk
  end
  4. Output fZk : k = 1;    ;Kg as the algorithmic regularization path.
=#

# Singular Value Treshold of A w/ constant gammma
function SVT(A::Array,gamma::Float64)
    fac = svdfact(A)
    fac[:U]*(diagm(max(0,fac[:S].-gamma)))*fac[:Vt]
end

function fit_arp_rrr{T<:Number}(X::Array{T,2}, Y::Array{T,2}; stepsize::Float64 =1.05, gamma_init::Float64 =0.01, max_itterations::Integer=200, min_active::Integer=0, verbose::Bool=false)
  @assert size(X,1) == size(Y,1)
  n = size(X,1)
  p = size(X,2)
  q = size(Y,2)
  gamma_list = zeros(T,max_itterations)
  chisq_list = zeros(T,max_itterations)
  Z_list = zeros(max_itterations,p,q)
  B = zeros(p,q)  # Preallocate working matrices
  Z = zeros(p,q)
  U = zeros(p,q)
  H = inv(X'*X/n+diagm(ones(p)))
  HXtY = H*X'*Y
  gamma = gamma_init
  last_itteration = max_itterations
  for k in 1:max_itterations
    gamma = gamma * stepsize
    B = HXtY+H*(Z-U)
    Z = SVT(B+U,gamma)
    gamma_list[k]= gamma
    Z_list[k,:,:] = Z
    U = U + B - Z
    num_active = length(find(Z))
    if verbose
      chisq = sumabs2(Y-X*Z)
      print("# ", k, " ", gamma,": ")
      print("Z norms= ", num_active, " ", sumabs(Z), " ", sumabs2(Z), " chi^2= ", chisq, " ")
      #print("B=",B," Z=",Z," U=",U," ")
      println("")
    end
    if sumabs2(Z) == 0. || num_active < min_active
      last_itteration = k
      break
    end
  end
  # resize to last_itteration?
  return (gamma_list,Z_list)
end

#= Test orig alg
 n = 25
 p = 4
 q = 3
 X = randn(n,p)
 B = make_B(p,q,frac_active = 0.1)
 Y = X*B + 0.01*randn(n,q)
 println("B norms: ", length(find(B)), "  ", sumabs(B), " ", sumabs2(B))

(gamma_list, Z_list)  = fit_arp_rrr(X,Y, max_itterations = 200)
 print_path(gamma_list, Z_list)
=#

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
function parsesentence(md, targetdir, notebookdir)
    return parsefigures(md, targetdir, notebookdir) |> parselinks |> parsebold |> parseitalics
end

function parsefigures(md, targetdir, notebookdir)
    parsedsentence = ""
    i = 1

    while i < length(md)
        pstart    = findnext("![", md, i)              !== nothing ? findnext("![", md, i)              : break
        pend      = findnext("]", md, pstart[end])    !== nothing ? findnext("]", md, pstart[end])    : break
        linkstart = findnext("(", md, pend[end])       !== nothing ? findnext("(", md, pend[end])       : break
        linkend   = findnext(")", md, linkstart[end]) !== nothing ? findnext(")", md, linkstart[end]) : break

        parsedsentence *= md[i:pstart[1] - 1]
        text = md[pstart[end] + 1:pend[1] - 1]
        link = md[linkstart[end] + 1:linkend[1] - 1]
        figurename = ""
        if link[end - 3:end] == ".svg"
            figuresvg = notebookdir * link
            figurename *= basename(figuresvg[1:end - 4])
            figurepdf = targetdir * "/figures/" * figurename * ".pdf"
            rsvg_convert() do cmd
               run(`$cmd $figuresvg -f pdf -o $figurepdf`)
           end
        else
            figurename *= basename(link)
            cp(notebookdir * link, targetdir * "/figures/" * figurename, force=true)
        end

        parsedsentence *= "\n\\begin{figure}[H]\n\t \\centering\n\t\\includegraphics[width=0.8\\textwidth]{./figures/" * figurename * "}\n\t\\caption{" * text * "}\n\t\\label{fig:" * figurename * "}\n\\end{figure}\n"
        i = linkend[end] + 1

    end
    parsedsentence *= md[i:end]
    return parsedsentence
end

function parselinks(md)
    parsedsentence = ""
    i = 1

    while i < length(md)
        pstart    = findnext("[", md, i)              !== nothing ? findnext("[", md, i)              : break
        pend      = findnext("]", md, pstart[end])    !== nothing ? findnext("]", md, pstart[end])    : break
        linkstart = findnext("(", md, pend[end])       !== nothing ? findnext("(", md, pend[end])       : break
        linkend   = findnext(")", md, linkstart[end]) !== nothing ? findnext(")", md, linkstart[end]) : break

        parsedsentence *= md[i:pstart[1] - 1]
        text = md[pstart[end] + 1:pend[1] - 1]
    link = md[linkstart[end] + 1:linkend[1] - 1]
        parsedsentence *= " \\href{" * link * "}{" * text * "}"
    i = linkend[end] + 1
    end
    parsedsentence *= md[i:end]
    return parsedsentence
end

function parseitalics(md)
    parsedsentence = ""
    i = 1
    while i < length(md)
        pstart    = findnext("*", md, i) !== nothing ? findnext("*", md, i)               : break
        pend      = findnext("*", md, pstart[end] + 1) !== nothing ? findnext("*", md, pstart[end] + 1)    : break
        parsedsentence *= md[i:pstart[1] - 1]
text = md[pstart[end] + 1:pend[1] - 1]
        parsedsentence *= "\\textit{" * text * "}"
        i = pend[end] + 1
    end
    parsedsentence *= md[i:end]
    return parsedsentence
end

function parsebold(md)
    parsedsentence = ""
    i = 1
    while i < length(md)
        pstart    = findnext("**", md, i) !== nothing ? findnext("**", md, i)               : break
        pend      = findnext("**", md, pstart[end]) !== nothing ? findnext("**", md, pstart[end])    : break
        parsedsentence *= md[i:pstart[1] - 1]
        text = md[pstart[end] + 1:pend[1] - 1]
        parsedsentence *= "\\textbf{" * text * "}"
        i = pend[end] + 1
        end
    parsedsentence *= md[i:end]
    return parsedsentence
end

function parseparagraph(paragraph, targetdir, notebookdir)
    parsedparagraph = ""
    for (i, sentence) in enumerate(split(paragraph, "\$"))
        if iseven(i)
            # It's a math sentence, i.e text is between $ $, so nothing should be parsed
            parsedparagraph *= "\$" * sentence * "\$"
                    else
            for (i, subsentence) in enumerate(split(sentence, "`"))
                if iseven(i)
                    # It's a code sentence, i.e text is between ` ` and is not a math sentence.
                    parsedparagraph *= "\\lstinline[style=julia]{" * subsentence * "}"
                else
                    parsedparagraph *= parsesentence(subsentence, targetdir, notebookdir)
                end
            end
        end
    end
    return parsedparagraph
end

"""
    configparser(config=nothing, template=:book)
This function contains the configuration to be used in the
Markdown to LaTeX parser.
"""
# function configparser(config=nothing, template=:book)
    # Note that the order is important! It defines the order
    # for the parser evaluation.
    # parserdictionary =
# end

function markdowntolatex(md, targetdir, notebookdir; template=:book)
    tag = false
    component = ""
    parsedtext = ""
    for line in split(md, "\n")
        l = strip(line)
        if l == ""
            l = "\n\n"
        end

        if startswith(l, "```math")
            tag = !tag
            parsedtext *= "\n\\begin{displaymath}\n"
            component = "displaymath"
            continue
        elseif startswith(l, "```") && tag
            tag = !tag
            parsedtext *= "\\end{" * component * "}\n"
            continue

        elseif startswith(l, "\$\$")
            if !tag
                tag = !tag
                parsedtext *= "\n\\begin{displaymath}\n"
            else
                tag = !tag
                parsedtext *= "\\end{displaymath}\n"
            end
            continue
        elseif startswith(l, "####")
            if template == :article || template == :matharticle
                parsedtext *= "\n\\textbf{" * l[6:end] * "}\n"
            else
                parsedtext *= "\n\\subsubsection{" * l[6:end] * "}\n"
            end
            continue
        elseif startswith(l, "###")
            if template == :article || template == :matharticle
                parsedtext *= "\n\\subsubsection{" * l[5:end] * "}\n"
            else
                parsedtext *= "\n\\subsection{" * l[5:end] * "}\n"
            end
            continue
        elseif startswith(l, "##")
            if template == :article || template == :matharticle
                parsedtext *= "\n\\subsection{" * l[4:end] * "}\n"
            else
                parsedtext *= "\n\\section{" * l[4:end] * "}\n"
            end
            continue
        elseif startswith(l, "#")
            if template == :article || template == :matharticle
                parsedtext *= "\n\\section{" * l[3:end] * "}\n"
            else
                parsedtext *= "\n\\chapter{" * l[3:end] * "}\n"
            end
            continue
    end

        if tag
            parsedtext *= "\t" * l * "\n"
        elseif !tag
            parsedtext *= parseparagraph(l, targetdir, notebookdir)
        end
    end
    return parsedtext
end


text:
using CompScienceMeshes, BEAST, Base.Test

fn = joinpath(dirname(@__FILE__),"assets","sphere35.in")

sol = 1.0
D, Δx = 1.0, 0.35
#Γ = meshsphere(D, Δx)
Γ = readmesh(fn)
X = raviartthomas(Γ)
Xm = subset(X,[1])
Xn = subset(X,[100])
Δt, Nt = 0.13/sol, 40
T = timebasisshiftedlagrange(Δt, Nt, 3)
U = timebasisdelta(Δt, Nt)

W = Xm ⊗ U
V = Xn ⊗ T

t = MWSingleLayerTDIO(sol,-1/sol,-sol,2,0)
Z1 = assemble(t, W, V)

## Run again, this time scales with another speedOfLight

sol = 3.0
D, Δx = 1.0, 0.35
#Γ = meshsphere(D, Δx)
Γ = readmesh(fn)
X = raviartthomas(Γ)
Xm = subset(X,[1])
Xn = subset(X,[100])
Δt, Nt = 0.13/sol, 40
T = timebasisshiftedlagrange(Δt, Nt, 3)
U = timebasisdelta(Δt, Nt)

W = Xm ⊗ U
V = Xn ⊗ T

t = MWSingleLayerTDIO(sol,-1/sol,-sol,2,0)
Z3 = assemble(t, W, V)

m = n = 1
@test find(Z1[m,n,:]) == find(Z3[m,n,:])

I = find(Z1[m,n,:])
@test all(sol*Z1[m,n,I] .≈ Z3[m,n,I])


text:
using Distributed
nprocs()==1 && addprocs()
using DISSEQT
using DISSEQT.AlignUtils
using SynapseClient
using DISSEQT.SynapseTools
using JLD
using Dates

function main()
# --- Synapse login ------------------------------------------------------------
    syn = SynapseClient.login()
    isdir("synapsecache") || mkdir("synapsecache")

# --- Setup --------------------------------------------------------------------
    # Set clean=true to rerun alignment. For clean=false, alignment will only be run if the log file is missing (i.e. no previous alignment was done).
    clean = false

    projectFolder = "syn11639899" # VignuzziLabPublic/Projects/FitnessLandscapes
    alignmentFolder = getchildbyname(syn, projectFolder, "Analysis", "Alignment")

    # Set uploadPath="synapseID" to upload files after running. Should point to "MyProject/Analysis/Alignment" folder. Set upload=nothing to skip uploading.
    uploadPath = alignmentFolder

    # Where to find sample .fastq files. Synapse Folder ID or local folder. Synapse ID should point to "MyProject/Raw Data/Sequencing".
    fastqPath = getchildbyname(syn, projectFolder, "Raw Data", "Sequencing")

    # Name of the run. (Should corrsespond to a subfolder of fastqPath.)
    runName = "H03UAAFXX"

    # This prefix will be added to all sample names. Normally same as runName. Set to "" if the .fastq files already have this prefix.
    namePrefix = runName

    # Where to find reference genomes. Synapse Folder ID or local folder. [Reference genomes will be uploaded if local.]
    # Synapse ID should point to "MyProject/Analysis/Alignment/ReferenceGenomes".
    referenceFolder = getchildbyname(syn, alignmentFolder, "ReferenceGenomes")

    # Rules for matching sample IDs to references.
    # The rule can either be a Regex or a function taking the sample ID and returning true/false.
    refs = [(r"-A1_", "Stop.fasta"), (r"-A2_", "Stop.fasta"), (r"-A3_", "Stop.fasta"), (r"-A4_", "Stop.fasta"), (r"-A5_", "WT.fasta"), (r"-A6_", "More.fasta"), (r"-A7_", "Less.fasta"), (r"-A8_", "Stop.fasta"), (r"-A9_", "WT.fasta"), (r"-A10_", "More.fasta"), (r"-A11_", "Less.fasta"), (r"-A12_", "Stop.fasta"), (r"-B1_", "Stop.fasta"), (r"-B2_", "Stop.fasta"), (r"-B3_", "Stop.fasta"), (r"-B4_", "Stop.fasta"), (r"-B5_", "WT.fasta"), (r"-B6_", "More.fasta"), (r"-B7_", "Less.fasta"), (r"-B8_", "Stop.fasta"), (r"-B9_", "WT.fasta"), (r"-B10_", "More.fasta"), (r"-B11_", "Less.fasta"), (r"-B12_", "Stop.fasta"), (r"-C1_", "Stop.fasta"), (r"-C2_", "Stop.fasta"), (r"-C3_", "Stop.fasta"), (r"-C4_", "Stop.fasta"), (r"-C5_", "WT.fasta"), (r"-C6_", "More.fasta"), (r"-C7_", "Less.fasta"), (r"-C8_", "Stop.fasta"), (r"-C9_", "WT.fasta"), (r"-C10_", "More.fasta"), (r"-C11_", "Less.fasta"), (r"-C12_", "Stop.fasta"), (r"-D1_", "Stop.fasta"), (r"-D2_", "Stop.fasta"), (r"-D3_", "Stop.fasta"), (r"-D4_", "Stop.fasta"), (r"-D5_", "WT.fasta"), (r"-D6_", "More.fasta"), (r"-D7_", "Less.fasta"), (r"-D8_", "Stop.fasta"), (r"-D9_", "WT.fasta"), (r"-D10_", "More.fasta"), (r"-D11_", "Less.fasta"), (r"-D12_", "Stop.fasta"), (r"-E1_", "Stop.fasta"), (r"-E2_", "Stop.fasta"), (r"-E3_", "Stop.fasta"), (r"-E4_", "Stop.fasta"), (r"-E5_", "WT.fasta"), (r"-E6_", "More.fasta"), (r"-E7_", "Less.fasta"), (r"-E8_", "Stop.fasta"), (r"-E9_", "WT.fasta"), (r"-E10_", "More.fasta"), (r"-E11_", "Less.fasta"), (r"-E12_", "Stop.fasta"), (r"-F1_", "Stop.fasta"), (r"-F2_", "Stop.fasta"), (r"-F3_", "Stop.fasta"), (r"-F4_", "Stop.fasta"), (r"-F5_", "WT.fasta"), (r"-F6_", "More.fasta"), (r"-F7_", "Less.fasta"), (r"-F8_", "Stop.fasta"), (r"-G1_", "Stop.fasta"), (r"-G2_", "Stop.fasta"), (r"-G4_", "Stop.fasta"), (r"-G5_", "WT.fasta"), (r"-G6_", "More.fasta"), (r"-G7_", "Less.fasta"), (r"-G8_", "Stop.fasta"), (r"-H1_", "Stop.fasta"), (r"-H2_", "Stop.fasta"), (r"-H3_", "Stop.fasta"), (r"-H4_", "Stop.fasta"), (r"-H5_", "WT.fasta"), (r"-H6_", "More.fasta"), (r"-H7_", "Less.fasta"), (r"-H8_", "Stop.fasta"),
            (r"Undetermined", "phiX174.fasta")] # Always keep this unless PhiX wasn't used in the sequencing.
    refs = getreferenceinfo(syn, refs, referenceFolder) # Get path/synapse id and local path info.

    # File with adapters. Synapse File ID or local file. [Uploaded if local.]
    adapters = getchildbyname(syn, alignmentFolder, "Adapters", "adapters.fa")

    # Local logFile.
    logFile = "AlignUtils.log"

    # Local file with sample info. Required for Uploading of a previously run alignment.
    sampleInfoFile = "sampleInfo.jld"

    # Minimum support for updating consensus during alignment. (Otherwise fallback to reference sequence.)
    consensusMinSupport = 100

    # Minimum support for updating consensus with indels during alignment. (Otherwise fallback to reference sequence.)
    consensusIndelMinSupport = 1000

    # Maximumum number of times a sample will be aligned with updated consensus.
    maxAlignIterations = 10

# --- Cleanup ------------------------------------------------------------------
    clean = clean || !isfile(logFile)

    if clean
        makecleanfolder("bam") || return
        makecleanfolder("temp") || return
    else
        println("Skipping Alignment [clean=false]")
    end

# --- Alignment ----------------------------------------------------------------

    samples = Sample[]
    if clean
        log = open(logFile,"w")
        isfile(sampleInfoFile) && rm(sampleInfoFile)

        start = time()
        println(log,"Starting alignment batch run at $(now())"); flush(log)
        samples = find_samples(syn, fastqPath, runName, namePrefix, log=log); flush(log)
        assign_reference!(samples,refs,log=log); flush(log)

        align_samples!(syn,samples,adapters,"bam","temp",log,consensusMinSupport=consensusMinSupport,consensusIndelMinSupport=consensusIndelMinSupport,maxAlignIterations=maxAlignIterations); flush(log)
        reference_sanity_check(samples, refs, log=log)

        duration = time()-start
        println(log,"Finished alignment batch run in $(duration)s."); flush(log)

        close(log)

        save(sampleInfoFile, "samples", samples, compress=true)
    end

# --- Upload -------------------------------------------------------------------
    if uploadPath != nothing
        if !clean # reload sample info if necessary
            try
                samples = load(sampleInfoFile, "samples")
            catch
                println("Sample info file is missing or corrupt. Please rerun alignment.")
            end
        end

        refPaths = [r[3] for r in refs]
        uploadaligned(syn, uploadPath, runName, @__FILE__, logFile, adapters, refPaths, samples)
    end

end

main()


text:
import BoteSalvatICX # For ionization crosssections

struct Bote2009 <: NeXLAlgorithm end

"""
    ionizationcrosssection(z::Int, shell::Int, energy::AbstractFloat, ::Type{Bote2009})

Compute the absolute ionization crosssection (in cm²) for the specified element, shell and
electon energy (in eV).
"""
ionizationcrosssection(z::Int, ss::Int, energy::AbstractFloat, ::Type{Bote2009}) =
    BoteSalvatICX.hasedge(z, ss) ?
    BoteSalvatICX.ionizationcrosssection(
        z,
        ss,
        energy,
        BoteSalvatICX.edgeenergy(z, ss),
    ) : 0.0


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libClp"], :libClp),
    LibraryProduct(prefix, ["libOsiClp"], :libOsiClp),
    LibraryProduct(prefix, ["libClpSolver"], :libClpSolver),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/Clp_jll.jl/releases/download/Clp-v1.17.6+4"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-gnu-libgfortran3-cxx03.tar.gz", "059439acc52693be3f02fa5af4a64178d29b8682df25371884bc19be19b9e786"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-gnu-libgfortran3-cxx11.tar.gz", "39d1b7528b328ddb3ab59ba28da8e0eccc9c9ae5af8ea07c45cf137b16705cef"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-gnu-libgfortran4-cxx03.tar.gz", "a98d527a2d0a29e8f2ea4158ab3ed1ffc825286427328cad1a7711ec441cc92d"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-gnu-libgfortran4-cxx11.tar.gz", "f435877ef8863296760285027a79e8d061836d21389d36f3713a1fa512793757"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-gnu-libgfortran5-cxx03.tar.gz", "ae7515a6849cff6818e5c0819ae8320535c5fcc797aca75ede4dd0486c7117dc"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-gnu-libgfortran5-cxx11.tar.gz", "3b46bacb16a3f561fbf47ce7079951b8a1cb003b750de72199568a168a5dfbf5"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-musl-libgfortran3-cxx03.tar.gz", "04f7518752e7195d7d725e19ee47e9d0184ae9273880424e60ff385bc6b94599"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-musl-libgfortran3-cxx11.tar.gz", "1583ea3f2567c3d88b5680b98901f75e0ce61af957c0c78c033fbc641649af88"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-musl-libgfortran4-cxx03.tar.gz", "62895d58c94272cc5060d0837ab78fc31f16e87629d5f7e95569c810af2d780a"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-musl-libgfortran4-cxx11.tar.gz", "39b4de5459d304e1db21817d3d90c85edb20630d5c6aa7c05b2a0b1fb4c4540d"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-musl-libgfortran5-cxx03.tar.gz", "cc4232538b8faf5a045ea397ed2e28951b6c1b7e74e030c0c9d79c20763ed560"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.aarch64-linux-musl-libgfortran5-cxx11.tar.gz", "b11bed9fb7f537b88040d68c4ce49fb5c83195e56ad6d466def251f992a292f3"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-gnueabihf-libgfortran3-cxx03.tar.gz", "30969dc5bac04677737167372e205586b154fe38e458ccba051f2ba4333e05b8"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-gnueabihf-libgfortran3-cxx11.tar.gz", "e28cc87c9c6b4c1b5a1161654505ad05c9048bf6e50699a477f56aa202cb1bab"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-gnueabihf-libgfortran4-cxx03.tar.gz", "d33c1e4c12af9ccec3e636f3b29d7538a7eae71e0d71cedffbec51e552165f20"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-gnueabihf-libgfortran4-cxx11.tar.gz", "7370687426d753d516e9855151e42b5a03c741e1bab7086ca63b033ee7abc16d"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-gnueabihf-libgfortran5-cxx03.tar.gz", "fa76a2342cae9d0b795d89b11f0475f7cce8726c26b3d95e86567ae711e803ee"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-gnueabihf-libgfortran5-cxx11.tar.gz", "e491590dbfb93578df18727d75b232c6d3fd183129e227ecb9d5b95b1487277c"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-musleabihf-libgfortran3-cxx03.tar.gz", "3676265097ab9c2e7f4acc66d72d0adc6c5035fb5c5c74245553e0f971934465"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-musleabihf-libgfortran3-cxx11.tar.gz", "673a83dd7c4416552fa735cc74253d63d47b0bf553818a3ad4af606ba5ee827d"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-musleabihf-libgfortran4-cxx03.tar.gz", "08acd28c93cef77b92c8624e716b98f4fb1fde49a5a9a1948eb63e11426079c9"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-musleabihf-libgfortran4-cxx11.tar.gz", "51fd91bb457b4383892f8b54cf70c1a337d542ae93853f51c2ec6debcd42bb43"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-musleabihf-libgfortran5-cxx03.tar.gz", "71c72009b18c5585dbcb4c77b54e1f64e10cd87d07924bd4e4d46d20ce3f1f94"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.armv7l-linux-musleabihf-libgfortran5-cxx11.tar.gz", "e88098e676b92a4feaee97cfb25ddc949c7c97f137420d52055089345363c4ba"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-gnu-libgfortran3-cxx03.tar.gz", "4bf40f8a12b821d6edc86a7dc55431c1ee73ec56d2875ef41c7b974bf6b85c6d"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-gnu-libgfortran3-cxx11.tar.gz", "809012f81214a22e5e6a6fd00b7f4eb3699a6965c7e87abe9fc0df47cf367e7d"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-gnu-libgfortran4-cxx03.tar.gz", "b0d35e6b29ed5c15bf62b6f2bd07c5911eb998f434919611253b8a8eac5dd4c2"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-gnu-libgfortran4-cxx11.tar.gz", "ec8324c8b230ed7befd57346f2285932906d927e874b403b7d87c1bd32fef9ea"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-gnu-libgfortran5-cxx03.tar.gz", "73f6bde94a838d595a70f5fb1d06485afe0adfcd248c119c98dbb64129af8fdb"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-gnu-libgfortran5-cxx11.tar.gz", "5dad9d1fab13faa83ebdd7bb7d690952a09ddb46cee428185a7abe8566ea3fad"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-musl-libgfortran3-cxx03.tar.gz", "ef21bfebbbeca91f0c3feac91b47ea087a5e229c0867217137ad3248e804004a"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-musl-libgfortran3-cxx11.tar.gz", "24e55b07918952a4320c08b2a2fff96625984438ebadcb6a6515a845fff1c5bb"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-musl-libgfortran4-cxx03.tar.gz", "6c1625eb19d2dd5f9e63c702a611520a11220d5e5448cf2e9c4d1af39596dc0e"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-musl-libgfortran4-cxx11.tar.gz", "29141af307f06d55c4646cad74d64644edc7185768a3d33cae647464bc04cf2f"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-musl-libgfortran5-cxx03.tar.gz", "f4c33d89f50a7287675165122e9c69289e2599cf78bcbe3cb041817ab2c50884"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-linux-musl-libgfortran5-cxx11.tar.gz", "271d40ca1f82e79a69f18526acb4eec12c18a12a1f6cde0b573f254e73c67f73"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-w64-mingw32-libgfortran3-cxx03.tar.gz", "9a18a4440afc4bfe501a5a14ae13a3d6a99d5c0c3974aec887edd172d1e6382d"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-w64-mingw32-libgfortran3-cxx11.tar.gz", "0129cfdb0018918a4aed88d5e4b7512fafc2ba7e837fb01ac7c1882611927ecc"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-w64-mingw32-libgfortran4-cxx03.tar.gz", "42a0314dc96d0f22c9d18c549a698c56b3b5caabf1c88390a12b2552c38c63b3"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-w64-mingw32-libgfortran4-cxx11.tar.gz", "450d21545f3b860ed05f636d122d9751992d94466c352ea4f86a86454db58cbc"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.i686-w64-mingw32-libgfortran5-cxx03.tar.gz", "e5151713849c9b274442dfd5585cfe38ead5706f764516f373c4b69aaf992a9a"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.i686-w64-mingw32-libgfortran5-cxx11.tar.gz", "ce6663cfb8cf471f119c7b01ca311656f09c710c5850795dc5b72d67422944c6"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-apple-darwin14-libgfortran3-cxx03.tar.gz", "96a580aa70f931196779d45722ab5aca18ee07b7c4944952550b600d253b0b4d"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-apple-darwin14-libgfortran3-cxx11.tar.gz", "9a5c31a93e87b804e2712a90ec15a4913fb126235445fa09b8d6047b74c3f364"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-apple-darwin14-libgfortran4-cxx03.tar.gz", "f9ca0649a69b0b4921db84a83dbc56bcd161842bac81eb6f31ce15e919073347"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-apple-darwin14-libgfortran4-cxx11.tar.gz", "53bcfb7ba0c9a5a5f45cbc033d09d10bf542559f2501167ad224a6aaf9b8c5c0"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-apple-darwin14-libgfortran5-cxx03.tar.gz", "4e5761e0edcddef0f32690451b0ee047e3f53095f7d12a5951fae29549b1926d"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-apple-darwin14-libgfortran5-cxx11.tar.gz", "5f16402c0f0451aa02aae7b131b6043819f29059d5be13c7021c4716b9e64fb5"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-gnu-libgfortran3-cxx03.tar.gz", "9985690f528d2b22e31184a725e643b808fb5d3a275d320b726a3baf70105eac"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-gnu-libgfortran3-cxx11.tar.gz", "991531cce8c7f62289f36aab551cc3318fa51867e3d5250a5b727c2a943be3cb"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-gnu-libgfortran4-cxx03.tar.gz", "1b158ee674aa6c2545eeed2dc301dc9a44efd1a3bc3a14d6be8886ff8ad1429f"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-gnu-libgfortran4-cxx11.tar.gz", "b47ff51552a4701f8cbbf29c8fea92a32a1ea279e9bb55973061741619d1c20e"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-gnu-libgfortran5-cxx03.tar.gz", "85e40e2e121f8d0c71a4cdc625b29348ff32457845c06ca8ccad9e3cf34569cf"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-gnu-libgfortran5-cxx11.tar.gz", "b3c7abe76b820bac309f89c437e422117e37cc75b2c6d518b796e03892238318"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-musl-libgfortran3-cxx03.tar.gz", "73acc89560dc5664444bd8f38f6173ae32f320d4893c6092503648a013075f0a"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-musl-libgfortran3-cxx11.tar.gz", "92c66b1014a4fa0f42e1e11568f9ac0ea592b8406f8df163ffbed17fd2ecf8f8"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-musl-libgfortran4-cxx03.tar.gz", "14334725b48f7c13cdde4f7c3a3adfa38b15b561e1ad9e6f1667729582b26af2"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-musl-libgfortran4-cxx11.tar.gz", "33fa27d4f91cda9b6913ff3d3eb591186c45249433b280f5bc523e0a98664be5"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-musl-libgfortran5-cxx03.tar.gz", "a437d3e04cb8cc5f672290e51b670d6f20f39086acbe87f48c3d236bff39c96b"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-linux-musl-libgfortran5-cxx11.tar.gz", "68edc01cd657d9ff979ed2329f2a6459dd69c4bc363f0305ecebf0539f3ba016"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-w64-mingw32-libgfortran3-cxx03.tar.gz", "f8f5e919145542e67ec5a38bc84f7eacff91f4a427079540eb87d74b6109884f"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-w64-mingw32-libgfortran3-cxx11.tar.gz", "ad8a66ecb1e1a19f33dc35346dbb9b9833b03c10a2002a5b9234a886e8e33348"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-w64-mingw32-libgfortran4-cxx03.tar.gz", "44b8d6c8d08f98274364f7008af4375b1a3c56a212792c49cfa6b44b17abfb95"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-w64-mingw32-libgfortran4-cxx11.tar.gz", "9385fbf05223267b36a3b5607aaac5bb87c378a79a7bdf769909f9603049c752"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Clp.v1.17.6.x86_64-w64-mingw32-libgfortran5-cxx03.tar.gz", "c5f936529b68613c72b6df816f305c98dbf4370dfcf82dc0c1ae3df2a0a40b12"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Clp.v1.17.6.x86_64-w64-mingw32-libgfortran5-cxx11.tar.gz", "d41cbab7481daf57a475dc445192085737e4768ff702402ba70c009e1bf1cd7a"),
)

# Install unsatisfied or updated dependencies:
# We added `, isolate=true` as otherwise, it would segfault when closing `OpenBLAS32`,
# probably because it is conflicting with Julia openblas.
unsatisfied = any(!satisfied(p; verbose=verbose, isolate=true) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end


text:
mutable struct ArrayBasedList{T} <: AbstractList{T}
    maxsize::Int
    size::Int
    curr::Int
    data::Vector{T}
    function ArrayBasedList{T}(maxsize::Int) where T
        size = 0
        curr = 1
        new(maxsize, size, curr, Vector{T}(undef, maxsize))
    end
end

ArrayBasedList{T}() where T = ArrayBasedList{T}(DEFAULT_MAX_SIZE)

isempty(l::ArrayBasedList) = l.size == 0
length(l::ArrayBasedList) = l.size

function insert!(l::ArrayBasedList, item)
    (l.size < l.maxsize) || throw(ErrorException("list capacity exceeded"))
    l.size += 1
    for i = l.size:-1:l.curr+1
        l.data[i] = l.data[i-1]
    end
    l.data[l.curr] = item
end

function append!(l::ArrayBasedList, item)
    (l.size < l.maxsize) || throw(ErrorException("list capacity exceeded"))
    l.size += 1
    l.data[l.size] = item
end

function remove!(l::ArrayBasedList)
    (1 <= l.curr <= l.size) || throw(ArgumentError("no current element"))
    # TODO unhandled case: movetoend -> remove -> remove
    item = l.data[l.curr]
    for i = l.curr:l.size-1
        l.data[i] = l.data[i+1]
    end
    l.size -= 1
    return item
end

function movetostart!(l::ArrayBasedList)
    l.curr = 1
end

function movetoend!(l::ArrayBasedList)
    l.curr = l.size + 1
end

function movetoposition!(l::ArrayBasedList, pos)
    (1 <= pos <= l.size+1) || throw(BoundsError(l, pos))
    l.curr = pos
end

prev!(l::ArrayBasedList) = (l.curr > 1) ? l.curr -= 1 : nothing
next!(l::ArrayBasedList) = (l.curr <= l.size) ? l.curr += 1 : nothing
currentposition(l::ArrayBasedList) = l.curr
function getvalue(l::ArrayBasedList)
    (1 <= l.curr <= l.size) || throw(ArgumentError("no current element"))
    return l.data[l.curr]
end


text:
using Pkg
Pkg.update(update_registry=false)

using Artifacts
using NixSourcerer
using Downloads
using CxxWrap.CxxWrapCore: libcxxwrap_julia_jll

function main()
    cd(@__DIR__) do
        meta = artifact_meta("libcxxwrap_julia", find_artifacts_toml(pkgdir(libcxxwrap_julia_jll)))
        url = nothing
        for dl in meta["download"]
            try
                # Test URL
                Downloads.download(dl["url"])
                url = dl["url"]
                break
            catch e
                @error e
            end
        end
        if url !== nothing
            p = NixSourcerer.read_project(joinpath(@__DIR__, "NixProject.toml"))
            p.specs["JlCxx"] = Dict("type" => "archive", "url" => url, "extraArgs" => Dict("stripRoot" => false))
            NixSourcerer.write_project(p, joinpath(@__DIR__, "NixProject.toml"))
            NixSourcerer.update_package()
        else
            error("Update failed")
        end
    end
end

main()


text:
using Weber
using Base.Test
include("find_timing.jl")

function cause_column_error()
  find_timing() do
    addtrial(moment(record,"test",unspecified_column_name="string"))
  end
end

const reserved_columns = [:reserved, :weber_version, :start_date,
                          :start_time, :offset, :trial, :time]

function cause_reserved_error(decl,use=decl)
  find_timing(columns=[decl]) do
    addtrial(moment(record,"test";[use => "test"]...))
  end
  nothing
end

@testset "Record Columns" begin
  @test_throws ErrorException cause_column_error()
  @test_throws ErrorException cause_reserved_error(:weber_version)
  @test_throws ErrorException cause_reserved_error(:start_date)
  @test_throws ErrorException cause_reserved_error(:start_time)
  @test_throws ErrorException cause_reserved_error(:offset)
  @test_throws ErrorException cause_reserved_error(:trial)
  @test_throws ErrorException cause_reserved_error(:time)
  @test_throws ErrorException cause_reserved_error(:joe,:bob)
end


text:
@testset "ClosestVector.jl" begin
  #=
  The quadratic triple QT = [Q, L, c] is an n-dimensional ellipsoid with rational coefficients.
  Where,
  Q is an nxn symmetric matrix of the quadratic form ,
  L is a rational column vector of length n &
  c is a rational number
  =#
  Q = matrix(QQ,3,3,[1 0 0; 0 1 0; 0 0 1]);
  L = matrix(fmpq[1,0,0]);
  c = fmpq(2//5);
  Lat = Zlattice(gram = matrix(QQ,3,3,[1,0,0,0,1,0,0,0,1]));
  v = matrix(fmpq[-1,0,0]);
  u = fmpq(3//5);
  @test Hecke.closest_vectors(Q, L,c, sorting=true)[1] == Hecke.closest_vectors(Lat, v, u, sorting=true)[1]

  function compare_functions(Q::MatrixElem{T}, K::MatrixElem{T}, d::T) where T <: RingElem
    List1 = Hecke.closest_vectors(Q, K, d, sorting=true);
    L, v, c = Hecke.convert_type(Q, K, d)
    List2 = Hecke.closest_vectors(L, v, c, sorting=true);
    for i in 1:length(List1)
      if List1[i] == List2[i]
        return true
      else
        error("The entry at position $i in the two lists don't match.")
      end
    end
  end

  function compare_functions(L::ZLat, v::MatrixElem{T} , upperbound::T) where T <: RingElem
    List1 = Hecke.closest_vectors(L, v, upperbound, sorting=true);
    Q, K, d = Hecke.convert_type(L, v, upperbound);
    List2 = Hecke.closest_vectors(Q, K, d, sorting=true);
    for i in 1:length(List1)
      if List1[i] == List2[i]
        return true
      else
        error("The entry at position $i in the two lists don't match.")
      end
    end
  end

  #EXAMPLE 1a: 3-dimensional sphere
  Q1 = matrix(QQ,3,3,[1 0 0; 0 1 0; 0 0 1]);
  L1 = matrix(fmpq[1,1,1]);
  c1 = fmpq(1);

  #EXAMPLE 1b: 4-dimensional sphere
  Q2 = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
  L2 = matrix(fmpq[1,1,1,1]);
  c2 = fmpq(3);

  #EXAMPLE 1c: 6-dimensional sphere
  Q3 = matrix(QQ,6,6,[1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1]);
  L3 = matrix(fmpq[1,0,1,0,2,0]);
  c3 = fmpq(4//3);

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

  #Examples for closest_vectors()
  #Example 1:
  Lat1 = Zlattice(gram = matrix(QQ,3,3,[1,0,0,0,1,0,0,0,1]));
  v1 = matrix(fmpq[-1,0,0]);
  u1 = fmpq(3//5);

  #Example 1:
  Lat2 = Zlattice(gram = matrix(QQ,4,4,[1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1]));
  v2 = matrix(fmpq[-1,-1,-1,-1]);
  u2 = fmpq(41//11);

  #Example 1:
  Lat3 = Zlattice(gram = matrix(QQ,6,6,[1,0,0,0,0,0, 0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1]));
  v3 = matrix(fmpq[-1,0,-1,0,-2,0]);
  u3 = fmpq(14//3);

  @testset "Comparing closest_vectors(Q, K, d) to closest_vectors(Lat, vec, upperbound) list" begin
    @test compare_functions(Q1, L1, c1) == true
    @test compare_functions(Q2, L2, c2) == true
    @test compare_functions(Q3, L3, c3) == true
  end

  @testset "Comparing closest_vectors(Lat, vec, upperbound) list to closest_vectors(Q, K, d)" begin
    @test compare_functions(Lat1, v1, u1) == true
    @test compare_functions(Lat2, v2, u2) == true
    @test compare_functions(Lat3, v3, u3) == true
  end
  @testset "Comparing closest_vectors() computed with respect to type fmpz and fmpq" begin
    Q0 = matrix(ZZ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
    L0 = matrix(fmpz[1,1,1,1]);
    c0 = fmpz(3);
    @test compare_functions(Q0, L0, c0) == true

    Q01 = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
    L01 = matrix(fmpq[1,1,1,1]);
    c01 = fmpq(3);
    @test compare_functions(Q01, L01, c01) == true

    List1 = Hecke.closest_vectors(Q0, L0, c0)
    List2 = Hecke.closest_vectors(Q01, L01, c01)
    for i in 1:length(List1)
      @test List1[i] == List2[i]
    end
  end

  @testset "Conversion from Quadratic Triple set to Closest Vector Enumeration input parameters" begin
    Q3 = matrix(QQ,6,6,[1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1]);
    L3 = matrix(fmpq[1,0,1,0,2,0]);
    c3 = fmpq(4//3);
    cvp = Hecke.convert_type(Q3,L3,c3);
    @test  det(basis_matrix(cvp[1])) == 1
    @test cvp[2] == matrix(fmpq[-1, 0, -1, 0, -2, 0])
    @test cvp[3] == 14//3
  end

  @testset "Conversion from Closest Vector Enumeration to Quadratic Triple set input parameters" begin
    Lat = Zlattice(gram = matrix(QQ,3,3,[1,0,0,0,1,0,0,0,1]));
    vec = matrix(fmpq[-1,-1,-1]);
    upbound = fmpq(2);
    qt = Hecke.convert_type(Lat, vec, upbound);
    @test qt[1] == matrix(QQ,3,3,[1 0 0; 0 1 0; 0 0 1])
    @test qt[2] == matrix(fmpq[1,1,1])
    @test qt[3] == 1
  end

  @testset "Quadratic Triple Tests" begin
    @testset "Projective Quadratic Triple" begin
      Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
      L = matrix(fmpq[1,1,1,1]);
      c = fmpq(3);
      @test Hecke.projective_quadratic_triple(Q, L, c, 1) == (Q[1, 1], L[1], 0)
      @test Hecke.projective_quadratic_triple(Q, L, c, 4) == (Q, L, c)
    end

    @testset "Finite set of integers for a 1-dimensional quadratic triple" begin
      Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
      L = matrix(fmpq[1,1,1,1]);
      c = fmpq(3);
      Q1 = Hecke.projective_quadratic_triple(Q, L, c, 1)[1];
      L1 = Hecke.projective_quadratic_triple(Q, L, c, 1)[2];
      c1 = Hecke.projective_quadratic_triple(Q, L, c, 1)[3];
      @test Hecke.range_ellipsoid_dim1(Q1, L1, c1) == -2:0
    end

    @testset "Different Types of Extended Quadratic Triple" begin
      @testset "Quadratic triple extension from a 1-dimensional quadratic triple" begin
        Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
        L = matrix(fmpq[1,1,1,1]);
        c = fmpq(3);
        @test Hecke.positive_quadratic_triple(fmpz(1), Q, L, c)[1] == Q[2:4, 2:4]
        @test Hecke.positive_quadratic_triple(fmpz(1), Q, L, c)[2] == L[2:4, 1:1]
        @test Hecke.positive_quadratic_triple(fmpz(1), Q, L, c)[3] == 6
      end

      @testset "Quadratic triple extension from a m-dimensional quadratic triple to m+n-1 quadratic triple" begin
        Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
        L = matrix(fmpq[1,1,1,1]);
        c = fmpq(3);
        @test Hecke.positive_quadratic_triple2(fmpz[-2,0,1], Q,L,c)[1][1] == 1
        @test Hecke.positive_quadratic_triple2(fmpz[-2,0,1], Q,L,c)[2][1] == 1
        @test Hecke.positive_quadratic_triple2(fmpz[-2,0,1], Q,L,c)[3] == 6
      end
    end
  end

  @testset "Finite set tests" begin
    Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
    Q1 = -Q;
    L = matrix(fmpq[1,1,1,1]);
    c = fmpq(3);
    @test Hecke.closest_vectors(Q, L, c, sorting=true)[1] == [-2, -1, -1, -1]
    @test size(Hecke.closest_vectors(Q, L, c), 1) == 9
    @test Hecke.closest_vectors(Q, L, c, sorting=true)[1] == Hecke.closest_vectors(Q1,L,c, sorting=true)[1]

    @test Hecke.closest_vectors(Q, L, c, equal=true, sorting=true)[1] == [-2, -1, -1, -1]
    @test size(Hecke.closest_vectors(Q, L, c, equal=true), 1) == 8
    @test Hecke.closest_vectors(Q, L, c, equal=true, sorting=true)[1] == Hecke.closest_vectors(Q1,L,c, sorting=true)[1]

    L1, L2, L3 = Hecke.convert_type(Q,L,c)
    @test Hecke.closest_vectors(L1, L2, L3, sorting=true)[1] == Hecke.closest_vectors(Q, L, c, sorting=true)[1]
    @test size(Hecke.closest_vectors(L1, L2, L3), 1) == size(Hecke.closest_vectors(Q, L, c), 1)
    @test Hecke.closest_vectors(L1, L2, L3, equal=true, sorting=true)[1] == Hecke.closest_vectors(Q, L, c, equal=true, sorting=true)[1]
    @test size(Hecke.closest_vectors(L1, L2, L3, equal=true), 1) == size(Hecke.closest_vectors(Q, L, c, equal=true), 1)
    t = Hecke.closest_vectors(L1, L2, L3, equal=true)
    L21 = Vector{typeof(L2[1])}()
    for i1 in L2
      push!(L21, i1)
    end
    for i in 1:size(t)[1]
      @test inner_product(ambient_space(L1), L21-t[i], L21-t[i]) == L3
    end

    x = matrix(QQ, 4, 1, Hecke.closest_vectors(Q, L, c)[5]);
    xt = transpose(matrix(QQ, 4, 1, Hecke.closest_vectors(Q, L, c)[5]));
    R = (xt * Q * x)[1] + (2 * xt * L)[1] + c;
    @test  R <= 0

    y = matrix(QQ, 4, 1, Hecke.closest_vectors(Q1, L, c)[5]);
    yt = transpose(matrix(QQ, 4, 1, Hecke.closest_vectors(Q1, L, c)[5]));
    S = (yt * Q * y)[1] + (2 * yt * L)[1] + c;
    @test  S <= 0
  end
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demogeom()

Demonstration of how to construct a 2D geometric shape from a piece-wise set of
NURBSs.

# Examples:
```julia
julia> demogeom()
```
"""
function demogeom()
coefs = [0.0 7.5 15.0 25.0 35.0 30.0 27.5 30.0;
         0.0 2.5  0.0 -5.0  5.0 15.0 22.5 30.0];
knots = vec([0.0 0.0 0.0 1/6 1/3 1/2 2/3 5/6 1.0 1.0 1.0]);

# Geometric boundary curve
geom = [
  nrbmak(coefs,knots);
  nrbline([30.0;30.0],[20.0;30.0]);
  nrbline([20.0;30.0],[20.0;20.0]);
  nrbcirc(10.0,[10.0;20.0],1.5*pi,0.0);
  nrbline([10.0;10.0],[0.0;10.0]);
  nrbline([0.0;10.0],[0.0;0.0]);
  nrbcirc(5.0,[22.5;7.5])
];

plot(title = "2D Geometry formed by a series of NURBS curves.",
     framestyle=:box,
     legend = false);
for i in geom
  nrbplot!(i,50);
end # for i
current();
end # demogeom


text:
using DeconvOptim
using Test
using FFTW, Noise, Statistics, Zygote
using Random
using TestImages, Noise
using Pkg
Pkg.add(url="https://github.com/emmt/OptimPackNextGen.jl")
using OptimPackNextGen

 # fix seed for reproducibility
Random.seed!(42)

@testset "Utils" begin
    @test gpu_or_cpu(randn((2,2))) == Array
    include("utils.jl")
end

include("analysis_tools.jl")
include("hessian_schatten_norm.jl")
include("conv.jl")
include("mappings.jl")

include("forward_models.jl")

include("lossfunctions.jl")

 # testing is rather hard, but include at least some basic testing
include("regularizer.jl")

include("main.jl")


text:
function palindrome(s)
    len = length(s)
    if(len==0 || len==1)
        return true
    end
    if(s[1] == s[len])
        return palindrome(SubString(s,2,len-1))
    end
    return false
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractCase end
export AbstractCase, EnumCase, SwitchCase,
       Leaf, Shaped

Continuation = Symbol
Branch = Pair{PatternInfo, Tuple{LineNumberNode, Continuation}}

"""for generating patterns with one-by-one checks
"""
struct EnumCase <: AbstractCase
    cases :: Vector{AbstractCase}
end

"""for generating patterns with a optimized switch
"""
struct SwitchCase <: AbstractCase
    cases :: Dict{TypeObject, AbstractCase}
end

"""specifying the code body of a case
"""
struct Leaf <: AbstractCase
    cont :: Continuation
end

"""checking the shape of the datum with the defined patterns
"""
struct Shaped <: AbstractCase
    pattern :: PatternInfo
    ln :: LineNumberNode
    case :: AbstractCase
end

@nospecialize

function case_split!(result::Vector{Branch}, branches :: Vector{Branch})
    for (p, (ln, branch)) in branches
        if p.pattern isa Or
            case_split!(result, Branch[info => (ln, branch) for info in p.pattern.ps])
        else
            push!(result, p => (ln, branch))
        end
    end
end

function build_dyn(top::TypeObject, branches::Vector{Branch})::AbstractCase
    length(branches) === 1 && return begin
        br = branches[1]
        Shaped(br.first, br.second[1], Leaf(br.second[2]))
    end
    @assert !isempty(branches)
    enum_cases = AbstractCase[]
    labels = Vector{TypeObject}(undef, length(branches))
    groups = [1]
    for i in eachindex(branches)
        current_type = branches[i].first.typetag
        @assert current_type <: top
        if current_type === top
            # Detected forcely specified matching order.
            # e.g:
            #    match val with
            #    | (_: Int) -> ...
            #    | _ -> ...
            #    | (_ : String) -> ..
            # Although Int and String are orthogonal,
            # the second pattern `_` breaks their merging
           push!(groups, i)
           push!(groups, i + 1)
           labels[i] = top
           continue
        end

        last_group_start = groups[end]
        non_orthogonal_types = TypeObject[]
        non_orthogonal_indices = Int[]
        for j in last_group_start : i - 1
            if typeintersect(current_type, labels[j]) !== Base.Bottom
            # Merging branches
            # e.g:
            #    match val with
            #   | (_ : Int) ->
            #   | (_ : String) ->
            #   | (_ : Float64) ->
            # We gives the 3 patterns such a label array:
            #   [Real, String, Real]
            # where, the 1st and 3rd branches merged
                push!(non_orthogonal_indices, j)
                push!(non_orthogonal_types, labels[j])
            end
        end
        merged_type = reduce(typejoin, non_orthogonal_types, init=current_type)
        for j in non_orthogonal_indices
            labels[j] = merged_type
        end
        labels[i] = merged_type
    end

    push!(groups, length(branches) + 1)
    n_groups = length(groups)
    for i_group in 1:n_groups-1
        start = groups[i_group]
        final = groups[i_group+1] - 1

        if start === final
            br = branches[start]
            push!(enum_cases, Shaped(br.first, br.second[1], Leaf(br.second[2])))
            continue
        elseif start > final
            continue
        end

        switch_map = Dict{TypeObject, Vector{Branch}}()
        for i in start:final
            vec = get!(switch_map, labels[i]) do
                Branch[]
            end
            push!(vec, branches[i])
        end
        switch = Pair{TypeObject, AbstractCase}[
            top => build_dyn(top, branches) for (top, branches) in switch_map
        ]
        push!(
            enum_cases,
            SwitchCase(Dict(switch))
        )
    end

    if length(enum_cases) === 1
        enum_cases[1]
    else
        EnumCase(enum_cases)
    end
end

function case_merge(branches::Vector{Branch})
    top = reduce(typejoin, TypeObject[pattern.typetag for (pattern, _) in branches])
    build_dyn(top, branches)
end

@specialize

max_stars_repo_name: KristofferC/AbstractPattern

max_issues_repo_name: KristofferC/AbstractPattern

max_forks_repo_name: KristofferC/AbstractPattern


text:
##### Difference operators
@inline δx⁺(i, j, k, c) = @inbounds c[i+1, j, k] - c[i, j, k]
@inline δy⁺(i, j, k, c) = @inbounds c[i, j+1, k] - c[i, j, k]
@inline δz⁺(i, j, k, c) = @inbounds c[i, j, k] - c[i, j, k+1]

@inline δx⁰(i, j, k, c) = @inbounds c[i, j, k] - c[i-1, j, k]
@inline δy⁰(i, j, k, c) = @inbounds c[i, j, k] - c[i, j-1, k]
@inline δz⁰(i, j, k, c) = @inbounds c[i, j, k-1] - c[i, j, k]

##### used only when halo points ≥ 2
@inline δx⁻(i, j, k, c) = @inbounds c[i-1, j, k] - c[i-2, j, k]
@inline δy⁻(i, j, k, c) = @inbounds c[i, j-1, k] - c[i, j-2, k]
@inline δz⁻(i, j, k, c) = @inbounds c[i, j, k-2] - c[i, j, k-1]
#####

##### Difference operators with functions
@inline δx⁺(i, j, k, grid, f::F, args...) where F<:Function = f(i+1, j, k, grid, args...) - f(i, j, k, grid, args...)
@inline δy⁺(i, j, k, grid, f::F, args...) where F<:Function = f(i, j+1, k, grid, args...) - f(i, j, k, grid, args...)
@inline δz⁺(i, j, k, grid, f::F, args...) where F<:Function = f(i, j, k, grid, args...) - f(i, j, k+1, grid, args...)

##### Derivative operators
@inline ∂x⁰(i, j, k, grid, c) = δx⁰(i, j, k, c) / ΔxF(i, j, k, grid)
@inline ∂y⁰(i, j, k, grid, c) = δy⁰(i, j, k, c) / ΔyF(i, j, k, grid)
@inline ∂z⁰(i, j, k, grid, c) = δz⁰(i, j, k, c) / ΔzF(i, j, k, grid)

##### Diffusive fluxes
@inline diffusive_flux_x(i, j, k, grid, κˣ::Number, c) = κˣ * Ax(i, j, k, grid) * ∂x⁰(i, j, k, grid, c)
@inline diffusive_flux_y(i, j, k, grid, κʸ::Number, c) = κʸ * Ay(i, j, k, grid) * ∂y⁰(i, j, k, grid, c)
@inline diffusive_flux_z(i, j, k, grid, κᶻ::Number, c) = κᶻ * Az(i, j, k, grid) * ∂z⁰(i, j, k, grid, c)

##### Laplacian diffusion operator
@inline function κ∇²(i, j, k, grid, κˣ, κʸ, κᶻ, c)
    return 1/volume(i, j, k, grid) * (δx⁺(i, j, k, grid, diffusive_flux_x, κˣ, c) +
                       δy⁺(i, j, k, grid, diffusive_flux_y, κʸ, c) +
                       δz⁺(i, j, k, grid, diffusive_flux_z, κᶻ, c))
end

##### calculate the tendency by diffusion for tracer c
@kernel function calc_diffusion!(Gc, grid, κˣ, κʸ, κᶻ, c, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + grid.Hx
    jj = j + grid.Hy
    kk = k + grid.Hz
    @inbounds Gc[ii, jj, kk] = Gc[ii, jj, kk] + κ∇²(ii, jj, kk, grid, κˣ, κʸ, κᶻ, c) * ΔT
end

function nut_diffusion!(Gcs, arch::Architecture, g, nutrients, κˣ, κʸ, κᶻ, ΔT)
    calc_diffusion_kernel! = calc_diffusion!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))

    events=[]
    for name in nut_names
        event = calc_diffusion_kernel!(Gcs[name].data, g, κˣ, κʸ, κᶻ, nutrients[name].data, ΔT, dependencies=barrier)
        push!(events, event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end
nut_diffusion!(Gcs, arch::Architecture, g, nutrients, κ, ΔT) =
    nut_diffusion!(Gcs, arch::Architecture, g, nutrients, κ, κ, κ, ΔT)


text:
using Base.Test

# This is a stand-alone test, showing how to call a function when setting obj.field = something.
# It would be useful for calling getters and setters on a wrapped C++ class
# The example here uses CppClass to stand-in for a real C++ class and allow a pure Julia demo

# type to wrap any field of a wrapped class, taking pointers in the type parameters for methods and data
immutable Property{ObjectT,ValueT,SetPtr,GetPtr,ObjPtr}
    set::Ptr{Void}
    get::Ptr{Void}
    obj::Ptr{Void}
    Property() = new(SetPtr,GetPtr,ObjPtr)
end

# Conversion from a field value type to a property uses the pointers in the parameters to call the setter
function Base.convert{ObjectT,ValueT,SetPtr,GetPtr,ObjPtr}(::Type{Property{ObjectT,ValueT,SetPtr,GetPtr,ObjPtr}}, v::ValueT)
    set_func = unsafe_pointer_to_objref(SetPtr)::Function
    cpp_obj = unsafe_pointer_to_objref(ObjPtr)::ObjectT
    set_func(cpp_obj, v)
    Property{ObjectT,ValueT,SetPtr,GetPtr,ObjPtr}()
end

# This is a place-holder for an actual C++ class
type CppClass
    a::Int
end

# getter and setter methods, which would normally call wrapped C++ methods
set_a(obj::CppClass, val::Int) = @show obj.a = val
get_a(obj::CppClass) = obj.a

# A wrapper for e.g. a C++ class, stored as a pointer. The pointer itself is a parameter
type WrappedCpp{CppPtr}
    cpp::Ptr{Void}
    a::Property{CppClass, Int, pointer_from_objref(set_a), pointer_from_objref(get_a), CppPtr}
    WrappedCpp() = new(CppPtr, Property{CppClass, Int, pointer_from_objref(set_a), pointer_from_objref(get_a), CppPtr}())
end

# Create a "C++" object
cppobj = CppClass(1);
# Wrap it
wrapped = WrappedCpp{pointer_from_objref(cppobj)}()

# Set a field on the wrapper
wrapped.a = 6

# Test that the value set on the wrapper propagated to the wrapped object
@test cppobj.a == 6


text:
# Plotting
function plotEKF(ekf, xtrue, n)
    ts   = ekf.txp[1:n]
    xhat = ekf.xhats[1:n, :]
    xt   = xtrue[1:n, :]
    es   = ekf.es[1:n, :]
    σ311 = 3*sqrt.(ekf.Ps[1:n,1])
    σ322 = 3*sqrt.(ekf.Ps[1:n,2])
    σ333 = 3*sqrt.(ekf.Ps[1:n,3])
    σ344 = 3*sqrt.(ekf.Ps[1:n,4])
    σ355 = 3*sqrt.(ekf.Ps[1:n,5])
    σ366 = 3*sqrt.(ekf.Ps[1:n,6])
    σ377 = 3*sqrt.(ekf.Ps[1:n,7])

    mat"""
    ts      = $ts;
    ts      = (ts - $(ekf.scSim.t0))/86400;
    xhat    = $xhat;
    xt      = $xt;
    es      = $es;
    s311    = $σ311;
    s322    = $σ322;
    s333    = $σ333;
    s344    = $σ344;
    s355    = $σ355;
    s366    = $σ366;
    s377    = $σ377;

    figure()
    subplot(3,2,1)
    plot(ts, es(:,1), 'k')
    hold on
    plot(ts, s311, 'r')
    plot(ts, -s311, 'r')
    xlim([0, 33.2])
    ylim([-25, 25])
    grid on
    ylabel("\$r_x\$ Error, km", "Interpreter", "latex")

    subplot(3,2,3)
    plot(ts, es(:,1), 'k')
    hold on
    plot(ts, s322, 'r')
    plot(ts, -s322, 'r')
    xlim([0, 33.2])
    ylim([-100, 100])
    grid on
    ylabel("\$r_y\$ Error, km", "Interpreter", "latex")

    subplot(3,2,5)
    plot(ts, es(:,3), 'k')
    hold on
    plot(ts, s333, 'r')
    plot(ts, -s333, 'r')
    xlim([0, 33.2])
    ylim([-50, 50])
    xlabel("Time, days", "Interpreter", "latex")
    ylabel("\$r_z\$ Error, km", "Interpreter", "latex")
    grid on

    subplot(3,2,2)
    plot(ts, es(:,4)*1e3, 'k')
    hold on
    plot(ts, s344*1e3, 'r')
    plot(ts, -s344*1e3, 'r')
    xlim([0, 33.2])
    ylim([-6, 6])
    grid on
    ylabel("\$v_x\$ Error, m/s", "Interpreter", "latex")
    leg1 = legend('Estimation Error', 'EKF \$3\\sigma\$', 'Interpreter', 'latex');

    subplot(3,2,4)
    plot(ts, es(:,5)*1e3, 'k')
    hold on
    plot(ts, s355*1e3, 'r')
    plot(ts, -s355*1e3, 'r')
    xlim([0, 33.2])
    ylim([-10, 10])
    grid on
    ylabel("\$v_y\$ Error, m/s", "Interpreter", "latex")

    subplot(3,2,6)
    plot(ts, es(:,6)*1e3, 'k')
    hold on
    plot(ts, s366*1e3, 'r')
    plot(ts, -s366*1e3, 'r')
    xlim([0, 33.2])
    ylim([-6, 6])
    grid on
    xlabel("Time, days", "Interpreter", "latex")
    ylabel("\$v_z\$ Error, m/s", "Interpreter", "latex")

    set(gca, "fontname", "Times New Roman", "fontsize", 10)
    set(leg1, "fontname", "Times New Roman", "fontsize", 8)
    set(gcf, "PaperUnits", "inches", "PaperPosition", [0.25, 0.25, 8.0, 8.0])
    set(gcf, "PaperPositionMode", "Manual")
    print("./figures/EKFPosVelError.eps", "-depsc", "-r300")

    figure()
    plot(ts, es(:,7), 'k')
    hold on
    plot(ts, s377, 'r')
    plot(ts, -s377, 'r')
    xlim([0, 33.2])
    xlabel("Time, days", "Interpreter", "latex")
    ylabel("m Error, kg", "Interpreter", "latex")
    leg2 = legend('Estimation Error', 'EKF \$3\\sigma\$', 'Interpreter', 'latex');

    set(gca, "fontname", "Times New Roman", "fontsize", 10)
    set(leg2, "fontname", "Times New Roman", "fontsize", 8)
    set(gcf, "PaperUnits", "inches", "PaperPosition", [0.25, 0.25, 4.0, 3.0])
    set(gcf, "PaperPositionMode", "Manual")
    print("./figures/EKFMass.eps", "-depsc", "-r300")

    figure()
    plot3(xhat(:,1),xhat(:,2), xhat(:,3), 'r')
    hold on
    plot3(xt(:,1),xt(:,2),xt(:,3), 'b')
    grid on
    axis equal
    """
end

plotEKF(ekf, xtrue, length(ekf.txp))

# Create GPS plot
gpsTempSim = GPSSim(2170, 0, 2170, 1)
plotGPS(gpsTempSim)

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
## ExoplanetsSysSim/examples/dr25_gaia_fgk/make_q1_q17_christiansen.jl
## (c) 2018 Danley C. Hsu & Eric B. Ford
# Script for creating the Q1-Q16 FGK subset stellar catalog in jld format
#    used in the Hsu et al. 2018 paper

using ExoplanetsSysSim
using JLD

include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","poisson_tests", "christiansen_func.jl"))

stellar_catalog_file_in = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_stellar.csv")
stellar_catalog_file_out = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1q17_dr25_kepler_fgk.jld")

stellar_catalog = setup_star_table_christiansen(ascii(stellar_catalog_file_in))

save(stellar_catalog_file_out,"stellar_catalog", stellar_catalog)


text:
using StatsMonitor
using ProgressMeter
using Dates

s = Stats(Second(1))

@showprogress for i in 1:2000
    sleep(.01)
    fit!(s, :id1 => randn() + i/1000)
    fit!(s, :id2 => i * randn())
    fit!(s, :id3 => rand(["A","A","B","C"]))

    # different data types, same id
    fit!(s, :id4 => rand() < .8)
    rand() < i/2000 && fit!(s, :id4 => nothing)
end

plot(s, size=(1000, 500))


text:
using AbstractTrees

struct Tree{T}
  value::T
  children::Vector{Tree{T}}
end

Tree{T}(x::T, xs::Tree{T}...) where T = Tree{T}(x, [xs...])
Tree{T}(x) where T = Tree(convert(T, x))

Tree(x::T, xs::Tree{T}...) where T = Tree{T}(x, xs...)

AbstractTrees.children(t::Tree) = t.children
AbstractTrees.printnode(io::IO, t::Tree) = show(io, t.value)

Base.show(io::IO, t::Type{Tree}) = print(io, "Tree")
Base.show(io::IO, t::Type{Tree{T}}) where T = print(io, "Tree{", @isdefined(T) ? T : :T, "}")

function Base.show(io::IO, t::Tree)
  println(io, typeof(t))
  print_tree(io, t)
end

Base.getindex(t::Tree, i::Integer) = t.children[i]
Base.getindex(t::Tree, i::Integer, is::Integer...) = t[i][is...]

# Utilities

isleaf(t) = isempty(children(t))

leaves(xs::Tree) = map(x -> x.value, Leaves(xs))

Base.map(f, t::Tree, ts::Tree...) =
  Tree{Any}(f(map(t -> t.value, (t, ts...))...),
            [map(f, chs...) for chs in zip(map(t -> t.children, (t, ts...))...)]...)


text:
#= ========== LevelDB ==================
LevelDb specific code.

Each DB should export the symbols in KVS.jl
=#
#using LevelDB

global path_splice = ""

mutable struct KVSLevelDb <: AbstractKVS
    name::String     # Name of db
    handle           # Handle to db object
end
KVSLevelDb(name::String) = KVSLevelDb(name, levelDb_open(name))

function levelDb_open(name::String)
    db_root::String = getindex(ENV, "HOME")
    db_dir::String = db_root * "/kinant/" * path_splice * "/"
    if stat(db_dir).inode == 0
        Base.mkdir(db_dir)
    end
    handle = LevelDB.open_db(db_dir * name, true)
    handle
end

function kvs_path_append(snippet::String)
    global path_splice = snippet
end

function kvs_open(hdl::KVSLevelDb)
    hdl.handle = levelDb_open(hdl.name)
    hld
end

function kvs_close(hdl::KVSLevelDb)
    close_db(hdl.handle)
    hdl.name = ""
    hdl.handle = nothing
end

function kvs_put(hdl::KVSLevelDb, key, value)
    b = byte_array(value)
    LevelDB.db_put(hdl.handle, byte_array(key), b, sizeof(b))
end

function kvs_get(hdl::KVSLevelDb, key)
    r = LevelDB.db_get(hdl.handle, byte_array(key))
    return length(r) > 0 ? array_to_type(r) : nothing
end

"""
    kvs_get_many(db::KVSLevelDb, first, last, max::Int; ktype::Type=Any; vtype::Type=Any)
Return *max* entries beginning from *first* and upto but not including *last*.
Optionally specify the type of key and value.
"""
function kvs_get_many(db::KVSLevelDb, first, last, max::Integer; ktype::Type=Any, vtype::Type=Any)
    #first = "\u00000000000000000000000000000000"
    #last  = "\uffffffffffffffffffffffffffffffff"
    keys = Vector{ktype}(max)
    vals = Vector{vtype}(max)
    i::Int = 0
    for (k, v) in LevelDB.db_range(db.handle, byte_array(first), byte_array(last))
        if array_to_type(k) == last break end # Don't include last
        i += 1
        if i > max # Fill only max entries
            i -= 1
            break
        end
        keys[i] = array_to_type(k)
        vals[i] = array_to_type(v)
    end
    (keys, vals, i)
end

#=
function kvs_list_channels(db::KVSLevelDb)
    @assert db.name == "properties"
    first = "\u00000000000000000000000000000000"
    last  = "\uffffffffffffffffffffffffffffffff"
    keys = Vector{id_t}(0)
    for (k, v) in LevelDB.db_range(db.handle, first, last)
        key = array_to_type(k, id_t)
        push!(keys, key)
    end
    keys
end
=#

"""
Update the stream in the db atomically. For level Db this is just a
plain update as each process holds a lock on the db.
"""
function kvs_atomic_update(db::KVSLevelDb, cid, s, op::atomic_op)
    @assert db.name == "recent"
    if op == ATOMIC_INC
        # if s.state != S_OPEN throw(StreamStateException("State != S_OPEN")) end
        s.ref = s.ref + 1
        s.seq = s.seq + 1
    elseif op == ATOMIC_DEC
        s.ref = s.ref - 1
    end
    # TODO: for consul, need to check if someone else is commiting (S_COMMITTING)
    # the stream. Then this process should back off. The best way to do that
    # would be to have a state transition checker throw an exception if the
    # transition is not right.

    val = byte_array(s)
    LevelDB.db_put(db.handle, byte_array(cid), val, sizeof(val))
    s
end

"""
Try to atomically update db. For level Db this is just a
plain update as each process holds a lock on the db.
"""
function kvs_atomic_try_update(db::KVSLevelDb, cid, s)
    @assert db.name == "recent"
    val = byte_array(s)
    LevelDB.db_put(db.handle, byte_array(cid), val, sizeof(val))
    true
end

"""
Warning! Extreme caution! Deletes db!
"""
function kvs_clobber!(db::KVSLevelDb)
    db_root::String = getindex(ENV, "HOME")
    db_path = db_root * "/thimble/" * db.name
    kvs_close(db)
    run(`mkdir -p $db_root/bak`)
    run(`mv $db_path $db_root/bak/`)
    run(`unlink $db_path`)
end

"""
    db_chunk_exists(fp::fp_t)

Returns true if chunk exists in db.
"""
function chunk_exists(db::KVSLevelDb, fp)
    c = LevelDB.db_get(db.handle, byte_array(fp))
    if (length(c) == 0) return false else return true end
end

"""
return true if chunk already exists, else put chunk and return false
"""
function kvs_put_chunk(db::KVSLevelDb, input::Vector{UInt8}, fp, first::Int64, last::Int64)
    if chunk_exists(db, fp)
        return true
    else
        val = byte_array(input[first:last])
        LevelDB.db_put(db.handle, byte_array(fp), val, length(val))
    end
end

"""
Atomically increment the chunk ref counts for all chunks in array.
The incrementing of each chunk's ref is atomic.

For levelDb increment is implemented by a get and put.
"""
function kvs_atomic_batch_increment(db::KVSLevelDb, chunks)
    @assert db.name == "chunk_ref"
    saved::UInt64 = 0        # Bytes saved for this extent
    prev_offset::UInt64 = 0
    index::UInt32 = 0
    exists::BitVector =  falses(length(chunks))
    for i in chunks
        ref::UInt64 = 0
        index += 1
        ret = LevelDB.db_get(db.handle, byte_array(i.fp)) # Get ref
        if (length(ret) != 0)
            ref = array_to_type(ret, UInt64)
        end

        # Compute return values
        if (ref > 0)
            exists[index] = true
            saved += i.offset - prev_offset
        end
        prev_offset = i.offset

        # Increment ref and save
        ref = ref + 1
        val = byte_array(ref)
        LevelDB.db_put(db.handle, byte_array(i.fp), val, sizeof(val)) # Put ref
    end
    (exists, saved)
end

# helper functions


text:
@testset "softmax" begin
    for dims in [(5,5), (5,)]
        x = randn(Float64, dims)
        y = softmax(x)
        dy = randn(Float64, dims)
        gputest(softmax, x)
        gputest(∇softmax, dy, x, y, checkgrad=false)
        y = logsoftmax(x)
        gputest(logsoftmax, x)
        gputest(∇logsoftmax, dy, x, y, checkgrad=false)
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/OpenBLAS_jll.jl
baremodule OpenBLAS_jll
using Base, Libdl, CompilerSupportLibraries_jll, Base.BinaryPlatforms
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libopenblas

# These get calculated in __init__()
const PATH = Ref("")
const LIBPATH = Ref("")
artifact_dir = ""
libopenblas_handle = C_NULL
libopenblas_path = ""

if Base.USE_BLAS64
    const libsuffix = "64_"
else
    const libsuffix = ""
end

if Sys.iswindows()
    const libopenblas = "libopenblas$(libsuffix).dll"
elseif Sys.isapple()
    const libopenblas = "@rpath/libopenblas$(libsuffix).dylib"
else
    const libopenblas = "libopenblas$(libsuffix).so"
end

function __init__()
    # make sure OpenBLAS does not set CPU affinity (#1070, #9639)
    if !haskey(ENV, "OPENBLAS_MAIN_FREE")
        ENV["OPENBLAS_MAIN_FREE"] = "1"
    end

    global libopenblas_handle = dlopen(libopenblas)
    global libopenblas_path = dlpath(libopenblas_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    LIBPATH[] = dirname(libopenblas_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libopenblas_path() = libopenblas_path

end  # module OpenBLAS_jll

max_stars_repo_name: BSnelling/julia

max_issues_repo_name: BSnelling/julia

max_forks_repo_name: BSnelling/julia


text:
"""
Vector of a judiVector and a judiWeight
"""
mutable struct judiVStack{vDT<:Number}
    m::Integer
    n::Integer
    components::Array{Any, 1}
end

function vcat(x::judiMultiSourceVector, y::judiMultiSourceVector)
    components = Array{Any}(undef, 2)
    components[1] = x
    components[2] = y
    m = length(x)+length(y)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::judiVStack, y::judiMultiSourceVector)
    components = Array{Any}(undef, length(x.components) + 1)
    for i=1:length(x.components)
        components[i] = x.components[i]
    end
    components[end] = y
    m = x.m+length(y)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::judiMultiSourceVector, y::judiVStack)
    components = Array{Any}(undef, length(y.components) + 1)
    components[1] = x
    for i=2:length(y.components)+1
        components[i] = y.components[i-1]
    end
    m = y.m + length(x)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::judiVStack, y::judiVStack)
    nx = length(x.components)
    ny = length(y.components)
    components = Array{Any}(undef, nx+ny)
    for i=1:nx
        components[i] = x.components[i]
    end
    for i=nx+1:nx+ny
        components[i] = y.components[i-nx]
    end
    m = x.m + y.m
    n = 1
    return judiVStack{Float32}(m, n, components)
end

for T ∈ [judiVector, judiWeights, judiWavefield]
    @eval begin
        function vcat(x::$T, y::$T)
            xn = deepcopy(x)
            push!(xn, y)
            return xn
        end
    end
end

function *(F::joAbstractLinearOperator, v::judiVStack)
    return sum(F.fop[i]*v[i] for i=1:length(v.components))
end

(msv::judiMultiSourceVector{T})(x::judiVStack{T}) where {T} = x
############################################################
## overloaded Base functions

# conj(jo)
conj(a::judiVStack{vDT}) where vDT =
    judiVStack{vDT}(a.m,a.n,a.components)

# transpose(jo)
transpose(a::judiVStack{vDT}) where vDT =
    judiVStack{vDT}(a.n,a.m,a.components)

# adjoint(jo)
adjoint(a::judiVStack{vDT}) where vDT =
    judiVStack{vDT}(a.n,a.m,a.components)

##########################################################
# Utilities

size(x::judiVStack) = (x.m, x.n)
size(x::judiVStack, ind::Integer) = (x.m, x.n)[ind]

length(x::judiVStack) = x.m

eltype(v::judiVStack{vDT}) where {vDT} = vDT

similar(x::judiVStack) = judiVStack{Float32}(x.m, x.n, 0f0 .* x.components)

similar(x::judiVStack,  element_type::DataType, dims::Union{AbstractUnitRange, Integer}...) = similar(x)

getindex(x::judiVStack, a) = x.components[a]

firstindex(x::judiVStack) = 1

lastindex(x::judiVStack) = length(x.components)

dot(x::judiVStack, y::judiVStack) = sum(dot(x[i],y[i]) for i=1:length(x.components))

function norm(x::judiVStack, order::Real=2)
    if order == Inf
        return max([norm(x[i], Inf) for i=1:length(x.components)]...)
    end
    sum(norm(x[i], order)^order for i=1:length(x.components))^(1/order)
end

iterate(S::judiVStack, state::Integer=1) = state > length(S.components) ? nothing : (S.components[state], state+1)

isfinite(S::judiVStack) = all(isfinite(c) for c in S)

##########################################################

# minus
-(a::judiVStack) = -1*a
+(a::judiVStack, b::judiVStack) = judiVStack{Float32}(a.m, a.n, a.components + b.components)
-(a::judiVStack, b::judiVStack) = judiVStack{Float32}(a.m, a.n, a.components - b.components)
+(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n,a.components .+ b)
-(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n,a.components .- b)
-(a::Number, b::judiVStack) = judiVStack{Float32}(b.m, b.n, a .- b.components)

*(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n, b .* a.components)
*(a::Number, b::judiVStack) = b * a
+(a::Number, b::judiVStack) = b + a

/(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n, a.components ./ b)

##########################################################

BroadcastStyle(::Type{judiVStack}) = Base.Broadcast.DefaultArrayStyle{1}()

broadcasted(::typeof(+), x::judiVStack, y::judiVStack) = x + y
broadcasted(::typeof(-), x::judiVStack, y::judiVStack) = x - y
broadcasted(::typeof(+), x::judiVStack, y::Number) = x + y
broadcasted(::typeof(-), x::judiVStack, y::Number) = x - y
broadcasted(::typeof(+), y::Number, x::judiVStack) = x + y
broadcasted(::typeof(-), y::Number, x::judiVStack) = x - y

function broadcasted(::typeof(*), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] = x.components[j] .* y.components[j]
    end
    return z
end

function broadcasted!(::typeof(*), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] = x.components[j] .* y.components[j]
    end
    return z
end

function broadcasted(::typeof(/), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] = x.components[j] ./ y.components[j]
    end
    return z
end

function broadcasted(::typeof(*), x::judiVStack, y::Number)
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] .*= y
    end
    return z
end

broadcasted(::typeof(*), y::Number, x::judiVStack) = x .* y

function broadcasted(::typeof(/), x::judiVStack, y::Number)
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] ./= y
    end
    return z
end

function materialize!(x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    for j=1:length(x.components)
        try
            x.components[j].data .= y.components[j].data
        catch e
            x.components[j].weights .= y.components[j].weights
        end
    end
    return x
end

function broadcast!(::typeof(identity), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    copy!(x,y)
end

function broadcasted(::typeof(identity), x::judiVStack)
    return x
end

function copy!(x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    for j=1:length(x.components)
        try
            x.components[j].data .= y.components[j].data
        catch e
            x.components[j].weights .= y.components[j].weights
        end
    end
end

function isapprox(x::judiVStack, y::judiVStack; rtol::AbstractFloat=sqrt(eps()), atol::AbstractFloat=0.0)
    x.m == y.m || throw("Shape error")
    all(isapprox(xx, yy; rtol=rtol, atol=atol) for (xx, yy)=zip(x.components, y.components))
end

############################################################

function A_mul_B!(x::judiMultiSourceVector, F::joCoreBlock, y::judiVStack)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    x.data .= z.data
end

function A_mul_B!(x::judiVStack, F::joCoreBlock, y::judiMultiSourceVector)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.components)
        x.components[j].data .= z.components[j].data
    end
end

mul!(x::judiMultiSourceVector, J::joCoreBlock, y::judiVStack) = A_mul_B!(x, J, y)
mul!(x::judiVStack, J::joCoreBlock, y::judiMultiSourceVector) = A_mul_B!(x, J, y)