text:
using LinearAlgebra

export transform,
       backtransform

"""
    backtransform(Rsets::ReachSolution, options::Options)

Undo a coordinate transformation.

### Input

- `Rsets`  -- flowpipe
- `option` -- problem options containing an `:transformation_matrix` entry

### Output

A new flowpipe where each reach set has been transformed.

### Notes

The transformation is implemented with a lazy `LinearMap`.
"""
function backtransform(Rsets, options::Options)
    transformation_matrix = options[:transformation_matrix]
    if transformation_matrix == nothing
        return Rsets
    end
    return project(Rsets, transformation_matrix)
end

"""
    transform(problem::InitialValueProblem, options::Options)

Interface function that calls the respective transformation function.

### Input

- `problem` -- discrete or continuous initial-value problem
- `option`  -- problem options

### Output

A tuple containing the transformed problem and the transformed options.

### Notes

The functions that are called in the background should return a the transformed
system components `A`, `X0`, and `U`, and also an inverse transformation matrix `M`.
If the system has an invariant, it is transformed as well.
"""
function transform(problem::InitialValueProblem, options::Options)
    method = options[:coordinate_transformation]
    if method == ""
        nothing # no-op
    elseif method == "schur"
        problem, T_inverse = schur_transform(problem)
        options[:transformation_matrix] = T_inverse
    else
        error("the transformation method $method is undefined")
    end
    return (problem, options)
end

"""
    schur_transform(problem::InitialValueProblem)

Applies a Schur transformation to a discrete or continuous initial-value problem.

### Input

- `problem` -- discrete or continuous initial-value problem

### Output

Transformed problem.

### Algorithm

We use Julia's default `schurfact` function to compute a
[Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition)
of the coefficients matrix ``A``.
"""
function schur_transform(problem::InitialValueProblem{PT, ST}
                         ) where {PT <: Union{ConstrainedLinearControlDiscreteSystem, ConstrainedLinearControlContinuousSystem}, ST<:LazySet}

    n = size(problem.s.A, 1)

    # full (dense) matrix is required by schur
    # result S is a struct such that
    # - S.Schur is in Schur form and
    # - A == S.vectors * S.Schur * S.vectors'
    S = schur(Matrix(problem.s.A))

    # recall that for Schur matrices, inv(T) == T'
    Z_inverse = copy(transpose(S.vectors))

    # obtain transformed system matrix
    A_new = S.Schur

    # apply transformation to the initial states
    X0_new = Z_inverse * problem.x0

    # apply transformation to the inputs
    B_new = Z_inverse * problem.s.B
    U_new = problem.s.U

    # matrix for reverting the transformation again
    T_inverse = S.vectors

    # apply transformation to the invariant
    if hasmethod(stateset, Tuple{typeof(problem.s)})
        invariant_new = T_inverse * problem.s.X
    else
        invariant_new = Universe(n)
    end

    system_new = _wrap_system(PT, A_new, B_new, invariant_new, U_new)
    problem_new = InitialValueProblem(system_new, X0_new)
    return problem_new, T_inverse
end

function _wrap_system(PT::Type{<:ConstrainedLinearControlDiscreteSystem},
                      A, B, invariant, U)
    return ConstrainedLinearControlDiscreteSystem(A, B, invariant, U)
end

function _wrap_system(PT::Type{<:ConstrainedLinearControlContinuousSystem},
                      A, B, invariant, U)
    return ConstrainedLinearControlContinuousSystem(A, B, invariant, U)
end


text:
using ReachabilityAnalysis

Dict([
    "opts" => Dict(:T=>5.0, :alg=>GLGM06(δ=0.005)),
    "X0"   => BallInf([1.0, 0.0, 0.0, 0.0, 0.0], 0.1),
    "info" => Dict(
        "name" => "five_dim_sys",
        "dim" => 5,
        "linear" => true,
        "hybrid" => false
    )
])

ext: jl

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
export View, homogenize, project, homography, bounds, camera_resolution

type View
  camera :: M34
  image :: Array{Float32, 3}
end

homogenize(v :: V3) = V4(v[1], v[2], v[3], 1.0)
*(m :: M4, v :: V3) = m * homogenize(v)
*(m :: M34, v :: V3) = m * homogenize(v)
project(v :: V4) = V3(v[1]/v[4], v[2]/v[4], v[3]/v[4])
project(v :: V3) = V2(v[1]/v[3], v[2]/v[3])

"""
Given 2 camera matrices, calculate the homography that takes pixels from camera 1 to
camera 2 through camera 1's z-depth image plane.
"""
homography(c1, c2, z) = M3((Array(c2)/[Array(c1) ; 0 0 0 z])*[eye(3) ; 0 0 1])

function corners(voi)
  (lo,hi) = voi

  [ V3(lo[1], lo[2], lo[3])
    V3(hi[1], lo[2], lo[3])
    V3(lo[1], hi[2], lo[3])
    V3(hi[1], hi[2], lo[3])
    V3(lo[1], lo[2], hi[3])
    V3(hi[1], lo[2], hi[3])
    V3(lo[1], hi[2], hi[3])
    V3(hi[1], hi[2], hi[3])
  ]
end

function bounds(cam, voi)

  function _proj(p)
    px = cam * p
    V3(px[1]/px[3], px[2]/px[3], px[3])
  end

  crnrs = corners(voi)
  xyzs = zip(map(_proj, crnrs)...)
  mns = map(minimum, xyzs)
  mxs = map(maximum, xyzs)

  (V3(mns), V3(mxs))
end

# smallest resolvable distance at unit depth (roughly, 1/f)
camera_resolution(cam) = minimum(svdvals([Array(cam) ; 0 0 0 1]\eye(4,2)))

max_stars_repo_name: twadleigh/aminda

max_issues_repo_name: twadleigh/aminda

max_forks_repo_name: twadleigh/aminda


text:
#!/usr/bin/env julia --project=@.
# module Starter
doc = """Web service starter

Usage:
  $(Base.basename(@__FILE__)) [--port=<num>] [--base_url=<url>] [--bind=<ip>]

Options:
  -h --help         Show this screen
  -p, --port=<num>  Port  [default: 8080]
  -b, --bind=<ip>   Bind address [default: 0.0.0.0]
  --base_url=<url>  Additional URL prefix for the service [default: /]

"""

using DocOpt  # import docopt function
using Pkg

args = docopt(doc, version=Pkg.project().version)
#@info args

using Sockets
BASE_URL = args["--base_url"]
HOST = Sockets.getaddrinfo(args["--bind"])
PORT = parse(Int, args["--port"])

app_server = Pkg.project().name
@info "Activating web service..."
@eval using $(Symbol(app_server))
m = getfield(Main, Symbol(app_server))

# endswith(PROGRAM_FILE, basename(@__FILE__)) && start_server()
m.AppServer.start_server(
    host = string(HOST),
    port = PORT,
    base_url = BASE_URL,
)
# end


text:
volumes = getallvolumes!(client)

@testset "List all Volumes" begin
    @test length(volumes) == 1
    @test volumes[1].id == "506f78a4-e098-11e5-ad9f-000f53306ae1"
    @test volumes[1].region.name == "New York 1"
    @test "metadata" in volumes[1].region.features
    @test "s-3vcpu-1gb" in volumes[1].region.sizes
    @test isempty(volumes[1].droplet_ids)
    @test volumes[1].name == "example"
    @test volumes[1].description == "Block store for examples"
    @test volumes[1].sizegigabytes == 10
    @test volumes[1].created_at == DateTime("2016-03-02T17:00:49")

end;

volume = createvolume!(client; size_gigabytes=10, name="example",
                       description="Block store for examples",
                       region="nyc1")

@testset "Create a new Volume" begin
    @test volume.id == "506f78a4-e098-11e5-ad9f-000f53306ae1"
    @test volume.region.name == "New York 1"
    @test "private_networking" in volume.region.features
    @test isempty(volume.droplet_ids)
    @test volume.name == "example"
    @test volume.description == "Block store for examples"
    @test volume.sizegigabytes == 10
    @test volume.created_at == DateTime("2016-03-02T17:00:49")
end;

volume = getvolume!(client, volume)

@testset "Retrieve an existing Volume" begin
    @test volume.id == "506f78a4-e098-11e5-ad9f-000f53306ae1"
    @test volume.region.name == "New York 1"
    @test "private_networking" in volume.region.features
    @test isempty(volume.droplet_ids)
    @test volume.name == "example"
    @test volume.description == "Block store for examples"
    @test volume.sizegigabytes == 10
    @test volume.created_at == DateTime("2016-03-02T17:00:49")
end;

volume_id = "82a48a18-873f-11e6-96bf-000f53315a41"
snapshots = getallvolumesnapshots!(client, volume_id)

@testset "List Snapshots for a Volume" begin
    @test length(snapshots) == 1
    @test snapshots[1].id == "8eb4d51a-873f-11e6-96bf-000f53315a41"
    @test snapshots[1].name == "big-data-snapshot1475261752"
    @test snapshots[1].regions[1] == "nyc1"
    @test snapshots[1].created_at == DateTime("2016-09-30T18:56:12")
    @test snapshots[1].resourceid == "82a48a18-873f-11e6-96bf-000f53315a41"
    @test snapshots[1].resourcetype == "volume"
    @test snapshots[1].mindisksize == 10
    @test snapshots[1].sizegigabytes == 0
end;

snapshot = snapshotvolume!(client, volume_id;
                                       name="big-data-snapshot1475261774")

@testset "Create a Snapshot from a Volume" begin
    @test snapshot.id == "8fa70202-873f-11e6-8b68-000f533176b1"
    @test snapshot.name == "big-data-snapshot1475261774"
    @test snapshot.regions[1] == "nyc1"
    @test snapshot.created_at == DateTime("2016-09-30T18:56:14")
    @test snapshot.resourceid == "82a48a18-873f-11e6-96bf-000f53315a41"
    @test snapshot.resourcetype == "volume"
    @test snapshot.mindisksize == 10
    @test snapshot.sizegigabytes == 0
end;

@testset "Delete a Volume" begin
    @test deletevolume!(client, volume)
end;

@testset "Delete a volume snapshot" begin
    @test deletevolumesnapshot!(client, snapshot)
end;

action = attachvolume!(client, "7724db7c-e098-11e5-b522-000f53304e51";
                       droplet_id=11612190)

@testset "Attach a Volume to a Droplet" begin
    @test action.id == 72531856
    @test action.status == "completed"
    @test action.actiontype == "attach_volume"
    @test action.started_at == DateTime("2015-11-12T17:51:03")
    @test action.completed_at == DateTime("2015-11-12T17:51:14")
    @test action.resourceid == nothing
    @test action.resourcetype == "volume"
end;

action = removevolume!(client, "7724db7c-e098-11e5-b522-000f53304e51";
                       droplet_id=11612190, region="nyc1")

@testset "Remove a Volume from a Droplet" begin
    @test action.id == 68212773
    @test action.status == "in-progress"
    @test action.actiontype == "detach_volume"
    @test action.started_at == DateTime("2015-10-15T17:46:15")
    @test action.completed_at == nothing
    @test action.resourceid == nothing
    @test action.resourcetype == "backend"
end;

action = resizevolume!(client, "7724db7c-e098-11e5-b522-000f53304e51";
                       size_gigabytes=10, region="nyc1")

@testset "Resize a Volume" begin
    @test action.id == 72531856
    @test action.status == "in-progress"
    @test action.actiontype == "resize"
    @test action.started_at == DateTime("2015-11-12T17:51:03")
    @test action.completed_at == DateTime("2015-11-12T17:51:14")
    @test action.resourceid == nothing
    @test action.resourcetype == "volume"
end;

actions = getallvolumeactions!(client, "7724db7c-e098-11e5-b522-000f53304e51")

@testset "List all Actions for a Volume" begin
    @test length(actions) == 1
    @test actions[1].id == 72531856
    @test actions[1].status == "completed"
    @test actions[1].actiontype == "attach_volume"
    @test actions[1].started_at == DateTime("2015-11-21T21:51:09")
    @test actions[1].completed_at == DateTime("2015-11-21T21:51:09")
    @test actions[1].resourceid == nothing
    @test actions[1].resourcetype == "volume"
end;

action = getvolumeaction!(client, "7724db7c-e098-11e5-b522-000f53304e51",
                           72531856)

@testset "Retrieve an existing Volume Action" begin
    @test action.id == 72531856
    @test action.status == "completed"
    @test action.actiontype == "attach_volume"
    @test action.started_at == DateTime("2015-11-12T17:51:03")
    @test action.completed_at == DateTime("2015-11-12T17:51:14")
    @test action.resourceid == nothing
    @test action.resourcetype == "volume"
end;


text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpScaleType.jl")
TabLevel=TabLevel*"    "

OpScaleType(Op::MajoranaOp) = return typeof(Op.scale)

function OpScaleType(Sp::MajoranaSp)
    ##Here we return the type one would get is all the terms were added together
    typeof(sum([one(OpScaleType(op)) for op in Sp.SpList]))
end

function OpScaleType(Arr::Array)
    AllScales=[one(OpScaleType(op)) for op in Arr]
    typeof(sum(AllScales))
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpScaleType.jl")


text:
using Parameters

#=
Parameters of the algorithm.

Some affect iteration length.
Some affect struct size.

Names from original C implementation in the comments.

Defaults are from original C.
=#

#=
TODO
Constants spread in the code.  Not parameters, here or in the original.

MAX_PASSES

IMAGE_SYNTH_MAX_NEIGHBORS was defined in the original, not here.
It was a hard limit on patchLength parameter

IMAGE_SYNTH_TERMINATE_FRACTION. 0.1

# MAX_PROBES_PER_PASS ??
=#

@with_kw struct ResynthesizerParameters
    maxPatchSize::UInt   = 30   # desired size of patches, some will be shorter
    maxProbeCount::UInt = 200
    withReplacement::Bool = true
    #b::Float64 = -1.1
    #c::UInt8
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
Base.FFTW.set_num_threads(8)
BLAS.set_num_threads(8)
include("../../src/MatlabPort.jl")

using GeometryTypes, MatlabPort, GLAbstraction
import MatlabPort: ISF, Normalize, PressureProject, land, Particles, StaggeredAdvect
import MatlabPort: VelocityOneForm, StaggeredSharp, SchroedingerFlow

function GeometryTypes.isinside(point, lower, upper)
    @inbounds for (p,l,u) in zip(point, lower, upper)
        (p > l && p < u) || return false
    end
    true
end

vol_size = (4,2,2)   # box size
vol_res = (64,32,32) # volume resolution
hbar = 0.1           # Planck constant
dt = 1/48            # time step
tmax = 50            # max time

jet_velocity = [1,0,0]; # jet velocity

nozzle_cen = [2-1.7, 1-0.034, 1+0.066]; # nozzle center
nozzle_len = 0.5;                   # nozzle length
nozzle_rad = 0.5;                   # nozzle radius

n_particles = 10;   # number of particles

isf = MatlabPort.ISF(vol_size, vol_res, hbar, dt)

# Set nozzle
isJet = land(
    abs.(isf.px .- nozzle_cen[1]) .<= nozzle_len./2,
    (isf.py .- nozzle_cen[2]) .^ 2 .+ (isf.pz .- nozzle_cen[3]) .^ 2 .<= nozzle_rad .^ 2
)

# initialize psi
psi1f = ones(size(isf.px))
psi2f = psi1f * 0.01
psi1f, psi2f = Normalize(psi1f, psi2f)

# constrain velocity
kvec = jet_velocity / isf.hbar
omega = sum(jet_velocity .^ 2) ./ (2*isf.hbar);
phase = kvec[1].*isf.px .+ kvec[2].*isf.py .+ kvec[3].*isf.pz;
# convert to complex
psi1 = (1.+0.0im) .* psi1f
psi2 = (1.+0.0im) .* psi2f
for iter = 1:10
    amp1 = abs.(psi1)
    amp2 = abs.(psi2)
    psi1[isJet] = amp1[isJet] .* exp.(1.0im*phase[isJet])
    psi2[isJet] = amp2[isJet] .* exp.(1.0im*phase[isJet])
    psi1, psi2 = PressureProject(isf, psi1, psi2)
end

max_particles = 20_000
max_history = 19

function removethedead!(particles, vol_size)
    filter!(particles.active) do i
        isinside(
            (particles.x[i], particles.y[i], particles.z[i]),
            (0,0,0), vol_size
        )
    end
end

## SET PARTICLES
particle = Particles(
    zeros(Float32, max_particles),
    zeros(Float32, max_particles),
    zeros(Float32, max_particles),
    0, Int[]
)

using GLVisualize; w = glscreen(); @async GLWindow.renderloop(w)

particle_vis = visualize(
    (Circle(Point2f0(0), 0.006f0), (particle.x, particle.y, particle.z)),
    boundingbox = nothing, # don't waste time on bb computation
    billboard = true, indices = particle.active
).children[]
_view(particle_vis, camera = :perspective)

function loop(
        N, isf, psi1, psi2, dt, kvec, omega, isJet, n_particles,
        nozzle_rad, nozzle_cen, particle, particle_vis
    )
    gpu_position_x = particle_vis[:position_x]
    gpu_position_y = particle_vis[:position_y]
    gpu_position_z = particle_vis[:position_z]
    for iter = 1:N
        t = iter * dt
        # incompressible Schroedinger flow
        psi1, psi2 = SchroedingerFlow(isf, psi1, psi2)
        psi1, psi2 = Normalize(psi1, psi2)
        psi1, psi2 = PressureProject(isf, psi1, psi2)

        # constrain velocity
        phase = kvec[1] .* isf.px .+ kvec[2] .* isf.py .+ kvec[3] .* isf.pz .- omega .* t
        amp1 = abs.(psi1)
        amp2 = abs.(psi2)
        psi1[isJet] = amp1[isJet] .* exp.(1.0im*phase[isJet])
        psi2[isJet] = amp2[isJet] .* exp.(1.0im*phase[isJet])
        psi1, psi2  = PressureProject(isf, psi1, psi2)

        # particle birth
        rt   = rand(Float32, n_particles)*2*pi
        newx = nozzle_cen[1] .* ones(Float32, size(rt))
        newy = nozzle_cen[2] .+ 0.9 .* nozzle_rad .* cos.(rt)
        newz = nozzle_cen[3] .+ 0.9 .* nozzle_rad .* sin.(rt)
        append!(particle, newx, newy, newz)
        removethedead!(particle, vol_size)
        # advect and show particles
        vx,vy,vz = VelocityOneForm(isf, psi1, psi2, isf.hbar)
        vx,vy,vz = StaggeredSharp(isf, vx,vy,vz)
        StaggeredAdvect(particle, isf, vx,vy,vz, isf.dt)

        # update gpu buffer
        r = 1:length(particle)
        gpu_position_x[r] = particle.x[r]
        gpu_position_y[r] = particle.y[r]
        gpu_position_z[r] = particle.z[r]
        set_arg!(particle_vis, :indices, map(x-> Cuint(x-1), particle.active))
        yield()
    end
end
@time loop(
    1000, isf, psi1, psi2, dt, kvec, omega, isJet, n_particles,
    nozzle_rad, nozzle_cen, particle, particle_vis
)


text:
function randomwavefields()
    modes = EigenAngle.(rand(TEST_RNG, ComplexF64, 10))
    heights = LMPParams().wavefieldheights

    v = rand(SVector{6,ComplexF64}, length(heights), length(modes))
    wavefields = LMP.Wavefields(v, heights, modes)

    return modes, wavefields
end

function test_Wavefields()
    _, wavefields = randomwavefields()
    @unpack wavefieldheights = LMPParams()

    @test LMP.numheights(wavefields) == length(wavefieldheights)
    @test LMP.nummodes(wavefields) == 10

    @test size(wavefields) == (length(wavefieldheights), 10)
    @test view(wavefields, 1:5) == view(wavefields.v, 1:5)

    newwavefields = similar(wavefields)
    @test size(newwavefields) == size(wavefields)
    @test LMP.heights(newwavefields) == LMP.heights(wavefields)
    @test LMP.eigenangles(newwavefields) == LMP.eigenangles(wavefields)

    newwavefields[1] = wavefields[1]*2
    @test newwavefields[1] == wavefields[1]*2

    # copy and ==
    cpwavefields = copy(wavefields)
    @test cpwavefields.heights == wavefields.heights
    @test cpwavefields.eas == wavefields.eas
    @test cpwavefields.v == wavefields.v
    @test cpwavefields == wavefields

    # isvalid
    heights = LMP.heights(wavefields)
    eas = LMP.eigenangles(wavefields)

    @test isvalid(wavefields)
    longerv = vcat(wavefields.v, rand(TEST_RNG, SVector{6,ComplexF64}, 1, length(eas)))
    badwavefields = LMP.Wavefields(longerv, heights, eas)
    @test !isvalid(badwavefields)

    widerv = hcat(wavefields.v, rand(TEST_RNG, SVector{6,ComplexF64}, length(heights)))
    badwavefields = LMP.Wavefields(widerv, heights, eas)
    @test !isvalid(badwavefields)
end

function test_WavefieldIntegrationParams(scenario)
    @unpack ea, tx, bfield, species = scenario

    params = LMPParams()
    topheight = params.topheight

    w1 = LMP.WavefieldIntegrationParams(topheight, LMP.BOTTOMHEIGHT, 0.0+0.0im, 1.0, 1.0, ea,
        tx.frequency, bfield, species, params)
    w2 = LMP.WavefieldIntegrationParams(topheight, ea, tx.frequency, bfield, species, params)

    @test w1 == w2
end

function test_integratewavefields_homogeneous(scenario)
    #==
    Check wavefields in homogeneous ionosphere are valid solutions to wave equation.
    Compares to Booker quartic solution.
    See, e.g. Pitteway 1965 pg 234; also Barron & Budden 1959 sec 10
    ==#

    @unpack ea, bfield, tx, ground, species = scenario
    params = LMPParams(earthcurvature=false)

    ionobottom = params.curvatureheight
    zs = 200e3:-500:ionobottom

    e = LMP.integratewavefields(zs, ea, tx.frequency, bfield, species; params=params)

    # Normalize fields so component 2 (Ey) = 1, as is used in Booker Quartic
    e1 = [s[:,1]/s[2,1] for s in e]
    e2 = [s[:,2]/s[2,2] for s in e]

    e1 = reshape(reinterpret(ComplexF64, e1), 4, :)
    e2 = reshape(reinterpret(ComplexF64, e2), 4, :)

    # Booker solution - single solution for entire homogeneous iono
    M = LMP.susceptibility(ionobottom, tx.frequency, bfield, species; params=params)
    booker = LMP.bookerwavefields(ea, M)

    @test isapprox(e1[:,end], booker[:,1]; atol=1e-6)
    @test isapprox(e2[:,end], booker[:,2]; atol=1e-6)
end

function test_wavefieldreflection(scenario)
    # Confirm reflection coefficients from wavefields match with dr/dz calculation.

    @unpack ea, bfield, tx, ground, species = scenario
    params = LMPParams()
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    modeequation = PhysicalModeEquation(ea, tx.frequency, waveguide)

    zs = params.wavefieldheights

    e = LMP.integratewavefields(zs, ea, tx.frequency, bfield, species; params=params)
    wavefieldRs = [LMP.bookerreflection(ea, s) for s in e]

    @unpack tolerance, solver = params.integrationparams

    Mtop = LMP.susceptibility(first(zs), tx.frequency, waveguide)
    Rtop = LMP.bookerreflection(ea, Mtop)
    prob = ODEProblem{false}(LMP.dRdz, Rtop, (first(zs), last(zs)), modeequation)
    sol = solve(prob, solver; abstol=tolerance, reltol=tolerance,
                saveat=zs, save_everystep=false)

    @test isapprox(wavefieldRs, sol.u, rtol=1e-3)
end

function test_wavefieldreflection_resonant(scenario)
    @unpack ea, bfield, tx, ground, species = scenario
    params = LMPParams()

    ztop = params.topheight
    zs = ztop:-100:0.0

    e = LMP.integratewavefields(zs, ea, tx.frequency, bfield, species; params=params)
    R = LMP.bookerreflection(ea, e[end])
    Rg = LMP.fresnelreflection(ea, ground, tx.frequency)

    # Ensure we are close to mode resonance with R
    f = LMP.modalequation(R, Rg)

    @test isroot(f)
end

function test_boundaryscalars(scenario)
    # Check if boundaryscalars match for isotropic=false and isotropic=true when both
    # actually are isotropic

    @unpack ea, bfield, tx, ground, species = scenario
    params = LMPParams()

    ztop = params.topheight
    zs = ztop:-100:0.0

    e = LMP.integratewavefields(zs, ea, tx.frequency, bfield, species; params=params)
    R = LMP.bookerreflection(ea, e[end])
    Rg = LMP.fresnelreflection(ea, ground, tx.frequency)

    a1, a2 = LMP.boundaryscalars(R, Rg, e[end], false)
    b1, b2 = LMP.boundaryscalars(R, Rg, e[end], true)

    # return a1, a2, b1, b2

    @test a1 ≈ b1
    @test a2 ≈ b2
end

function test_fieldstrengths(scenario)
    @unpack bfield, tx, ground, species = scenario
    modes, wavefields = randomwavefields()

    modes = LMP.eigenangles(wavefields)

    LMP.fieldstrengths!(view(wavefields,:,1), LMPParams().wavefieldheights, modes[1],
                         tx.frequency, bfield, species, ground)
end

@testset "wavefields.jl" begin
    @info "Testing wavefields"

    test_Wavefields()

    for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario)
        test_WavefieldIntegrationParams(scn)
        @test_nowarn test_fieldstrengths(scn)
    end

    @testset "Integration" begin
        for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario)
            test_wavefieldreflection(scn)
        end
        for scn in (resonant_scenario, )
            test_wavefieldreflection_resonant(scn)
        end
        for scn in (homogeneousiono_scenario, )
            test_integratewavefields_homogeneous(scn)
        end

        @test_skip test_boundaryscalars(isotropicB_resonant_scenario)
    end
end


text:
dir = "./src"
push!(LOAD_PATH,dir)

module physicalConstants

export k, g, p, MM

# Boltzmann constant (k) SI units: m2 kg s-2 K-1
const k = 1.38064852e-23

# Standard gravity (g) SI units: m s-2
const g = 9.80665

# Air at atmospheric pressure (Pa) SI units: kg s-2 m-2
const p = 101325;

# Molar mass of air (kg mol^-1) SI units: kg mol-1
const MM = 28.97/(1000*6.02214179e23);

end


text:
module HaftaViz
using Hafta
using Gadfly
using Formatting

# writemime functions for HFB
function Base.writemime(io, ::MIME"text/html", state::Hafta.HFB.HFBState)
    width,height = 10cm,8cm

    write(io, "
")
    write(io, "
HFBState $(size(state.U))
")

    # Table of energies and other values
    E, Ek, Ei, Ep = energy(state)
    Aest = trace(state.rho)
    html = """

        Energy
        Free particles E
        Interaction energy (Γ)
        Pairing energy (Δ)
        Particle number
        Chem.pot.

        {:.5f}
        {:.5f}
        {:.5f}
        {:.5f}
        {:.5f}
        {:.5f}

    """
    write(io,format(html, E, Ek, Ei, Ep, Aest, state.lambda))

    # rho and kappa matrices
    rho,kappa = abs(state.rho), abs(state.kappa)
    maxz = max(maximum(rho), maximum(kappa))
    scale = Scale.color_continuous(minvalue=0, maxvalue=maxz)

    write(io, "
")
    write(io, "")
    p = spy(rho, Guide.title("rho matrix"), scale)
    draw(SVG(io,width,height,false), p)
    write(io, "	")
    write(io, "")
    p = spy(kappa, Guide.title("kappa matrix"), scale)
    draw(SVG(io,width,height,false), p)
    write(io, "	")
    write(io, "
")

    # U and V matrices
    U,V = abs(state.U), abs(state.V)
    maxz = max(maximum(U), maximum(V))
    scale = Scale.color_continuous(minvalue=0, maxvalue=maxz)

    write(io, "
")
    write(io, "")
    p = spy(abs(state.U), Guide.title("U matrix"), scale)
    draw(SVG(io,width,height,false), p)
    write(io, "	")
    write(io, "")
    p = spy(abs(state.V), Guide.title("V matrix"), scale)
    draw(SVG(io,width,height,false), p)
    write(io, "	")
    write(io, "
")

    # the end
    write(io, "
")
end

function Base.writemime(io, ::MIME"text/html", hfbi::Hafta.HFB.HFBIterator)
    width,height = 20cm,6cm

    write(io, "
")
    write(io, """

            HFBIterator ($(length(hfbi)) iterations)

""")

    write(io, "
")
    p = plot(
        x=1:length(hfbi.es), y=hfbi.es,
        Geom.line,# Geom.point,
        Guide.title("Absolute value of energy"),
        Guide.xlabel("Iteration"), Guide.ylabel("E")
    )
    draw(SVG(io,width, height,false), p)
    #write(io, "
")
    write(io, "<br /")

    #write(io, "
")
    logdiffs = log10(abs(diff(hfbi.es)))
    p = plot(
        x=1:length(logdiffs), y=logdiffs,
        yintercept=[-10, -15], Geom.hline(color="red"),
        Geom.line,# Geom.point,
        Guide.title("Differences of consecutive energies"),
        Guide.xlabel("Iteration"), Guide.ylabel("log10(ΔE)")
    )
    draw(SVG(io,width, height,false), p)
    write(io, "
")

    # Output the first and the last state
    if length(hfbi.states) > 1
        write(io, """

        Final state

""")
        write(io, "
")
        writemime(io, "text/html", hfbi.states[end])
        write(io, "
")
    end

    if false
        write(io, """

                Initial state

""")
        write(io, "
")
        writemime(io, "text/html", hfbi.states[1])
        write(io, "
")
    end

    write(io, "
")
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineIdentityUserAssignedIdentities <: SwaggerModel
    principalId::Any # spec type: Union{ Nothing, String } # spec name: principalId
    clientId::Any # spec type: Union{ Nothing, String } # spec name: clientId

    function VirtualMachineIdentityUserAssignedIdentities(;principalId=nothing, clientId=nothing)
        o = new()
        validate_property(VirtualMachineIdentityUserAssignedIdentities, Symbol("principalId"), principalId)
        setfield!(o, Symbol("principalId"), principalId)
        validate_property(VirtualMachineIdentityUserAssignedIdentities, Symbol("clientId"), clientId)
        setfield!(o, Symbol("clientId"), clientId)
        o
    end
end # type VirtualMachineIdentityUserAssignedIdentities

const _property_map_VirtualMachineIdentityUserAssignedIdentities = Dict{Symbol,Symbol}(Symbol("principalId")=>Symbol("principalId"), Symbol("clientId")=>Symbol("clientId"))
const _property_types_VirtualMachineIdentityUserAssignedIdentities = Dict{Symbol,String}(Symbol("principalId")=>"String", Symbol("clientId")=>"String")
Base.propertynames(::Type{ VirtualMachineIdentityUserAssignedIdentities }) = collect(keys(_property_map_VirtualMachineIdentityUserAssignedIdentities))
Swagger.property_type(::Type{ VirtualMachineIdentityUserAssignedIdentities }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineIdentityUserAssignedIdentities[name]))}
Swagger.field_name(::Type{ VirtualMachineIdentityUserAssignedIdentities }, property_name::Symbol) =  _property_map_VirtualMachineIdentityUserAssignedIdentities[property_name]

function check_required(o::VirtualMachineIdentityUserAssignedIdentities)
    true
end

function validate_property(::Type{ VirtualMachineIdentityUserAssignedIdentities }, name::Symbol, val)
end


text:
# tensorcopy.jl
#
# Method for copying one tensor to another according to the
# specified labels, thereby possibly having to permute the
# data.

# Simple method
# --------------
function tensorcopy(A::StridedArray,labelsA,outputlabels=labelsA)
    dims=size(A)
    C=similar(A,dims[indexin(outputlabels,labelsA)])
    tensorcopy!(A,labelsA,C,outputlabels)
end

# In-place method
#-----------------
function tensorcopy!(A::StridedArray,labelsA,C::StridedArray,labelsC)
    NA=ndims(A)
    perm=indexin(labelsC,labelsA)
    length(perm) == NA || throw(LabelError("invalid label specification"))
    isperm(perm) || throw(LabelError("labels do not specify a valid permutation"))
    for i = 1:NA
        size(A,perm[i]) == size(C,i) || throw(DimensionMismatch("destination tensor of incorrect size"))
    end
    NA==0 && (C[1]=A[1]; return C)
    perm==collect(1:NA) && return copy!(C,A)
    tensorcopy_native!(A,C,perm)
end

# Implementation
#----------------
const PERMUTEGENERATE=[1,2,3,4,5,6,7,8]

@ngenerate N typeof(C) function tensorcopy_native!{TA,TC,N}(A::StridedArray{TA,N},C::StridedArray{TC,N},perm)
    @nexprs N d->(stridesA_{d} = stride(A,perm[d]))
    @nexprs N d->(stridesC_{d} = stride(C,d))
    @nexprs N d->(dims_{d} = size(C,d))

    startA = 1
    local Alinear::Array{TA,N}
    if isa(A, SubArray)
        startA = A.first_index
        Alinear = A.parent
    else
        Alinear = A
    end
    startC = 1
    local Clinear::Array{TC,N}
    if isa(C, SubArray)
        startC = C.first_index
        Clinear = C.parent
    else
        Clinear = C
    end

    if length(C)<=4*TBASELENGTH
        @stridedloops(N, i, dims, indA, startA, stridesA, indC, startC, stridesC, @inbounds Clinear[indC]=Alinear[indA])
    else
        @nexprs N d->(minstrides_{d} = min(stridesA_{d},stridesC_{d}))

        # build recursive stack
        depth=ceil(Integer, log2(length(C)/TBASELENGTH))+2 # 2 levels safety margin
        level=1 # level of recursion
        stackpos=zeros(Int,depth) # record position of algorithm at the different recursion level
        stackpos[level]=0
        stackblength=zeros(Int,depth)
        stackblength[level]=length(C)
        @nexprs N d->begin
            stackbdims_{d} = zeros(Int,depth)
            stackbdims_{d}[level] = dims_{d}
        end
        stackbstartA=zeros(Int,depth)
        stackbstartA[level]=startA
        stackbstartC=zeros(Int,depth)
        stackbstartC[level]=startC

        stackdC=zeros(Int,depth)
        stackdA=zeros(Int,depth)
        stackdmax=zeros(Int,depth)
        stacknewdim=zeros(Int,depth)
        stackolddim=zeros(Int,depth)

        while level>0
            pos=stackpos[level]
            blength=stackblength[level]
            @nexprs N d->(bdims_{d} = stackbdims_{d}[level])
            bstartA=stackbstartA[level]
            bstartC=stackbstartC[level]

            if blength<=TBASELENGTH || level==depth # base case
                @stridedloops(N, i, bdims, indA, bstartA, stridesA, indC, bstartC, stridesC, @inbounds Clinear[indC]=Alinear[indA])
                level-=1
            elseif pos==0
                # find which dimension to divide
                dmax=0
                maxval=0
                newdim=0
                olddim=0
                dC=0
                dA=0
                @nexprs N d->begin
                    newmax=bdims_{d}*minstrides_{d}
                    if bdims_{d}>1 && newmax>maxval
                        dmax=d
                        olddim=bdims_{d}
                        newdim=olddim>>1
                        dC=stridesC_{d}
                        dA=stridesA_{d}
                        maxval=newmax
                    end
                end
                stackolddim[level]=olddim
                stacknewdim[level]=newdim
                stackdmax[level]=dmax
                stackdC[level]=dC
                stackdA[level]=dA

                stackpos[level+1]=0
                @nexprs N d->(stackbdims_{d}[level+1] = (d==dmax ? newdim : bdims_{d}))
                stackblength[level+1]=div(blength,olddim)*newdim
                stackbstartA[level+1]=bstartA
                stackbstartC[level+1]=bstartC

                stackpos[level]+=1
                level+=1
            elseif pos==1
                olddim=stackolddim[level]
                newdim=stacknewdim[level]
                dmax=stackdmax[level]
                dC=stackdC[level]
                dA=stackdA[level]

                stackpos[level+1]=0
                @nexprs N d->(stackbdims_{d}[level+1] = (d==dmax ? olddim-newdim : bdims_{d}))
                stackblength[level+1]=div(blength,olddim)*(olddim-newdim)
                stackbstartA[level+1]=bstartA+dA*newdim
                stackbstartC[level+1]=bstartC+dC*newdim

                stackpos[level]+=1
                level+=1
            else
                level-=1
            end
        end
    end
    return C
end


text:
# A variety of means

# Geometric mean
function geomean(a::RealArray)
    s = 0.0
    n = length(a)
    for i = 1 : n
        @inbounds s += log(a[i])
    end
    return exp(s / n)
end

# Harmonic mean
function harmmean(a::RealArray)
    s = 0.0
    n = length(a)
    for i in 1 : n
        @inbounds s += inv(a[i])
    end
    return n / s
end

# Trimmed mean

function trimmean(x::RealArray, p::FloatingPoint)
    n = length(x)
    n > 0 || error("x can not be empty.")
    rn = min(iround(n * p), n-1)
    if rn == n - 1
        return median(x)
    else
        sx = sort(x)
        nl = rn >> 1
        nh = (rn - nl)
        s = 0.0
        @inbounds for i = (1+nl) : (n-nh)
            s += x[i]
        end
        return s / (n - rn)
    end
end

# Weighted mean
function wmean{T<:Number,W<:Real}(v::AbstractArray{T}, w::AbstractArray{W}; wsum::W=NaN)
	# wsum is the pre-computed sum of weights

	n = length(v)
	length(w) == n || raise_dimerror()
    sv = zero(T)

    if isnan(wsum)
    	sw = zero(W)
    	for i = 1 : n
    		@inbounds wi = w[i]
        	@inbounds sv += v[i] * wi
        	sw += wi
    	end
    else
    	for i = 1 : n
    		@inbounds sv += v[i] * w[i]
    	end
    	sw = wsum
    end

    return sv / sw
end

# Weighted mean (fast version for contiguous arrays using BLAS)
function wmean{T<:BlasReal}(v::Array{T}, w::Array{T}; wsum::T=NaN)
	sw = (isnan(wsum) ? sum(w) : wsum)
	dot(v, w) / sw
end

mean{T<:Number}(v::AbstractArray{T}, w::WeightVec) = wmean(v, values(w); wsum=sum(w))


text:
############################################################################
# imports.jl
#

import Base: *

############################################################################
# NDTensors
#

# Take the square root of a diagonal order 2 tensor
# TODO: make a AbstractDiagTensor supertype?
# TODO: make this work for non-diagonal tensors.
# TODO: define SecondOrderTensor?
function Base.sqrt(T::NDTensors.Tensor{<:Number,2}; atol=1e-15)
  @assert isdiag(T)
  sqrtT = copy(T)
  for n in 1:mindim(T)
    Tnn = T[n, n]
    if Tnn < 0 && abs(Tnn) < atol
      sqrtT[n, n] = 0
    else
      sqrtT[n, n] = sqrt(Tnn)
    end
  end
  return sqrtT
end

############################################################################
# tagset.jl
#

# Used for findfirst
Base.keys(ts::TagSet) = Base.OneTo(length(ts))

import ITensors: Tag, commontags

macro tag_str(s)
  return Tag(s)
end

commontags(ts::TagSet, i::Index) = commontags(ts, tags(i))

commontags(t::Vector) = foldl(commontags, t)

maxlength(tag::Tag) = length(tag.data)

function Base.length(tag::Tag)
  n = 1
  while n <= maxlength(tag) && tag[n] != zero(eltype(tag))
    n += 1
  end
  return n - 1
end

Base.lastindex(tag::Tag) = length(tag)

Base.getindex(tag::Tag, r::UnitRange) = Tag([tag[n] for n in r])

# TODO: make this work directly on a Tag, without converting
# to String
function Base.parse(::Type{T}, tag::Tag) where {T<:Integer}
  return parse(T, string(tag))
end

function Base.startswith(tag::Tag, subtag::Tag)
  for n in 1:length(subtag)
    tag[n] ≠ subtag[n] && return false
  end
  return true
end

function tag_starting_with(ts::TagSet, prefix)
  x = findfirst(t -> startswith(t, Tag(prefix)), ts)
  if isnothing(x) #in case the function is called one a link leg. Can be probably improved
    return x
  end
  return ts[findfirst(t -> startswith(t, Tag(prefix)), ts)]
end

ITensors.replacetags(ts::TagSet, p::Pair) = replacetags(ts, first(p), last(p))

############################################################################
# qn.jl
#

#
# New QN interface.
# Allows constructors like:
# Index((("Sz", 0) × ("SzParity", 0, 2) => 1) ⊕
#       (("Sz", 1) × ("SzParity", 1, 2) => 2))
#

function LinearAlgebra.cross(
  qn1::Tuple{String,Int,Vararg{Int}}, qn2::Tuple{String,Int,Vararg{Int}}
)
  return QN(qn1, qn2)
end
function LinearAlgebra.cross(qn1::QN, qn2::Tuple{String,Int,Vararg{Int}})
  return ITensors.addqnval(qn1, ITensors.QNVal(qn2...))
end

⊕(qn1::Pair{QN,Int}, qn2::Pair{QN,Int}) = [qn1, qn2]
function ⊕(qn1::Pair{QN,Int64}, qn2::Pair{<:Tuple{String,Int,Vararg{Int}},Int})
  return [qn1, QN(first(qn2)) => last(qn2)]
end
function ⊕(qn1::Pair{<:Tuple{String,Int,Vararg{Int}},Int}, qn2::Pair{QN,Int})
  return [QN(first(qn1)) => last(qn1), qn2]
end
function ⊕(
  qn1::Pair{<:Tuple{String,Int,Vararg{Int}},Int},
  qn2::Pair{<:Tuple{String,Int,Vararg{Int}},Int},
)
  return [QN(first(qn1)) => last(qn1), QN(first(qn2)) => last(qn2)]
end
⊕(qns::Vector{Pair{QN,Int}}, qn::Pair{QN,Int}) = push!(copy(qns), qn)

############################################################################
# indexset.jl
#

ITensors.IndexSet(is::IndexSet...) = unioninds(is)
ITensors.IndexSet(is::Tuple{Vararg{<:IndexSet}}) = unioninds(is...)

Base.copy(is::IndexSet) = IndexSet(copy.(ITensors.data(is)))

ITensors.noncommoninds(is::IndexSet) = is

############################################################################
# itensor.jl
#

using ITensors.NDTensors

function LinearAlgebra.ishermitian(T::ITensor, pairs=0 => 1; kwargs...)
  Tᴴ = swapprime(dag(T), pairs)
  return isapprox(Tᴴ, T; kwargs...)
end

# Helpful for making sure the ITensor doesn't contract
ITensors.sim(A::ITensor) = ITensors.setinds(A, sim(inds(A)))

LinearAlgebra.isdiag(T::ITensor) = isdiag(tensor(T))

function eigendecomp(T::ITensor; ishermitian=true, kwargs...)
  @assert ishermitian
  D, U = eigen(T; ishermitian=ishermitian, kwargs...)
  return U', D, dag(U)
end

*(A::ITensor) = A

ITensors.noncommoninds(A::ITensor) = inds(A)

# TODO: implement something like this?
#function sqrt(::Order{2}, T::ITensor)
#end

# Take the square root of T assuming it is Hermitian
# TODO: add more general index structures
function Base.sqrt(T::ITensor; ishermitian=true, atol=1e-15)
  @assert ishermitian
  if isdiag(T) && order(T) == 2
    return itensor(sqrt(tensor(T)))
  end
  U′, D, Uᴴ = eigendecomp(T; ishermitian=ishermitian)
  # XXX: if T is order 2 and diagonal, D may just be a view
  # of T so this would also modify T
  #D .= sqrt.(D)
  sqrtD = D
  for n in 1:mindim(D)
    Dnn = D[n, n]
    if Dnn < 0 && abs(Dnn) < atol
      sqrtD[n, n] = 0
    else
      sqrtD[n, n] = sqrt(Dnn)
    end
  end
  return U′ * sqrtD * Uᴴ
end

############################################################################
# mps.jl
#

# TODO: make this definition AbstractMPS
# Handle orthogonality center correctly
Base.getindex(ψ::MPS, r::UnitRange{Int}) = MPS([ψ[n] for n in r])


text:
@testset "Utilities" begin
    t = collect(-2:0.01:2)
    U = [cos.(t).*exp.(-t.^2) sin.(2*t)]
    S = Diagonal([2.; 3.])
    V = [sin.(t).*exp.(-t) cos.(t)]
    A = U*S*V'
    σ = 0.5
    Â = A + σ*randn(401, 401)
    n_1 = norm(A-Â)
    B = optimal_shrinkage(Â)
    optimal_shrinkage!(Â)
    @test norm(A-Â) < n_1
    @test norm(A-B) == norm(A-Â)

    X = randn(3, 100)
    Y = randn(3, 100)
    k = 3

    @test AIC(k, X, Y) == 2*k-2*log(sum(abs2, X- Y))
    @test AICC(k, X, Y) == AIC(k, X, Y)+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)
    @test BIC(k, X, Y) == -2*log(sum(abs2, X -Y)) + k*log(size(X)[2])
    @test AICC(k, X, Y, likelihood = (X,Y)->sum(abs, X-Y)) == AIC(k, X, Y, likelihood = (X,Y)->sum(abs, X-Y))+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)

    # Numerical derivatives
    function lorenz(u,p,t)
        x, y, z = u
        ẋ = 10.0*(y - x)
        ẏ = x*(28.0-z) - y
        ż = x*y - (8/3)*z
        return [ẋ, ẏ, ż]
    end
    u0 = [1.0;0.0;0.0]
    tspan = (0.0,50.0)
    dt = 0.005
    prob = ODEProblem(lorenz,u0,tspan)
    sol = solve(prob, Tsit5(), saveat = dt)

    X = Array(sol)
    DX = similar(X)
    for (i, xi) in enumerate(eachcol(X))
        DX[:,i] = lorenz(xi, [], 0.0)
    end

    windowSize, polyOrder = 9, 4
    halfWindow = Int(ceil((windowSize+1)/2))
    DX_sg = similar(X)
    DX_sg_cropped = similar(X[:,halfWindow+1:end-halfWindow])
    X_cropped = similar(X[:,halfWindow+1:end-halfWindow])
    for i =1:size(X,1)
        DX_sg[i,:] = savitzky_golay(X[i,:], windowSize, polyOrder, deriv=1, dt=dt, crop=false)
        X_cropped[i,:], DX_sg_cropped[i,:] = savitzky_golay(X[i,:], windowSize, polyOrder, deriv=1, dt=dt)
    end
    DX_sg2 = savitzky_golay(X, windowSize, polyOrder, deriv=1, dt=dt, crop=false)
    X_cropped2, DX_sg_cropped2 =  savitzky_golay(X, windowSize, polyOrder, deriv=1, dt=dt)
    @test(DX_sg2 == DX_sg)
    DX_sg = DX_sg[:,halfWindow+1:end-halfWindow]
    @test X_cropped == X[:,halfWindow+1:end-halfWindow]
    @test DX_sg_cropped == DX_sg
    @test X_cropped2 == X[:,halfWindow+1:end-halfWindow]
    @test DX_sg_cropped2 == DX_sg

    DX = DX[:,halfWindow+1:end-halfWindow]
    @test isapprox(DX_sg, DX, rtol=1e-2)

    # Sampling
    X = randn(Float64, 2, 100)
    t = collect(0:0.1:9.99)
    Y = randn(size(X))
    xt = burst_sampling(X, 5, 10)
    @test 10 <= size(xt)[end] <= 60
    @test all([any(xi .≈ X) for xi in eachcol(xt)])
    xt, tt = burst_sampling(X, t, 5, 10)
    @test all(diff(tt) .> 0.0)
    @test size(xt)[end] == size(tt)[end]
    @test all([any(xi .≈ X) for xi in eachcol(xt)])
    @test !all([any(xi .≈ Y) for xi in eachcol(xt)])
    xs, ts = burst_sampling(X, t, 2.0, 1)
    @test all([any(xi .≈ X) for xi in eachcol(xs)])
    @test size(xs)[end] == size(ts)[end]
    @test ts[end]-ts[1] ≈ 2.0
    X2n = subsample(X, 2)
    t2n = subsample(t, 2)
    @test size(X2n)[end] == size(t2n)[end]
    @test size(X2n)[end] == Int(round(size(X)[end]/2))
    @test X2n[:, 1] == X[:, 1]
    @test X2n[:, end] == X[:, end-1]
    @test all([any(xi .≈ X) for xi in eachcol(X2n)])
    xs, ts = subsample(X, t, 0.5)
    @test size(xs)[end] == size(ts)[end]
    @test size(xs)[1] == size(X)[1]
    @test all(diff(ts) .≈ 0.5)
    # Loop this a few times to make sure it's right
    @test_nowarn for i in 1:20
        xs, ts = burst_sampling(X, t, 2.0, 1)
        xs, ts = subsample(X, t, 0.5)
    end
end


text:
module VirtualMachine

export QVM, showstate, step!, execute!

using ..Utils
using ..QuantInstr
using LinearAlgebra

⊗(x, y) = kron(x, y)

const GATES = Dict(:I => () -> [1 0;
                                0 1],

                   :X => () -> [0 1;
                                1 0],
                   :RX => γ -> exp(im * γ * GATES[:X]()),

                   :Y => () -> [0 -im;
                                im 0],
                   :RY => β -> exp(im * β * GATES[:Y]()),

                   :Z => () -> [1  0;
                                0 -1],
                   :RZ => α -> exp(im * α * GATES[:Z]()),

                   :H => () -> [1  1;
                                1 -1] / sqrt(2),

                   :S => () -> [1 0;
                                0 im],

                   :T => () -> [1 0;
                                0 exp(im * π / 4)],

                   :PHASE => θ -> [1 0; 0 exp(im * θ)],

                   :SWAP => () -> [1 0 0 0;
                                   0 0 1 0;
                                   0 1 0 0;
                                   0 0 0 1],
                  )

# Type aliases
const Quip = Vector{QInstr}

"""
The Quantum Virtual Machine,
    a Turing machine that holds the instruction set
    as well as the current state of the wave function
"""
struct QVM
    # Stacks, each element is a stack frame
    stack::Vector{Quip}
    pc_stack::Vector{Int}

    # Classical memory
    reg::Vector{Bool}

    # Wavefunction
    wfn::Vector{Complex{Float64}}

    # A dictionary of primitive operations
    dict::Dict{Symbol, Function}

    function QVM(quip::Vector{QInstr}, nbit::Int;
            dict=GATES, nreg=32)

        wfn = zeros(Complex{Float64}, 2^nbit)
        wfn[1] = 1

        registers = zeros(Int, nreg)
        new([quip], [1], registers, wfn, dict)
    end
end

"""
pretty print the current stack
for internal use by showstate
"""
function showstack(stack::Vector{QInstr}, pc::Int, n::Int; io::IO=stdout)
    len = length(stack)

    # Checking boundaries
    if pc - n < 1
        s = 1
        e = 1 + 2n > len ? len : 1 + 2n
    elseif pc + n > len
        e = len
        s = 1 - 2n < 1 ? 1 : len - 2n
    else
        s = pc - n
        e = pc + n
    end

    # Print stack
    if s != 1
        println(io, "...")
    end

    for i in s:e
        print(io, "$i: $(stack[i])")
        if i == pc
            println(io, "   <--- pc")
        else
            println(io)
        end
    end

    if e != len
        println(io, "...")
    end

    if pc > len
        println(io, "$(len + 1):  <- pc")
    end
end

function showstate(vm::QVM, n=2; io::IO=stdout)
    println(io, "##########################")
    for i in eachindex(vm.stack)
        stack = vm.stack[i]
        pc = vm.pc_stack[i]

        println(io, "$i ---------------------------")
        showstack(stack, pc, n, io=io)
        println(io)
    end
end

function step!(vm::QVM)
    dict = vm.dict
    quip = vm.stack[end]
    pc = vm.pc_stack[end]

    ################################################
    # PC reached the end of the current stack frame
    #
    if pc > length(quip)
        # VM reached the end of program do nothing
        if length(vm.stack) == 1
            return false
        end

        # pop stack frame
        pop!(vm.stack)
        pop!(vm.pc_stack)
        vm.pc_stack[end] += 1 # increment pc
        return true
    end

    instr = quip[pc]
    op = instr.op
    func = instr.func
    param = instr.param
    target = instr.target
    control = instr.control
    cbit_ctrl = instr.cbit_ctrl

    #####################################
    # Function call push new stack frame
    #
    if op == :funcall
        new_quip = func(param...)
        push!(vm.stack, new_quip)
        push!(vm.pc_stack, 1)

        ##########################################################
        # Measurement, collapse wave function and save to register
        #
    elseif op == :measure
        qbit = target[1]
        cbit = target[2]
        # Bases
        bas0 = [1; 0]
        bas1 = [0; 1]

        wfn0 = deepcopy(vm.wfn)
        applygate!(wfn0, bas0 * bas0', [qbit])
        P0 = vm.wfn' * wfn0

        # collapse to 0 state
        if rand() < norm(P0)
            vm.reg[cbit] = false
            vm.wfn .= wfn0 / sqrt(P0)

            # collapse to 1 state
        else
            wfn1 = vm.wfn - wfn0
            P1 = 1 - P0

            vm.reg[cbit] = true
            vm.wfn .= wfn1 / sqrt(P1)
        end
        vm.pc_stack[end] += 1 # increment pc

        ##################################
        # Gate operation, call applygate!
        #
    else
        operator = get(dict, op) do
            throw("$op not in dict, this should not happen")
        end

        # cBit control
		if cbit_ctrl
			conds = [vm.reg[abs(i)] ⊻ (i < 0) for i in control]
			if reduce((a, b) -> a && b, conds)
				applygate!(vm.wfn, operator(param...), target)
			end

            # qBit control
        else
            applygate!(vm.wfn, operator(param...), target, control=control)
        end
        vm.pc_stack[end] += 1 # increment pc
    end

    true
end

"""
Run step!(vm) until pc reaches the end of the program strip
"""
function execute!(vm::QVM; verbose::Bool=false, io::IO=stdout)
    if verbose
        showstate(vm, io=io)
    end
    while step!(vm)
        if verbose
            showstate(vm, io=io)
        end
    end
end

end


text:
function Ulp()
    return Ulp(())
end

function value(obj::Ulp, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
# using PyPlot

function fit_weibul(data)

    y = log(sort(data))
    rank = 1:length(data)  # ranks = {1, 2, 3, ... 10}
    median_rank = (rank - 0.3)/(length(rank) + 0.4)
    x = log(-log(1 - median_rank))

    intercept, slope = linreg(x,y)

    shape = 1/slope
    x_intercept = - intercept / shape
    scale = exp(-x_intercept/slope)

#     @show(shape, scale)
    return Weibull(shape, scale)
end

function f2(f)
    @sprintf("%.2f", f)
end

## Array indexes

≈(x,y::Real) = isapprox(x, y)
≈(x,yt::Tuple) = isapprox(x, yt[1], atol=yt[2])
≈(xs::Array,y::Int) = find(x-> trunc(x) == y, xs)
≈(xs::Array,y::Real) = find(x->isapprox(x, y), xs)
≈(xs::Array,yt::Tuple) = find(x->isapprox(x, yt[1], atol=yt[2]), xs)

## Dictionary Handling ##

markdownsplitter(x::String) = @pipe x |> strip(_) |> strip(_,'|') |> split(_,'|') |> map(strip,_)
dictFromArray(fields) = x -> OrderedDict(zip(fields,x))
arrayFromDictKeyArray(d::Dict, keys::Array, exclude...) = [ d[k] for k in filter(x-> !(x in exclude),keys) ]

⊂(d::Dict, keys::Array) = Dict( [ k=>v for (k,v) in filter((x,y)-> in(x,Set(keys)),d) ] )
∪(u::Dict, v::Dict) = Dict( concat(u ⊂ collect(keys(v)),v) )

## Plotting Configs ##

function configure_matplotlib_publication_style(rc)

    ## NSF Grant Settings
    rc("axes", linewidth=2)                                # but this works

    rc("font", size=18, family="serif", variant="normal", weight="heavy", serif=["Times New Roman"])
    rc("text", color="k")
    rc("axes", facecolor="white", color_cycle="k", labelsize= 16, labelcolor= "k", linewidth= 1, edgecolor="k")

    rc("figure", facecolor="white")
    rc("figure", figsize=(6,4))
    rc("savefig", facecolor="white", edgecolor="white", dpi=300)

    rc("figure.subplot", wspace=0.4)

    rc("xtick.major", pad=5)
    rc("ytick.major", pad=5)

    rc("legend", frameon=true, handlelength=1.5)
    rc("text", usetex=false)

    rc("grid", linestyle="-", color="lightgrey" )
    rc("axes", grid = true)
    rc("axes", labelweight="normal" )
    rc("axes.formatter", use_mathtext=false )

    rc("font", weight="normal")
    rc("xtick", labelsize=15)
    rc("ytick", labelsize=15)
    rc("legend", scatterpoints = 1)

	rc("svg.fonttype", "none")

    return;
end

using Formatting

f1(f) = fmt(FormatSpec("0.1f"), f)
f2(f) = fmt(FormatSpec("0.2f"), f)
f3(f) = fmt(FormatSpec("0.3f"), f)
e2(f) = fmt(FormatSpec(".2e"), f)
e1(f) = fmt(FormatSpec(".1e"), f)

function ff(f::Real, n::Int=1)
    abs(f) < 1/10^n || abs(f) > 10^n ? fmt(FormatSpec(".$(n)e"), f) : fmt(FormatSpec(".$(n)f"), f)
end

type HTML
   s::String
end

import Base.writemime

function writemime(io::IO, ::MIME"text/html", x::HTML)
    write(io, x.s);
end

simplify_mat(d::Dict{ASCIIString, Any}) = Dict{Symbol, Any}( [ symbol(string(k)) => simplify_mat(v) for (k,v) in d ] )
simplify_mat(arr::Array{Float64,2}) = arr[:] # slicing???
simplify_mat(o::Any) = o

function convert_slices(matdata)
    # @show(keys(matdata))
    @assert Set(map(x->size(x,1), matdata[:data] |> values)) |> length == 1

    map(matdata[:indexes] |> keys) do key
        L = matdata[:data][key] |> length

        map(matdata[:indexes][key]) do idx
            @match idx begin
                (:idx_neg_1, [i j k]) => (:idx_end, i:k:L)
                (         s, [i -1 k]) => (s, i+1:k:L)
                (         s, [i  j k]) => (s, i+1:k:j+1)
            end
        end |> sort |> OrderedDict
    end |> Dict
end

function load_test_data_raw(info::Dict, testspath::String; stage="norm", method="cycles")
    folder = "$testspath/$name"

    # find the files
    trends_file = glob("data/data (*$stage*trends*$method*).mat", folder) |> first;
    tracks_file = glob("data/data (*$stage*tracking*$method*).mat", folder) |> first;

    # read the files
    trends_mat = matread(trends_file);
    tracks_mat = matread(tracks_file);

    return { :info => info, :trends => trends_mat, :tracks => tracks_mat }
end

function simplify_test_data(testdata::Dict)
    Dict( [ symbol(string(k)) => simplify_mat(v) for (k,v) in testdata])
end

function load_test_data(name::String, tests::String; stage="norm", method="cycles")
    test_data = load_test_data_raw(name, tests) |> simplify_test_data

    trends_slices = mat_slices(test_data[:trends])
    tracks_slices = mat_slices(test_data[:tracks])

    return test_data ∪ {
        :trends => ( test_data[:trends] ∪ {:indexes => trends_slices} ),
        :tracks => ( test_data[:tracks] ∪ {:indexes => tracks_slices} ),
    }
end

function simple_fit(x, y)

    fit = GLM.lm([ones(x) x], y)
    b,a = fit.pp.beta0;

    fit_r2 = r2(y, fit);

    return { :fit => fit, :r2 => fit_r2, :a => a, :b => b, :f => (x -> a*x + b),
        :fstr => "\$$( ff(a,1)) \\,x + $(ff(b,1))\$"  }
end

function r2(ydata, fit)
    1-sumabs2(GLM.residuals(fit))/sumabs2(ydata - mean(ydata)) # or 1-var(residuals(mod))/var(y)
end

## HTML Text ###
th(i) = "$i	"
th(s::Array) = map(th, s)
th(s...) = map(th, s)

td(i) = "$i	"
td(s::Array) = map(td, s)
td(s...) = map(td, s)

tr(i) = "
$i
"
tr(s::Dict) = [ tr(k)*tr(v) for (k,v) in s ]
tr(s::Array) = map(tr, s)
tr(s...) = map(tr, s)

## Test Helpers ##
join_nt(s) = "\t"*join(s, "\n\t")*"\n"
join_ss(s) = join(s, " ")

maketablerows(s::Array) = tr( [ join( [ td(i) for i in j ], "") for j in s ] ) |> joint_nt

tabulate(s; headers::Array=[]) = HTML("""

$( isempty(headers) ? "" : tr( [ th(i...) for i in headers] |> join_ss ) |> joint_nt )
$( maketablerows(s) )

""")

## Data Handling ##

## Research Dir
load_tests_from_md(f::Function, testfile) = load_tests_from_md(testfile, f=f)

function load_md_table(testfile; f::Function=(x)->x)
    tests = open(testfile) do f readlines(f) end
    headers = tests[1] |> markdownsplitter |> xs->map(lowercase, xs) |> xs->map(symbol, xs)
    @show headers

    @pipe( tests[3:end]
        |> map(markdownsplitter,_)
        |> filter( l -> any(  x -> !isempty(x), l), _)
        |> map(dictFromArray(headers),_)
        |> map(f,_)
    )
end


text:
include("../common.jl")

# Offset to add to the version number.  Remember to always bump this.
version_offset = v"0.2.2"

# Minimum Julia version supported: this is important to decide which versions of
# the dependencies to use, in particular the JLL stdlibs.
min_julia_version = v"1.8"

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, configure(version_offset, min_julia_version)...; julia_compat="1.8")


text:
module MOTProject

# Remember to ]dev MOTProject to "add" it while developing
# ]activate MOTProject - add pkgs to module
# Base.compilecache(Base.PkgId(MOTProject)) - recompile without restarting julia

# Using pkgs for the whole module
using StaticArrays
using VideoIO, ImageView
using Plots: Gray, @layout
using Distributions, LinearAlgebra
using DataStructures, ImageFiltering, SpecialFunctions
using Clustering
using Distances

# Include my support files with custome functions and objects and stuff
include("MOTobjects.jl")
include("ImageSim.jl")
include("Preprocess.jl")
include("BackgroundEstimation.jl")
include("ClusterMethods.jl")
include("Filters.jl")
include("Pruning.jl")

# Export specific function (dont have to call MOTProject.myfun to use)
export
    Detection,
    ImageSim,
    SigmaClipping, PixelThreshold, HierarchicalAgglomerative,
    Preprocess,
    KF, KF!, TOMHT, TOMHT!,
    MBest

end # end module


text:
export RanklibRanker

@with_kw mutable struct RanklibRanker <: BaseEstimator
    ranklib::Vector{UInt8} = UInt8[]
    ranker::Int = 0
    metric2t::String = "MAP"
end

is_classifier(::RanklibRanker) = false
is_ranker(::RanklibRanker) = true

function gridparams(m::RanklibRanker)
    grid = [
        "ranker" => [0:4; 6:8],
        "metric2t" =>
        ["MAP"; ["$s@$k" for k in [20, 50, 100]
        for s in split("NDCG DCG P RR ERR")]]
    ]
    filter(gridparams(grid)) do d
        d["metric2t"] == "MAP" || d["ranker"] ∈ [3, 4, 6]
    end
end

function fit!(m::RanklibRanker, x, y, w = nothing; group = nothing, columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    if isnothing(group) && ndims(x) == 3
        group = repeat([size(x, 2)], size(x, 3))
    end
    @unpack ranker, metric2t = m
    dst = to_svm(x, y, w, group = group)
    run(`$RANKLIB -train $dst -save ranklib -ranker $ranker -metric2t $metric2t`)
    m.ranklib = read("ranklib")
    foreach(rm, ["ranklib", dst])
    BSON.bson("model.bson", model = m)
    return m
end

function predict(m::RanklibRanker, x)
    @unpack ranklib = m
    dst = to_svm(x)
    write("ranklib", ranklib)
    run(`$RANKLIB -load ranklib -rank $dst -score preds`)
    ŷ = readdlm("preds", Float32)[:, 3]
    foreach(rm, ["ranklib", "preds", dst])
    return ŷ
end


text:
using Base.Test

mesh = generateMesh(Mesh1D, 0, 5, 50)

ns = nodes(4)
𝜒 = atlas(mesh)
𝜒⁻¹ = invatlas(mesh)
k = 25
@test_approx_eq map(x->𝜒(k,𝜒⁻¹(k,x)),ns) ns

# 2D tests

for p,n in zip(AC274_2D.P1,AC274_2D.N1)
    @test fapply(p,n) == 1
end

for p,n in zip(AC274_2D.P2,AC274_2D.N2)
    @test fapply(p,n) == 1
end

for p = (1,2,3,4)
    @test AC274_2D.do_quad(x->1,p) == 0.5
end


text:
module CellDataTests

using Test

@testset "CellFields" begin include("CellFieldsTests.jl") end

@testset "CellQuadratures" begin include("CellQuadraturesTests.jl") end

@testset "QPointCellFields" begin include("QPointCellFieldsTests.jl") end

@testset "CellDofBases" begin include("CellDofBasesTests.jl") end

@testset "AttachDirichlet" begin include("AttachDirichletTests.jl") end

@testset "AttachConstraints" begin include("AttachConstraintsTests.jl") end

@testset "Law" begin include("LawTests.jl") end

end # module


text:
## This file auto-generated by Yggdrasil/build_hashes.jl
shard_hash_table = Dict{CompilerShard,String}(
    CompilerShard("Rootfs", v"2018.9.18", Linux(:x86_64, libc=:glibc), :targz)
      => "c1af8465a7269439e01ca7f6a94c45232bcf72604f1914fe753bcba327c9921e",
    CompilerShard("Rootfs", v"2018.9.18", Linux(:x86_64, libc=:glibc), :squashfs)
      => "ea97920cb26dedf55dc09195968e4ecfb2b5e5b01126ace761d30ab87b20bdb5",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "99d07a3b5c5beacf99225a97e9967595e10af1ba746a111e70ae2ad58d7bad58",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "a4c1da447504500ddfdddb549a9d501b038d63dc0eb51f9428a46cec78137c6a",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "c6b17620f51bf2858930e77def0ba94695a1859161a77973f68dd88fe65f026e",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "479a63174fed3fb1896935d195f987c2a7476cfd95e65461125895fa9784dc49",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "066ad44bdb4b64f4383332f60fd2d25cf8205a0fc4669b5383e90f65ee7c0be6",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "764433184b49aa91b3e31a19af545a552fa3f47c9f9958cd3391915e002021ea",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "fcc5316cd7b9705d7458220e39117f663ba5a8538415c0dab7dcf0a2bb683e38",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "a7f8c427c990e3f98cef865781cd79df06e80fe1b6133f1fbea9363a69940bfb",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "f1752354ec2c0d12208afb3d24d66f8aba2c81e95eab4256dff14e24785d48b5",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "64d38d6bc74c26518a10f6b79b4eff97277e27439976f7ddb08e71aeb4c68138",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "f49fb42525cfb39973074c25cb97e40815965bf771fb0c4b8428506272e09de0",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "79be7ed29f6146cc0a733801396361a7ce46e33411cf4926fa2b35d8ffdeee64",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "cc4c558d813ce04e754727945a86d982755aabfae6ef43b617f8643801d7887f",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "0a42123eba096786a6129c049f70076d095586ef9360cfd63d3f25a4709193c4",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "956e1b8822cb22e3d4d7dbf8d09c1ed24378925608ce23e62a82d1bbeadbbb03",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "2ea88f39229da11b6b1d26a35b2cbcefcf9d9d22ca17c7a159a74773822c7107",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "4c6971687a0337e639b51ffd5205b1914afb06e559bac4876b4ebefb5747bae7",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "4f8d7e1563143037ec526237c22687836bb4972cce324037d2fa625a4bfa2d87",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "8f454b7d7e13853ef0e6e371fe1f0f3a525fafd66527c46ce7fa15e5d615fc71",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "6b8635430269b4d30652b0f8cd14cabf53389587c462379b4b42b2a22dece5db",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "46aac87721fea5dd45b436133093bacce0e11dbeae3754fa19811fca6088a030",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "c50da72b2658888e465916ba85b8a68f9201ac972634a3f2b146305d310b7848",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "1532cacb9b5fda1fb880dd8ce6bcfb507380be966fe2a4bb089e02f57762a3f1",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "f9d2f07aa5fcd73fea43ac59407b871d96be2b87368a1c866ee9d2326f046785",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "36505a472f5c82facc555ce444825d9e2dcee923acbf7cbb49b321e6768bbca9",
    CompilerShard("BaseCompilerShard", v"2018.10.10", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "9acc139d4e561e259b22cb3e9299b4b8335d412ea72264aab40645a4d5db9d8d",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "7dc166b739b9aa587ce2d4cbf7fa020704489942ad4f2526e232ba30a6d3435d",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "4ad459a0375b4c6c923e979bf0e0c6f65f93b7bd43f6f4916639bbf99b4e8d51",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "5b55cc6ee306f93bfc4cb47f4a680e3ca84997e611b1d9d12c612f6d9eeaef70",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "0fc66a43c744dc24b8e8a0dbcb5143af6eb2d6eaa3780a48dfc32254ccc0d38d",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-gnu"))
      => "c806a846d89dd7801ee561e14d1169e2e0f90c461b70de6beaa01347d7ef8c21",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-gnu"))
      => "139c48a5b6bc2f32ce4db7612b23d5ea44878a0ab436521a54c6820c863d77c7",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "3d3ce2af22169fc4626b30c2faa87fcbca2c6a178051c7acdf3d5874b927ea9a",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "d02784b6cb34ebe6ff2049294a4012e56977f6cdda717bc860eea1b9a84b66cd",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "0eef33759de7915b18c62880f29da77ece8ac296d2ba4e4f6acd1df081c1c2d8",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "199a52fed6b56c20f3615bbe3f884a532803862d602eb9ed797d93e273159980",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("aarch64-linux-musl"))
      => "7009e1e1e46074e8415ae303cb177783e376e277e959bef2de8dda75768d254b",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("aarch64-linux-musl"))
      => "551398cc28d32eea9651b77055baae05ce387f8396e09761c34052d7b19c6849",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "12c156f32eec4a295b7005263258faad0ee78131a6fe96a86e16e8a97993174e",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "886e765c8da85da3df6ae85cea7ba2d0636818c22cc7e9dd727dccbb0fe92bdf",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "1a548f235068e991aae2dd6c48463d6ff4907939fe8805c32946e3407bd0324e",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "c3e7635405651f56151497415ee9d5dda93903e63d603bf2b67aae100691cb32",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-gnueabihf"))
      => "48de5284b5b49492bf4a93dd34875038b75cc6ed99a1967944af91d9b12b0d78",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-gnueabihf"))
      => "bb8f931df2253e617f71a612ff1f5c6be4b1f2fc6f2d546a2ea11b5f5198b571",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "7850cb66369c8f76946dea0261cbd201140bfdcb2baae455a41f204a03075f92",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "ae487494f2df74f534d442ffb8733ef45d8c2385f41150062fa4d18166449216",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "6657446b0614928a5c45d704864c6d4ff6e1695f96d6ecdae97eadd2b72cf995",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "7ab1c8af2e615e417f3cb6d2956a9190deb9305f0e29560ab58bae11b5765e01",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("arm-linux-musleabihf"))
      => "8780bad241eff84c17772277c083d26a47ab4590417467d978a2d289e7e3cd7f",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("arm-linux-musleabihf"))
      => "f92cdf344f688529c3d7b4979054a3a5a524c93dbc48246d59455628dcb20708",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "748cf8fdf995b9b8455f911cdd724d3a5eaee2186aec8ef3d796360559fb3ff6",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "196bbaea8ad69640ec84fc23c6a6dfc1a876ede81b0ada93ede6a81229fa8775",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "019612f298c42fc678dd8f2311ea4829e5a20f6bd064396de82729822e8e7ab4",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "d01027dcd2449408dca68a60333302ff84517a05d48d1b0547dd97a20852a4ac",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-gnu"))
      => "28774a18cc00c353ac0fb2d0041337de46426cdc8ce88bbdb969c38eb666e0b8",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-gnu"))
      => "dce0beb546bd4060e6737b0a5d997049f0f386c173708934b4de73cecd691251",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "a61f94db0252c7a45fad4d296bc24ab662e19425aded183fcc9d0f0ace55383d",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "cd1543254059ca984d77d597e097375e65318c76b55680d7118bdbeee5667300",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "10eeff7ae6adff390050a4cee9dd162ae574a93af3273207860fade15643d835",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "19f3255667baf9dc1ccbac7cf521af2a5f0d050e426fc1ad4fc640544f1b25f5",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-linux-musl"))
      => "c2d0fd3a594036ccdc5b67eaf9b9d150f7bcf7a76096208c38681a7cabe2960e",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-linux-musl"))
      => "d79987b21b85fddca127115df9ae7ff128b450e9408fa40b7f146a1b9e99916a",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "a6df16225fe350b14445e11c4e9b87f9ff813e5286cece8c86f24fa14afc3e24",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "5c2c0e51bac622b6c738ac6d9f6389ed67264e37e7e9b1bc5c9ef20fbbcf8841",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "8c3232dd0c1cc579342b391e9756c6f33891eb0bc3dd82eb1eb26c658824bfe1",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "711873d2999d3aa8fb9b61cd9ad0e5bee92de54e16928ecf0dde135e04d2f544",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("i686-w64-mingw32"))
      => "961f369d9d4d89ca328ed3ff072b744243fe2b7e63e0aed52a139ce1c64b9dab",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("i686-w64-mingw32"))
      => "9c696daf5a9411df1e9d3ed0b66a2d16f28de85a351a6eaf1127a7bff329efd5",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "dab6635da13f65823882b0eba7c07922993977d1e8860f7cb30e31028ff0a37f",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "f339fb925614b1525dc347127af57f933c83d9428b1211fe4a84a62a1752fb80",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "985f3161eecdb4d9c9822292e126aaeab4ef4fd723f13b9774025cb33407fdd2",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "f5914829c645060358a914c410b4aa68c1cbf6110cd55d74f93c6378d3d9ec9b",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "e0b9d8d844385501583d660f8c63b34cd2ba1615696301b136a68c59330a80c7",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("powerpc64le-linux-gnu"))
      => "3ac831c81077e19ee62dbef57efd4eb9e1d36489e2f0b2a602578a4e1fd845ef",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "988acdd754238b732c09ca6e401efa44f21ceacbdccfb8945df6f25d2be3e024",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "7bfbf57cd1053fe502165602024ac7218d7112bcb416ce7113f5783e371416d2",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "fed1e6da9c2e1726bcd2943482343bdd3e7374f3cb3580d42a87e47061bef69a",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "7fc25bdc08a94e2c05f065436fbe23612dff37b2d2aaa276348afe2b52bc6aca",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-apple-darwin14"))
      => "eca90a29631f07e64f8fcc8ebbf72f6f99518c95d449173ee3db43ff7f494fa2",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-apple-darwin14"))
      => "878deaa5632d4954f638c90c544855c64b7b1d409541568361e2da5b3f4a7a8c",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "4566b42fc012f7c409c7b84c228448181554eb2a9c6ab4f6e11134d5020f2a1e",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "a6227eac0630907fd2c5f809276eed9836a75e3d32fbd55b9475b266ab89cd49",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "e7616725db345c8cb55273b8273a71106945b1fa5ded9b3111f8bdcc9096efb5",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "173770508429468b54b7c10bbb3903ca244c7842956f5b1dbbe88a476451c2ea",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-gnu"))
      => "de416b08ee85a08c0c0de32dab9c77946d536b2b30a4bc9bfe96acfac0424f00",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-gnu"))
      => "e9584274e172c1dd96fbc20936bd303de12a7061b091315795465eac1709dffe",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "cfd8b9b5ece88bb4b31ccec23f842dcf1b42baa793bb5578aac7e3dee3dadd76",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "d8352a9097c7154caec3198d4d331381f2732d292a544d83d995d8c443ac34bf",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "2c599807efb6a61d2fd03077fd51ebc439edac7dabc813a547cdef44423119cc",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "6db08582a959772ae6593305407bda521b527e8cd4c6ea4843bb61b17ca28a54",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-linux-musl"))
      => "199a430d2165a0c6cb6835399300b5ec4bba9e85a2e1d6ed4ed5c7167af6ecc7",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-linux-musl"))
      => "7fb291b9b5c9761a5f478f558a20d8e1af1734435d321e5e116d956a46094577",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "3cf20599c1ed19bd80148e36dfe241b808b35cdbe85db251de24211b5d5dc2ff",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "6755b76f31531352f29e0859aef6426ed0b529b2d466b629cdbf205abc100f87",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "dc33e79115108bd8474aa84a6879e54f67d2ab3798a18e23c0c73afc236952b3",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "b12a4f982b80d619afcf7680e8f37d3d546f7b4d97e367abe3503f81f1ff5646",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "f85ce7085f811cc7e1fe64104e448e5c01f6f388bbf92af3d4af107c641abc8e",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-unknown-freebsd11.1"))
      => "cb7f44bd30968e86fde8db36ae0e70154180d24f43c187cc73b4d3b44e581bfb",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "9c3f9e060d293746b5124197e91dbbd561b5b72f684a1c3fdf308c0e001166fb",
    CompilerShard("GCC", v"4.8.5", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "4d3781a7769923f5d67c4ccdc99497d7baef5216bc370423395fb508cd6ed19c",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "779bfcf40c920e3cd744d20684af4cfbf9e553a0eb5522bc4451b1e3dd64571c",
    CompilerShard("GCC", v"7.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "283dd88174df8be0c20b35aeb4ca7ee1fc143dbc239100b40cd0146075ada4e8",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :targz; target=platform_key_abi("x86_64-w64-mingw32"))
      => "5e9bc0b0970f3852f33e3360ed658a9b538b5b09e6b87d1f0dd0293402c17761",
    CompilerShard("GCC", v"8.1.0", Linux(:x86_64, libc=:glibc), :squashfs; target=platform_key_abi("x86_64-w64-mingw32"))
      => "b24ba4cb2a3749fe1f3d6503533fe9f504a24baf9f16ba37ecbb1fc7c44f05db",
    CompilerShard("LLVM", v"6.0.1-0", Linux(:x86_64, libc=:glibc), :targz)
      => "d53a63d79c29b7a3d0eb3e516e9b8b08cb6d2dfb692a67242323b3737d7e0e9a",
    CompilerShard("LLVM", v"6.0.1-0", Linux(:x86_64, libc=:glibc), :squashfs)
      => "80578c2f73b6d3cd4a5138bd732f4b2aff28968e0f2cf1b35e290f2a202de9f8",
)


text:
@testset "Ewald Summation" begin

function test_ewald_NaCl()
    lat_vecs = [1.0 0   0;
                0   1.0 0;
                0   0   1.0]
    b_vecs = [zeros(3)]
    latsize = [2, 2, 2]
    lattice = Sunny.Lattice(lat_vecs, b_vecs, latsize)
    sys = ChargeSystem(lattice)
    sys.sites .= reshape([1, -1, -1, 1, -1, 1, 1, -1], (1, 2, 2, 2))

    ewald_result = Sunny.ewald_sum_monopole(lattice, sys.sites, extent=30)

    answer = -1.7475645946331822
    @test isapprox(answer, ewald_result / 4; rtol=1e-7)
end

test_ewald_NaCl()

function test_ewald_CsCl()
    lat_vecs = [1.0 0   0;
                0   1.0 0;
                0   0   1.0]
    b_vecs = [[0.,  0.,  0.],
              [0.5, 0.5, 0.5]]
    latsize = [1, 1, 1]
    lattice = Sunny.Lattice(lat_vecs, b_vecs, latsize)
    sys = ChargeSystem(lattice)
    sys.sites .= reshape([1, -1], (2, 1, 1, 1))

    ewald_result = Sunny.ewald_sum_monopole(lattice, sys.sites, extent=30)

    # Madelung constants are reported relative to the
    #  nearest-neighbor distance in the crystal
    ewald_result *= √(3/4)

    answer = -1.76267477307099
    @test isapprox(answer, ewald_result; rtol=1e-7)
end

test_ewald_CsCl()

function test_ewald_ZnS()
    lat_vecs = [0.0 0.5 0.5;
                  0.5 0.0 0.5;
                  0.5 0.5 0.0]
    b_vecs = [[0.,  0.,  0.],
              [0.25, 0.25, 0.25]]
    latsize = [1, 1, 1]
    lattice = Sunny.Lattice(lat_vecs, b_vecs, latsize)
    sys = ChargeSystem(lattice)
    sys.sites .= reshape([1, -1], (2, 1, 1, 1))

    ewald_result = Sunny.ewald_sum_monopole(lattice, sys.sites, extent=30)

    # Madelung constants are reported relative to the
    #  nearest-neighbor distance in the crystal
    ewald_result *= √(3/16)

    answer = -1.63805505338879
    @test isapprox(answer, ewald_result; rtol=1e-7)
end

test_ewald_ZnS()

function test_ewald_ZnSB4()
    a = 1.
    c = √(8/3) * a
    u = 3/8

    lat_vecs = [ 0.5a    0.5a    0.0;
                -0.5*√3a 0.5*√3a 0.0;
                 0.0     0.0       c]
    b_vecs = [[0.5a,  0.5/√3*a,  0.],
              [0.5a, -0.5/√3*a, 0.5c,],
              [0.5a, 0.5/√3*a, u*c],
              [0.5a, -0.5/√3*a, (0.5+u)*c]]
    b_vecs = [[1/3, 2/3, 0],
              [2/3, 1/3, 0.5],
              [1/3, 2/3, 3/8],
              [2/3, 1/3, 7/8]]
    latsize = [1, 1, 1]
    lattice = Sunny.Lattice(lat_vecs, b_vecs, latsize)
    sys = ChargeSystem(lattice)
    sys.sites .= reshape([1, 1, -1, -1], (4, 1, 1, 1))

    ewald_result = Sunny.ewald_sum_monopole(lattice, sys.sites, extent=30)

    # Madelung constants are reported relative to the
    #  nearest-neighbor distance in the crystal
    ewald_result *= u * c

    answer = -1.64132162737
    @test isapprox(answer, ewald_result / 2; rtol=1e-7)
end

test_ewald_ZnSB4()

"Self-energy of a physical dipole with moment p, and displacement d=2ϵ"
function _dipole_self_energy(; p::Float64=1.0, ϵ::Float64=0.1)
    d, q = 2ϵ, p/2ϵ
    return -q^2 / d
end

"""Approximates a dipolar `SpinSystem` by generating a monopolar `ChargeSystem` consisting of
    opposite charges Q = ±1/(2ϵ) separated by displacements d = 2ϵp centered on the original
    lattice sites.
"""
function _approx_dip_as_mono(sys::SpinSystem; ϵ::Float64=0.1) :: ChargeSystem
    lattice = sys.lattice
    sites = sys.sites

    # Need to expand the underlying unit cell to the entire system size
    new_lat_vecs = lattice.size' .* lattice.lat_vecs
    new_latsize = [1, 1, 1]

    frac_transform = inv(new_lat_vecs)

    new_nbasis = 2 * prod(size(sites))
    new_sites = zeros(new_nbasis, 1, 1, 1)
    new_basis = empty(lattice.basis_vecs)
    sizehint!(new_basis, new_nbasis)

    ib = 1
    for idx in eachindex(lattice)
        @inbounds r = lattice[idx]
        @inbounds p = sites[idx]

        # Add new charges as additional basis vectors
        push!(new_basis, frac_transform * (r + ϵ * p))
        push!(new_basis, frac_transform * (r - ϵ * p))

        # Set these charges to ±1/2ϵ
        new_sites[ib, 1, 1, 1]   =  1 / (2ϵ)
        new_sites[ib+1, 1, 1, 1] = -1 / (2ϵ)

        ib += 2
    end

    new_basis = map(b->mod.(b, 1), new_basis)

    new_lattice = Sunny.Lattice(new_lat_vecs, new_basis, new_latsize)

    return ChargeSystem(new_lattice, new_sites)
end

"""
Tests that `ewald_sum_monopole` and `ewald_sum_dipole` give consistent results
 up to approximation error.
TODO: This version uses existing functions, but suffers from catastrophic
 cancellation issues due to the self-energy exploding. Should write a
 separate monopole ewald which skips these self-energy interactions.
"""
function test_mono_dip_consistent()
    lat_vecs = [1.0 0   0;
                0   1.0 0;
                0   0   1.0]
    b_vecs = [[0.,  0.,  0.],
              [0.5, 0.5, 0.5]]
    latsize = [2, 2, 2]
    cryst = Crystal(lat_vecs, b_vecs)
    sys = SpinSystem(cryst, Sunny.AbstractInteraction[], latsize)
    rand!(sys)

    dip_ewald = Sunny.ewald_sum_dipole(sys.lattice, sys.sites; extent=15)

    csys = _approx_dip_as_mono(sys; ϵ=0.001)
    mono_ewald = Sunny.ewald_sum_monopole(csys.lattice, csys.sites; extent=15)
    dip_self_en = _dipole_self_energy(; ϵ=0.001)
    num_spins = prod(size(sys))

    @test isapprox(dip_ewald, mono_ewald - num_spins * dip_self_en; rtol=1e-3)
end

test_mono_dip_consistent()

end

# Beck claims this should be independent of η, and converge to -2.837297
# I find that it only does for large η?
function test_ξ_sum(;extent=2, η=1.0) :: Float64
    extent_ixs = CartesianIndices(ntuple(_->-extent:extent, Val(3)))

    real_space_sum = 0.0
    recip_space_sum = 0.0

    for JKL in extent_ixs
        JKL = convert(SVector, JKL)
        k = 2π * JKL

        if all(JKL .== 0)
            continue
        end

        dist = norm(JKL)
        kdist = 2π * dist

        real_space_sum += erfc(η * dist) / dist
        recip_space_sum += exp(-kdist^2 / (4η^2)) / kdist^2
    end

    return real_space_sum + 4π * recip_space_sum - 2η/√π - π/η^2
end


text:
@testset "646.maximum-length-of-pair-chain.jl" begin
    @test find_longest_chain([(1, 2), (2, 3), (3, 4), (1, 3)]) == 2
    @test find_longest_chain([(1, 2), (1, 2), (1, 2)]) == 1
    @test find_longest_chain([(1, 2), (2, 3)]) == 1
    @test find_longest_chain([(1, 2), (2, 3), (3, 4)]) == 2
    @test find_longest_chain(Tuple{Int, Int}[]) == 0
end


text:
#=
surface_test:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-19
=#
using Distributed
if nprocs() == 1
    addprocs(3)
end
# # using Revise
using Test
using Logging
using SparseArrays
# ENV["JULIA_DEBUG"] = "surface_extraction_parallel"
# ENV["JULIA_DEBUG"] = "all"
# using ViewerGL
# using Plasm
@everywhere using LarSurf
# @everywhere using Distributed
# global_logger(SimpleLogger(stdout, Logging.Debug))
# # set logger on all workers
# for wid in workers()
#     @spawnat wid global_logger(SimpleLogger(stdout, Logging.Debug))
# end
@testset "Init and deinit" begin
    block_size = [2, 2, 2]
    # println("Channel type apriori $(typeof(LarSurf._ch_block))")
    LarSurf.lsp_setup(block_size)
    # println("Channel type aposteriori: $(typeof(LarSurf._ch_block))")
    # println("channel aposteriori $(LarSurf._ch_block)")
    LarSurf.lsp_deinit_workers()
end

@testset "Setup parallel surface computation" begin
    block_size = [2, 2, 2]
    LarSurf.lsp_setup(block_size)
    for wid in workers()
        # println("testing on $wid")
        ftr = @spawnat wid LarSurf._single_boundary3
        @test fetch(ftr) != nothing
    end
end

#
@testset "parallel surface extraction " begin
    block_size = [2, 2, 2]
    segmentation = LarSurf.data234()

    LarSurf.lsp_setup(block_size)

    @debug "Setup done"
    larmodel = LarSurf.lsp_get_surface(segmentation)
    @test LarSurf.check_surface_euler(larmodel[2])
    FVtri = LarSurf.triangulate_quads(larmodel[2])
    # @test LarSurf.check_surface_euler(FVtri)
    # LarSurf.check_LARmodel(larmodel)
    # Plasm.view()
    # Plasm.view( Plasm.numbering(.6)(larmodel) )
	V, FV = larmodel
	# FV1 = [[f[1], f[3], f[4]] for f in FV]
	# FV2 = [[f[1], f[3], f[2]] for f in FV]
	# FVtri = vcat(FV1, FV2)
# ViewerGL.VIEW([
#     ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.1))
# 	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
# ])
# ViewerGL.VIEW([
#     ViewerGL.GLGrid(V,FV,ViewerGL.Point4d(1,1,1,0.1))
# 	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
# ])

end

@testset "run on local" begin

    block_size = [2, 2, 2]
    segmentation = LarSurf.data234()
    block_id = 3

    data_size = LarSurf.size_as_array(size(segmentation))

    LarSurf.lsp_setup(block_size)
    # b3, larmodel = LarSurf.get_boundary3(block_size)
    # LarSurf.set_single_boundary3(b3, block_size)
    n, bgetter = LarSurf.block_getter(segmentation, block_size;fixed_block_size=true)
    block = LarSurf.get_block(block_id, bgetter...)
    # outdata, offset, sz =block
    data_for_channel = (block..., block_id)

    # tm_put = @elapsed put!(LarSurf._ch_block, data_for_channel)
    # fbl = take!(LarSurf._ch_block)
    # alternativelly can be channel step skipped with:
    fbl = data_for_channel
    faces = LarSurf.code_multiply_decode(data_size, fbl...)
    @test length(faces) == 16

end
#

# @testset "parallel surface extraction big" begin
#     block_size = [64, 64, 64]
#     # segmentation = LarSurf.data234()
#     segmentation = LarSurf.generate_cube(512)
#
#     LarSurf.lsp_setup(block_size)
#     # for wid in workers()
#     #     # println("testing on $wid")
#     #     ftr = @spawnat wid LarSurf._single_boundary3
#     #     @test fetch(ftr) != nothing
#     # end
#
#     @debug "Setup done"
#     tt = @elapsed LarSurf.lsp_get_surface(segmentation)
#     println("time: $tt")
#
#
# end
# ch = RemoteChannel(()->Channel{Int}(32));

# @testset "Job enquing parallel surface computation" begin
#     block_size = [8, 8, 8]
#     LarSurf.lsp_job_enquing(block_size)
#

    # for wid in workers()
    #     println("testing on $wid")
    #     ftr = @spawnat wid LarSurf._single_boundary3
    #     @test fetch(ftr) != nothing
    # end

# end


text:
struct OraExecMode
    val::UInt32
end

UInt32(exec_mode::OraExecMode) = exec_mode.val

Base.convert(::Type{OraExecMode}, val::UInt32) = OraExecMode(val)
Base.convert(::Type{UInt32}, exec_mode::OraExecMode) = exec_mode.val

OraExecMode(val::Integer) = OraExecMode(UInt32(val))

const ORA_MODE_EXEC_DEFAULT             = OraExecMode(0)
const ORA_MODE_EXEC_DESCRIBE_ONLY       = OraExecMode(16)
const ORA_MODE_EXEC_COMMIT_ON_SUCCESS   = OraExecMode(32)
const ORA_MODE_EXEC_BATCH_ERRORS        = OraExecMode(128)
const ORA_MODE_EXEC_PARSE_ONLY          = OraExecMode(256)
const ORA_MODE_EXEC_ARRAY_DML_ROWCOUNTS = OraExecMode(1048576)

# `OraExecMode` can be "ored" together
Base.:|(a::OraExecMode, b::OraExecMode) = OraExecMode(a.val | b.val)
Base.:&(a::OraExecMode, b::OraExecMode) = OraExecMode(a.val & b.val)


text:
using MatrixNetworks

function partition_and_map_subgraph(A,subset,minsize)
  map = collect(subset)
  A1 = A[map,map]
  sets = recursive_spectral(A1,minsize)
  for s in sets
    for i = 1:length(s)
      s[i] = map[s[i]]
    end
  end
  return sets
end

function recursive_spectral(A,minsize)
  n = size(A,1)
  if n < minsize
    sets = Vector{Vector{Int}}()
    push!(sets,Vector{Int}(1:size(A,1)))
    return sets
  else
    sc = spectral_cut(A)
    s1 = sc.set
    s2 = setdiff(1:n, s1)
    println("partitioning $n nodes into $(length(s1)) and $(length(s2)),lam $(sc.lam2) cond $(minimum(sc.sweepcut_profile.conductance))\n")
    # this is always in the largest component
    sets = partition_and_map_subgraph(A,s1,minsize)
    s2 = partition_and_map_subgraph(A,s2,minsize)
    for s in s2
      push!(sets,s)
    end
    return sets
  end
end

function set2clusvec(sets,n)
  c = zeros(n)
  next = 1
  for s in sets
    c[s] .= next
    next += 1
  end
  return c
end
## This is just a test
using MAT
mat = matread("../data/emailEUcore.mat")
A = mat["A"]            # Largest connected component of email EU
ASim = mat["Asim"]      # Undirected version of same graph: ASim = spones(A+A')

sets = recursive_spectral(ASim,100)

n = size(A,1)
c = set2clusvec(sets,n)

@show maximum(c)

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
module GMT

using Printf

# Need to know what GMT version is available or if none at all to warn users on how to
# install GMT.
try
	# Due to a likely Julia bug next command fails when this file called with 'using'
	#ver_s = @capture_out run(`gmt --version`);
	#@show(length(ver_s))			# Prints 0 length

	# So resort to write result to file on disk and read from it
	write("gmtversion__.txt", read(`gmt --version`))
	global const GMTver = Meta.parse(read("gmtversion__.txt", String)[1:3])
	rm("gmtversion__.txt")
	global foundGMT = true
catch
	global foundGMT = false
	global const GMTver = 5.0		# Don't want to raise an error in libgmt_h.jl due to a missing var
end

#const GMTver = 6.0
const FMT = "ps"

export
	GMTver, FMT, gmt,
	bar, bar!,
	basemap, basemap!, blockmean, blockmedian, blockmode, coast, coast!, colorbar, colorbar!, contour,
	contour!, filter1d, filter2d, fitcircle, gmt2kml,  gmtconnect, gmtconvert, gmtinfo, gmtregress,
	gmtread, gmtselect, gmtsimplify, gmtspatial, gmtvector, gmtwrite, gmtwich,
	grd2cpt, grd2kml, grd2xyz, grdblend, grdclip, grdcontour,
	grdcontour!, grdcut, grdedit, grdfft, grdfilter, grdgradient, grdhisteq, grdimage, grdimage!, grdinfo,
	grdlandmask, grdpaste, grdproject, grdsample, grdtrack, grdtrend, grdview, grdview!, grdvolume,
	grid_type, histogram, histogram!, image, image!, imshow, imshow!, kml2gmt, logo, logo!, makecpt,
	mapproject, nearneighbor, plot, plot!, plot3d, plot3d!, project,
	pscontour, pscontour!, psconvert, psbasemap, psbasemap!, psclip, psclip!, pscoast, pscoast!,
	pshistogram, pshistogram!, psimage, psimage!, psmask, psmask!, psrose, psrose!, psscale, psscale!,
	pssolar, pssolar!, psternary, psternary!, pstext, pstext!, pswiggle, pswiggle!, psxy, psxy!, psxyz,
	psxyz!, rose, rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, spectrum1d,
	sphdistance, sphinterpolate, sphtriangulate,
	surface, ternary, ternary!, text, text!, text_record, trend1d, trend2d, triangulate, splitxyz,
	wiggle, wiggle!, xy, xy!, xyz2grd

include("common_docs.jl")
include("libgmt_h.jl")
include("libgmt.jl")
include("gmt_main.jl")
include("common_options.jl")
include("blocks.jl")
include("filter1d.jl")
include("fitcircle.jl")
include("gmt2kml.jl")
include("gmtconnect.jl")
include("gmtconvert.jl")
include("gmtinfo.jl")
include("gmtlogo.jl")
include("gmtreadwrite.jl")
include("gmtset.jl")
include("gmtselect.jl")
include("gmtsimplify.jl")
include("gmtspatial.jl")
include("gmtregress.jl")
include("gmtvector.jl")
include("grd2cpt.jl")
include("grd2kml.jl")
include("grd2xyz.jl")
include("grdblend.jl")
include("grdclip.jl")
include("grdcut.jl")
include("grdedit.jl")
include("grdcontour.jl")
include("grdfft.jl")
include("grdfilter.jl")
include("grdhisteq.jl")
include("grdinfo.jl")
include("grdimage.jl")
include("grdgradient.jl")
include("grdlandmask.jl")
include("grdpaste.jl")
include("grdproject.jl")
include("grdsample.jl")
include("grdtrack.jl")
include("grdtrend.jl")
include("grdview.jl")
include("grdvolume.jl")
include("imshow.jl")
include("kml2gmt.jl")
include("makecpt.jl")
include("mapproject.jl")
include("nearneighbor.jl")
include("project.jl")
include("psbasemap.jl")
include("pscoast.jl")
include("pscontour.jl")
include("psconvert.jl")
include("pshistogram.jl")
include("psimage.jl")
include("psscale.jl")
include("psrose.jl")
include("pssolar.jl")
include("pstext.jl")
include("psxy.jl")
include("pswiggle.jl")
include("plot.jl")
include("plot3d.jl")
include("sample1d.jl")
include("spectrum1d.jl")
include("sphdistance.jl")
include("sphinterpolate.jl")
include("sphtriangulate.jl")
include("splitxyz.jl")
include("surface.jl")
include("triangulate.jl")
include("trend1d.jl")
include("trend2d.jl")
include("xyz2grd.jl")

if (!foundGMT)
	println("\n\nYou don't seem to have GMT installed and I don't currently install it automatically,")
	println("so you will have to do it yourself. Please follow instructions bellow but please note")
	println("that since GMT is migrating to Github (https://github.com/GenericMappingTools/gmt),")
	println("some of the links may change in a near future.\n\n")
	if (Sys.iswindows() && Sys.WORD_SIZE == 64)
		println("1) Download and install the official version at (the '..._win64.exe':")
		println("\t\t https://gmt.soest.hawaii.edu/projects/gmt/wiki/Download")
		println("\n2) Or even better, download and install the GMT6dev version at:")
		println("\t\t http://w3.ualg.pt/~jluis/downloads/gmt.html")
	elseif (Sys.iswindows() && Sys.WORD_SIZE == 32)
		println("Download and install the official version at (the '..._win32.exe':")
		println("\t\t https://gmt.soest.hawaii.edu/projects/gmt/wiki/Download")
	else
		println("https://gmt.soest.hawaii.edu/projects/gmt/wiki/BuildingGMT")
	end
end

end # module


text: run(`python3 mand.py`)