text:
module Toolips
#==
Toolips.jl, a module created for light-weight interaction with Javascript via
bursts of HTML and Javascript working in tandem. There might be entirely new
intentions for this design, and dramatic changes. There are numerous features
coming to spruce up this code base to better facilitate things like
authentication.
~ TODO LIST ~ If you want to help out, you can try implementing the following:
=========================
- TODO Incoming user information would be nice.
- TODO Secret keys.
- TODO Authentication. (UUID's would be nice)
- TODO Parsing envvariables and CLI from main in the generated file (this file,
create_serverdeps(name::String)
- TODO Production vs dev environments
==#
using Crayons
using Sockets, HTTP, Pkg
include("interface/servables.jl")
# Server
export Route, ServerTemplate, Logger, stop!
# Components
export Page
export html, html_file, file, css, css_file, js, js_file, fn
export Button, Form, TextArea, TextBox, RadioSet, Slider
export Canvas
export List, UnorderedList, A, DropDown
# Structure servables (frontend)
export Header, Navbar, Body, Columns
# Methods
export getargs, getarg, getpost, write_file, lists

function create_serverdeps(name::String)
    Pkg.generate(name)
    dir = pwd() * "/"
    src = dir * name * "/src"
    public = dir * name * "/public"
    logs = dir * name * "/logs"
    mkdir(public)
    mkdir(logs)
    touch(name * "/dev.jl")
    touch(name * "/prod.jl")
    touch(logs * "/log.txt")
    rm(src * "/$name.jl")
    touch(src * "/$name.jl")
    open(src * "/$name.jl", "w") do io
        write(io, """
# Welcome to your new Toolips server!
using Main.Toolips\n
function main()
        # Essentials
    global LOGGER = Logger()
    routes = make_routes()
    server_template = ServerTemplate(IP, PORT, routes)
    global TLSERVER = server_template.start()
    return(TLSERVER)
end
\n
# Routes
function make_routes()
        # Pages
    four04 = html("
404, Page not found!
")
    index = html("
Hello world!
")
        # Routes
    routes = []
    homeroute = Route("/", index)
    four04route = Route("404", four04)
    push!(routes, homeroute)
    push!(routes, four04route)
    routes
end
\n
main()
        """)
    end
    open(name * "/dev.jl", "w") do io
        write(io, """
        IP = "127.0.0.1"
        PORT = 8000
        PUBLIC = "../public"
        include("src/$name.jl")
        """)
    end
    open(name * "/prod.jl", "w") do io
        write(io, """
        IP = "127.0.0.1"
        PORT = 8000
        PUBLIC = "../public"
        include("src/$name.jl")
        """)
    end
end
function new_webapp(name::String = "ToolipsApp")
    create_serverdeps(name)
end

export new_webapp
# --

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Xorg_libxcb"
version = v"1.13"

# Collection of sources required to build libxcb
sources = [
    ArchiveSource("https://www.x.org/archive/individual/xcb/libxcb-$(version.major).$(version.minor).tar.bz2",
                  "188c8752193c50ff2dbe89db4554c63df2e26a2e47b0fa415a70918b5b851daa"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libxcb-*/
CPPFLAGS="-I${prefix}/include"
# When compiling for things like ppc64le, we need newer `config.sub` files
update_configure_scripts
./configure --prefix=${prefix} --host=${target} --enable-malloc0returnsnull=no
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p) || Sys.isfreebsd(p)]

products = [
    LibraryProduct("libxcb-composite", :libxcb_composite),
    LibraryProduct("libxcb-damage", :libxcb_damage),
    LibraryProduct("libxcb-dpms", :libxcb_dpms),
    LibraryProduct("libxcb-dri2", :libxcb_dri2),
    LibraryProduct("libxcb-dri3", :libxcb_dri3),
    LibraryProduct("libxcb-glx", :libxcb_glx),
    LibraryProduct("libxcb-present", :libxcb_present),
    LibraryProduct("libxcb-randr", :libxcb_randr),
    LibraryProduct("libxcb-record", :libxcb_record),
    LibraryProduct("libxcb-render", :libxcb_render),
    LibraryProduct("libxcb-res", :libxcb_res),
    LibraryProduct("libxcb-screensaver", :libxcb_screensaver),
    LibraryProduct("libxcb-shape", :libxcb_shape),
    LibraryProduct("libxcb-shm", :libxcb_shm),
    LibraryProduct("libxcb", :libxcb),
    LibraryProduct("libxcb-sync", :libxcb_sync),
    LibraryProduct("libxcb-xf86dri", :libxcb_xf86dri),
    LibraryProduct("libxcb-xfixes", :libxcb_xfixes),
    LibraryProduct("libxcb-xinerama", :libxcb_xinerama),
    LibraryProduct("libxcb-xinput", :libxcb_xinput),
    LibraryProduct("libxcb-xkb", :libxcb_xkb),
    LibraryProduct("libxcb-xtest", :libxcb_xtest),
    LibraryProduct("libxcb-xvmc", :libxcb_xvmc),
    LibraryProduct("libxcb-xv", :libxcb_xv),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_util_macros_jll"),
    BuildDependency("Xorg_xcb_proto_jll"),
    Dependency("XSLT_jll"),
    Dependency("Xorg_libXau_jll"),
    Dependency("Xorg_libXdmcp_jll"),
    Dependency("Xorg_libpthread_stubs_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# This example demonstrates basic manipulation of a visualized
# mechanism in a standalone window.

#  For a similar demo rendered in Jupyter instead, see the
#  `interactive_manipulation.ipynb` notebook.

# Import our packages
using MeshCatMechanisms
using RigidBodyDynamics

# Blink provides the standalone window we'll use to display
# the visualizer and controls
using Blink
# Install the Blink.AtomShell if it isn't already installed
AtomShell.isinstalled() || AtomShell.install()

# Create a random mechanism and its associated visualizer

# seed the random number generator so we get repeatable results
import Random
Random.seed!(75)
mechanism = rand_chain_mechanism(Float64,
    [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1:5]]...)
mvis = MechanismVisualizer(mechanism, Skeleton(randomize_colors=true))

# Open the visualizer in a new window
#
# We don't open windows when we're running this test
# on the Travis CI build servers
if !haskey(ENV, "CI")
    open(mvis, Window())
end

# Create sliders to manipulate the visualizer's configuration
widget = manipulate!(mvis)

# Render those sliders in a new window
#
# We don't open windows when we're running this test
# on the Travis CI build servers
if !haskey(ENV, "CI")
    body!(Window(), widget)
end


text:
include("external/spglib.jl")

"""Bring ``k``-point coordinates into the range [-0.5, 0.5)"""
function normalize_kpoint_coordinate(x::Real)
    x = x - round(Int, x, RoundNearestTiesUp)
    @assert -0.5 ≤ x < 0.5
    x
end
normalize_kpoint_coordinate(k::AbstractVector) = normalize_kpoint_coordinate.(k)

"""Construct the coordinates of the ``k``-points in a (shifted) Monkorst-Pack grid"""
function kgrid_monkhorst_pack(kgrid_size; kshift=[0, 0, 0])
    kgrid_size = Vec3{Int}(kgrid_size)
    kshift     = Vec3{Rational{Int}}(kshift)

    start   = -floor.(Int, (kgrid_size .- 1) ./ 2)
    stop    = ceil.(Int, (kgrid_size .- 1) ./ 2)
    kcoords = [(kshift .+ Vec3([i, j, k])) .// kgrid_size
               for i=start[1]:stop[1], j=start[2]:stop[2], k=start[3]:stop[3]]
    vec(normalize_kpoint_coordinate.(kcoords))
end

@doc raw"""
    bzmesh_uniform(kgrid_size; kshift=[0, 0, 0])

Construct a (shifted) uniform Brillouin zone mesh for sampling the ``k``-points.
Returns all ``k``-point coordinates, appropriate weights and the identity SymOp.
"""
function bzmesh_uniform(kgrid_size; kshift=[0, 0, 0])
    kcoords = kgrid_monkhorst_pack(kgrid_size; kshift=kshift)
    kcoords, ones(length(kcoords)) ./ length(kcoords), [one(SymOp)]
end

@doc raw"""
     bzmesh_ir_wedge(kgrid_size, symmetries; kshift=[0, 0, 0])

Construct the irreducible wedge of a uniform Brillouin zone mesh for sampling ``k``-points,
given the crystal symmetries `symmetries`. Returns the list of irreducible ``k``-point
(fractional) coordinates, the associated weights and the new `symmetries` compatible with
the grid.
"""
function bzmesh_ir_wedge(kgrid_size, symmetries; kshift=[0, 0, 0])
    all(isequal.(kgrid_size, 1)) && return bzmesh_uniform(kgrid_size; kshift)

    # Transform kshift to the convention used in spglib:
    #    If is_shift is set (i.e. integer 1), then a shift of 0.5 is performed,
    #    else no shift is performed along an axis.
    kshift = Vec3{Rational{Int}}(kshift)
    is_shift = map(kshift) do ks
        ks in (0, 1//2) || error("Only kshifts of 0 or 1//2 implemented.")
        convert(Int, 2 * ks)
    end

    # Filter those symmetry operations that preserve the MP grid
    kcoords_mp = kgrid_monkhorst_pack(kgrid_size; kshift)
    symmetries = symmetries_preserving_kgrid(symmetries, kcoords_mp)

    # Give the remaining symmetries to spglib to compute an irreducible k-point mesh
    # TODO implement time-reversal symmetry and turn the flag to true
    Ws = [symop.W for symop in symmetries]
    _, mapping, grid = spglib_get_stabilized_reciprocal_mesh(
        kgrid_size, Ws; is_shift, is_time_reversal=false
    )
    # Convert irreducible k-points to DFTK conventions
    kgrid_size = Vec3{Int}(kgrid_size)
    kirreds = [(kshift .+ grid[ik + 1]) .// kgrid_size for ik in unique(mapping)]
    kirreds = normalize_kpoint_coordinate.(kirreds)

    # Find the indices of the corresponding reducible k-points in `grid`, which one of the
    # irreducible k-points in `kirreds` generates.
    k_all_reducible = [findall(isequal(elem), mapping) for elem in unique(mapping)]

    # Number of reducible k-points represented by the irreducible k-point `kirreds[ik]`
    n_equivalent_k = length.(k_all_reducible)
    @assert sum(n_equivalent_k) == prod(kgrid_size)
    kweights = n_equivalent_k / sum(n_equivalent_k)

    # This loop checks for reducible k-points, which could not be mapped to any irreducible
    # k-point yet even though spglib claims this can be done.
    # This happens because spglib actually fails for some non-ideal lattices, resulting
    # in *wrong results* being returned. See the discussion in
    # https://github.com/spglib/spglib/issues/101
    for (iks_reducible, k) in zip(k_all_reducible, kirreds), ikred in iks_reducible
        kred = (kshift .+ grid[ikred]) .// kgrid_size

        found_mapping = any(symmetries) do symop
            # If the difference between kred and W' * k == W^{-1} * k
            # is only integer in fractional reciprocal-space coordinates, then
            # kred and S' * k are equivalent k-points
            all(isinteger, kred - (symop.S * k))
        end

        if !found_mapping
            error("The reducible k-point $kred could not be generated from " *
                  "the irreducible kpoints. This points to a bug in spglib.")
        end
    end

    kirreds, kweights, symmetries
end

# TODO Maybe maximal spacing is actually a better name as the kpoints are spaced
#      at most that far apart
@doc raw"""
Selects a kgrid size to ensure a minimal spacing (in inverse Bohrs) between kpoints.
A reasonable spacing is `0.13` inverse Bohrs (around ``2π * 0.04 \AA^{-1}``).
"""
function kgrid_from_minimal_spacing(lattice, spacing)
    lattice       = austrip.(lattice)
    spacing       = austrip(spacing)
    recip_lattice = compute_recip_lattice(lattice)
    @assert spacing > 0
    isinf(spacing) && return [1, 1, 1]

    [max(1, ceil(Int, norm(vec) / spacing)) for vec in eachcol(recip_lattice)]
end
function kgrid_from_minimal_spacing(model::Model, args...)
    kgrid_from_minimal_spacing(model.lattice, args...)
end

@doc raw"""
Selects a kgrid size which ensures that at least a `n_kpoints` total number of ``k``-points
are used. The distribution of ``k``-points amongst coordinate directions is as uniformly
as possible, trying to achieve an identical minimal spacing in all directions.
"""
function kgrid_from_minimal_n_kpoints(lattice, n_kpoints::Integer)
    lattice = austrip.(lattice)
    n_dim   = count(!iszero, eachcol(lattice))
    @assert n_kpoints > 0
    n_kpoints == 1 && return [1, 1, 1]
    n_dim == 1 && return [n_kpoints, 1, 1]

    # Compute truncated reciprocal lattice
    recip_lattice_nD = 2π * inv(lattice[1:n_dim, 1:n_dim]')
    n_kpt_per_dim = n_kpoints^(1/n_dim)

    # Start from a cubic lattice. If it is one, we are done. Otherwise the resulting
    # spacings in each dimension bracket the ideal k-point spacing.
    spacing_per_dim = [norm(vec) / n_kpt_per_dim for vec in eachcol(recip_lattice_nD)]
    min_spacing, max_spacing = extrema(spacing_per_dim)
    if min_spacing ≈ max_spacing
        return kgrid_from_minimal_spacing(lattice, min_spacing)
    else
        number_of_kpoints(spacing) = prod(vec -> norm(vec) / spacing, eachcol(recip_lattice_nD))
        @assert number_of_kpoints(min_spacing) + 0.05 ≥ n_kpoints
        @assert number_of_kpoints(max_spacing) - 0.05 ≤ n_kpoints

        # TODO This is not optimal and sometimes finds too large grids
        spacing = Roots.find_zero(sp -> number_of_kpoints(sp) - n_kpoints,
                                  (min_spacing, max_spacing), Roots.Bisection(),
                                  xatol=1e-4, atol=0, rtol=0)

        # Sanity check: Sometimes root finding is just across the edge towards
        # a larger number of k-points than needed. This attempts a slightly larger spacing.
        kgrid_larger = kgrid_from_minimal_spacing(lattice, spacing + 1e-4)
        if prod(kgrid_larger) ≥ n_kpoints
            return kgrid_larger
        else
            return kgrid_from_minimal_spacing(lattice, spacing)
        end
    end
end
function kgrid_from_minimal_kpoints(model::Model, args...)
    kgrid_from_minimal_kpoints(model.lattice, args...)
end


text:
using ITensorNetworkAD
using AutoHOOT, ITensors, Zygote
using ITensorNetworkAD.ITensorNetworks:
  TreeTensor, uncontract_inds_binary_tree, tree_approximation, inds_binary_tree
using ITensorNetworkAD.ITensorNetworks: inds_network, project_boundary, Models
using ITensorNetworkAD.ITensorAutoHOOT: SubNetwork, batch_tensor_contraction

const itensorah = ITensorNetworkAD.ITensorAutoHOOT

@testset "test TreeTensor" begin
  i = Index(2, "i")
  j = Index(3, "j")
  k = Index(2, "k")
  l = Index(4, "l")
  m = Index(5, "m")

  A = randomITensor(i, j, k)
  B = randomITensor(k, l, m)
  C = randomITensor(i, j, l, m)
  tree_A = TreeTensor(A)
  tree_B = TreeTensor(B)
  tree_C = TreeTensor(C)

  out = A * B
  network = [tree_A, tree_B]
  nodes, dict = itensorah.generate_einsum_expr([network])
  out_list = itensorah.compute_graph(nodes, dict; cutoff=1e-15, maxdim=1000)
  @test isapprox(out, ITensor(out_list[1]))

  out = A * B * C
  out2 = contract(tree_A, tree_B, tree_C; cutoff=1e-15, maxdim=1000)
  @test isapprox(out, ITensor(out2))
end

@testset "test batch_tensor_contraction" begin
  i = Index(2, "i")
  j = Index(3, "j")
  k = Index(2, "k")
  l = Index(2, "l")
  m = Index(2, "m")
  A = randomITensor(i, j)
  B = randomITensor(j, k)
  C = randomITensor(k, l)
  D = randomITensor(l, m)
  E = randomITensor(m, i)

  function network(A)
    tensor_network = SubNetwork(SubNetwork(A, B, C), D, E)
    out = itensorah.batch_tensor_contraction(
      TreeTensor, [tensor_network], A; cutoff=1e-15, maxdim=1000, optimize=false
    )
    return sum(out)[]
  end
  grad_A = gradient(network, A)
  @test isapprox(grad_A[1], B * C * D * E)
end

@testset "test uncontract_inds_binary_tree" begin
  i = Index(2, "i")
  j = Index(2, "j")
  k = Index(2, "k")
  l = Index(2, "l")
  m = Index(2, "m")
  A = randomITensor(i)
  B = randomITensor(j)
  C = randomITensor(k)
  D = randomITensor(l)
  E = randomITensor(m)

  path = [[[A, B], [C, D]], E]
  uncontract_inds = [i, j, k, l, m]
  btree = uncontract_inds_binary_tree(path, uncontract_inds)
  @test btree == [[[[i], [j]], [[k], [l]]], [m]]
  out = tree_approximation([A, B, C, D, E], btree)
  @test isapprox(contract(out), A * B * C * D * E)
end

@testset "test MPS times MPO" begin
  N = (10, 3)
  linkdim = 3
  cutoff = 1e-15
  tn_inds = inds_network(N...; linkdims=linkdim)
  tn = map(inds -> randomITensor(inds...), tn_inds)
  state = 1
  tn = project_boundary(tn, state)
  x, A = tn[:, 1], tn[:, 2]
  out_true = contract(MPO(A), MPS(x); cutoff=cutoff, maxdim=linkdim * linkdim)
  out2 = batch_tensor_contraction(
    TreeTensor,
    [SubNetwork(SubNetwork(x), SubNetwork(A))];
    cutoff=cutoff,
    maxdim=linkdim * linkdim,
  )
  tsr_true = contract(out_true...)
  tsr_nrmsquare = (tsr_true * tsr_true)[1]
  @test isapprox(tsr_true, ITensor(out2[1]))

  maxdims = [2, 4, 6, 8]
  for dim in maxdims
    out = contract(MPO(A), MPS(x); cutoff=cutoff, maxdim=dim)
    out2 = batch_tensor_contraction(
      TreeTensor, [SubNetwork(SubNetwork(x), SubNetwork(A))]; cutoff=cutoff, maxdim=dim
    )
    residual1 = tsr_true - contract(out...)
    residual2 = tsr_true - ITensor(out2[1])
    error1 = sqrt((residual1 * residual1)[1] / tsr_nrmsquare)
    error2 = sqrt((residual2 * residual2)[1] / tsr_nrmsquare)
    print("maxdim, ", dim, ", error1, ", error1, ", error2, ", error2, "\n")
  end
end

@testset "test inds_binary_tree" begin
  i = Index(2, "i")
  j = Index(2, "j")
  k = Index(2, "k")
  l = Index(2, "l")
  m = Index(2, "m")

  T = randomITensor(i, j, k, l, m)
  M = MPS(T, (i, j, k, l, m); cutoff=1e-5, maxdim=5)
  network = M[:]

  out = inds_binary_tree(network, [i, j, k, l, m]; algorithm="mincut")
  @test length(out) == 2
  out = inds_binary_tree(network, [i, j, k, l, m]; algorithm="mps")
  @test length(out) == 2
end

@testset "test PEPS" begin
  N = (8, 8) #(12, 12)
  linkdim = 2
  cutoff = 1e-15
  tn_inds = inds_network(N...; linkdims=linkdim)
  tn = map(inds -> randomITensor(inds...), tn_inds)
  state = 1
  tn = project_boundary(tn, state)

  function get_contracted_peps(maxdim)
    x = tn[:, 1]
    for i in 2:(N[2] - 1)
      A = tn[:, i]
      x = contract(MPO(A), MPS(x); cutoff=cutoff, maxdim=maxdim)[:]
    end
    out = contract(x..., tn[:, N[2]]...)
    network = SubNetwork(tn[:, 1])
    for i in 2:(N[2])
      network = SubNetwork(network, SubNetwork(tn[:, i]))
    end
    out2 = batch_tensor_contraction(TreeTensor, [network]; cutoff=cutoff, maxdim=maxdim)
    return out[], ITensor(out2[1])[]
  end

  ITensors.set_warn_order(100)
  out_true, out2 = get_contracted_peps(linkdim^N[2])
  print(out_true, out2)
  @test abs((out_true - out2) / out_true) < 1e-3

  maxdims = [10, 11, 12, 13, 14, 15, 16] #[2, 4, 8, 16, 24, 32, 40, 48, 56, 64]
  for dim in maxdims
    out, out2 = get_contracted_peps(dim)
    error1 = abs((out - out_true) / out_true)
    error2 = abs((out2 - out_true) / out_true)
    print("maxdim, ", dim, ", error1, ", error1, ", error2, ", error2, "\n")
  end
end


text:
# contains functions for file i/o

"""
    aswsav(file)
    Saves ASWING geometry state to a file
"""
function aswsav(file)
    # split filepath and filename
    filepath, filename = splitdir(file)

    # create path to file and change working directory
    wd = pwd()
    if !isempty(filepath)
        if !isfile(filepath)
            mkpath(filepath)
        end
        cd(filepath)
    end

    # add .asw extension if not already present
    if length(filename) <= 4 || filename[end-3:end] != ".asw"
        filename *= ".asw"
    end
    # write file
    ccall((:outputjl_,libaswing), Nothing,(Ptr{Cchar},Ref{Int32}),filename,length(filename))
    # change back to working directory
    cd(wd)
    return nothing
end

"""
    pntsav(file, ipnt1=1, ipnt2=ASWING.NPOINT)
    Saves ASWING operating point state to a file
"""
function pntsav(file::String, ipnt1=1, ipnt2=ASWING.NPOINT)
    # split filepath and filename
    filepath, filename = splitdir(file)
    # create path to file and change working directory
    wd = pwd()
    if !isempty(filepath)
        if !isfile(filepath)
            mkpath(filepath)
        end
        cd(filepath)
    end
    # add .pnt extension if not already present
    if length(filename) <= 4 || filename[end-3:end] != ".pnt"
        filename *= ".pnt"
    end
    # write file
    ccall((:pntsavjl_, libaswing), Nothing, (Ptr{Cchar}, Ref{Int32}, Ref{Int32},
        Ref{Int32}), filename, length(filename), ipnt1, ipnt2)
    # change back to working directory
    cd(wd)
    return nothing
end

"""
    setsav(file)
    Saves ASWING operating point state to a file
"""
function setsav(file)
    # split filepath and filename
    filepath, filename = splitdir(file)
    # create path to file and change working directory
    wd = pwd()
    if !isempty(filepath)
        if !isfile(filepath)
            mkpath(filepath)
        end
        cd(filepath)
    end
    # add .set extension if not already present
    if length(filename) <= 4 || filename[end-4:end] != ".set"
        filename *= ".set"
    end
    # write file
    ccall((:setsavjl_,libaswing), Nothing,(Ptr{Cchar},Ref{Int32}),filename,length(filename))
    # change back to working directory
    cd(wd)
    return nothing
end

"""
    readasw(filename::AbstractString, R::Type{<:Real}=Float64)

    Reads a .asw file
"""
function readasw(filename::AbstractString, R::Type{<:Real}=Float64)
    local aswname, si_units, constants, refvals
    sensors = Sensor{R}[]
    engines = Engine{R}[]
    weights = Weight{R}[]
    struts = Strut{R}[]
    joints = Joint{R}[]
    jangles = Jangle{R}[]
    grounds = Ground{R}[]
    beams = Beam{R}[]
    # Open File
    f = open(filename)
    line = 0
    nbeam = 0
    njang = 0
    # Loop through file line by line
    while !eof(f)
        str = readline(f)
        line = line+1
        if isinline(lowercase(str),"end",1:3)
            break
        elseif isinline(str,"!",1:1) || isinline(str,"#",1:1) || isinline(str,"%",1:1)
        elseif strip(str) == ""
            # Check for and read name block
        elseif isinline(lowercase(str),"name",1:4)
            # println("Reading Name")
            while !eof(f)
                str = readline(f)
                line = line+1
                if isinline(lowercase(str),"end",1:3)
                    break
                elseif isinline(str,"!",1:1) || isinline(str,"#",1:1) || isinline(str,"%",1:1)
                else
                    aswname = chomp(str)
                end
            end
            # Check for and read unit block
        elseif isinline(lowercase(str),"unit",1:4)
            # println("Reading Units")
            while !eof(f)
                str = readline(f)
                line = line+1
                if isinline(lowercase(str),"end",1:3)
                    break
                elseif isinline(str,"!",1:1) || isinline(str,"#",1:1) || isinline(str,"%",1:1)
                elseif strip(str) == ""
                else
                    substr = split(str)
                    if lowercase(substr[1]) == "l" && lowercase(substr[3]) == "ft"
                        si_units = false
                    elseif lowercase(substr[1]) == "m"
                    elseif lowercase(substr[1]) == "f" && lowercase(substr[3]) == "lb"
                        si_units = false
                    elseif lowercase(substr[1]) == "t"
                    else
                        si_units = true
                    end
                end
            end
        # Check for and read constant block
        elseif isinline(lowercase(str),"cons",1:4)
            # println("Reading Constants")
            nvalue = 3
            skip = 0
            constin, line = parseblock(f, line, nvalue, skip)
            constants = Constants{R}(;
                gravity = constin[1],
                sea_level_air_density = constin[2],
                sea_level_speed_of_sound = constin[3])
            # Check for and read reference block
        elseif isinline(lowercase(str),"refe",1:4)
            # println("Reading Reference Values")
            nvalue = 3
            skip = 0
            refin,line = parseblock(f, line, nvalue, skip)
            area = refin[1,1]
            chord = refin[1,2]
            span = refin[1,3]
            moments = zeros(R,3)
            accelerations = zeros(R,3)
            velocities = zeros(R,3)
            for i = 1:3
                moments[i] = refin[2,i]
                accelerations[i] = refin[3,i]
                velocities[i] = refin[4,i]
            end
            refvals = ReferenceValues{R}(area, chord, span, moments, accelerations,
                velocities)
            # Check for and read sensor block
        elseif isinline(lowercase(str), "sens", 1:4)
            # println("Reading Sensor Values")
            nvalue = 12
            skip = 2
            sensin,line = parseblock(f, line, nvalue, skip)
            for i = 1:size(sensin,1)
                sensor = Sensor{R}(Int(sensin[i,1]), Int(sensin[i,2]), sensin[i,3],
                    [sensin[i,4], sensin[i,5], sensin[i,6]],
                    [sensin[i,7], sensin[i,8], sensin[i,9]],
                    [sensin[i,10], sensin[i,11], sensin[i,12]])
                push!(sensors, sensor)
            end
            # Check for and read engine block
        elseif isinline(lowercase(str), "engi", 1:4)
            # println("Reading Engine Values")
            nvalue = 25
            skip = 3
            enginein, line = parseblock(f, line, nvalue, skip)
            for i = 1:size(enginein, 1)
                engine = Engine{R}(Int(enginein[i,1]),
                    Int(enginein[i,2]), Int(enginein[i,3]),
                    enginein[i,4], [enginein[i,5], enginein[i,6], enginein[i,7]],
                    [enginein[i,8], enginein[i,9], enginein[i,10]],
                    enginein[i,11], enginein[i,12], enginein[i,13], enginein[i,14],
                    enginein[i,15], enginein[i,16], enginein[i,17], enginein[i,18],
                    enginein[i,19], enginein[i,20], enginein[i,21], enginein[i,22],
                    enginein[i,23], enginein[i,24], enginein[i,25])
                push!(engines,engine)
            end
            # Check for and read weight block
        elseif isinline(lowercase(str), "weig", 1:4)
            # println("Reading Weight Values")
            nvalue = 17
            skip = 1
            weightin, line = parseblock(f, line, nvalue, skip)
            for i = 1:size(weightin, 1)
                weight = Weight{R}(Int(weightin[i,1]), weightin[i,2],
                    [weightin[i,3], weightin[i,4], weightin[i,5]],
                    weightin[i,6], weightin[i,7], weightin[i,8],
                    [weightin[i,9], weightin[i,10], weightin[i,11]],
                    weightin[i,12], weightin[i,13], weightin[i,14],
                    weightin[i,15], weightin[i,16], weightin[i,17])
                push!(weights, weight)
            end
            # Check for and read strut block
        elseif isinline(lowercase(str), "stru", 1:4)
            # println("Reading Strut Block")
            nvalue = 10
            skip = 1
            strutin,line = parseblock(f, line, nvalue, skip)
            for i = 1:size(strutin, 1)
                strut = Strut{R}(strutin[i,1], strutin[i,2],
                    [strutin[i,3], strutin[i,4], strutin[i,5]],
                    [strutin[i,6], strutin[i,7], strutin[i,8]],
                    strutin[i,9], strutin[i,10])
                push!(struts,strut)
            end
            # Check for and read joint block
        elseif isinline(lowercase(str), "join", 1:4)
            # println("Reading Joint Block")
            nvalue = 5
            skip = 2
            jointin, line = parseblock(f, line, nvalue, skip)
            for i = 1:size(jointin,1)
                joint = Joint{R}(Int(jointin[i,1]), Int(jointin[i,2]), jointin[i,3],
                    jointin[i,4], Int(jointin[i,5]))
                push!(joints, joint)
            end
            # Check for and read jangle block
        elseif isinline(lowercase(str), "jang", 1:4)
            # println("Reading Jangle Block")
            nvalue = 4
            njang = njang+1
            while !eof(f)
                str = readline(f)
                line = line + 1
                if lowercase(str[1:3]) == "end"
                    break
                elseif (str[1] == "!" || str[1] == "#" || str[1] == "%")
                elseif strip(str) == ""
                else
                    substr = split(str)
                    rdat = zeros(nvalue)
                    for i = 1:length(substr)
                        rdat[i] = Meta.parse(substr[i])
                    end
                    Njoint = rdat[1]
                    hinge_axis = [rdat[2], rdat[3], rdat[4]]
                    nvalue = 2
                    skip = 0
                    jangle_data, line = parseblock(f, line, nvalue, skip)
                    for i = 1:size(jangle_data,1)
                        hinge_moment[i] = jangle_data[i,1]
                        hinge_angle[i] = jangle_data[i,2]
                    end
                    jangle = Jangle{R}(Njoint, hinge_axis, hinge_moment, hinge_angle)
                    push!(jangles, jangle)
                    break
                end
            end
            # Check for and read ground block
        elseif lowercase(str[1:4]) == "grou"
            # println("Reading Ground Block")
            nvalue = 3
            skip = 1
            groundin,line = parseblock(f, line, nvalue, skip)
            for i = 1:size(groundin,1)
                ground = Ground{R}(groundin[i,1], groundin[i,2], groundin[i,3])
                push!(grounds, ground)
            end
            # Check for and read new beam
        elseif lowercase(str[1:4]) == "beam"
            # println("Reading Beam Blocks")
            nbeam = nbeam+1
            beamvars = [Array{R}(undef, 0, 2) for idx in 1:length(beamvarval)]
            # Read beam number and physical index
            substr = split(chomp(str))
            beam_number = Meta.parse(substr[2])
            if length(substr) > 2
                physical_index = Meta.parse(substr[3])
            else
                physical_index = beam_number
            end
            # Read beam name
            str = readline(f)
            line = line+1
            beam_name = chomp(str)
            # Read beam data
            str = readline(f)
            line = line+1
            beam_data_end = false
            while !eof(f) && !beam_data_end
                if isinline(lowercase(str),"end",1:3)
                    break
                elseif isinline(str,"!",1:1) || isinline(str,"#",1:1) || isinline(str,"%",1:1)
                    str = readline(f)
                elseif strip(str) == ""
                    str = readline(f)
                else
                    bvar = split(chomp(str))
                    if bvar[1] != "t"
                        error("First variable must be \"t\"")
                    end
                    for i = 2:length(bvar)
                        if ~haskey(beamvarval, bvar[i])
                            @warn("Variable ", bvar[i] , " is not a valid variable", )
                        end
                    end
                    nvalue = length(split(str))
                    qadd = zeros(nvalue)
                    qfac = ones(nvalue)
                    iline = 1
                    while !eof(f)
                        str = readline(f)
                        idx = findfirst(isequal('!'), str[2:end])
                        if idx != nothing
                            str = str[1:idx]
                        end
                        line = line + 1
                        if occursin("end", lowercase(str))
                            beam_data_end = true
                            break
                        elseif isinline(strip(str), "t", 1:1)
                            break
                        elseif isinline(str,"!",1:1) || isinline(str,"#",1:1) || isinline(str,"%",1:1)
                        elseif strip(str) == ""
                        elseif isinline(str,"*",1:1)
                            substr = split(str[2:end])
                            for i = 1:length(substr)
                                qfac[i] = Meta.parse(substr[i])
                            end
                        elseif isinline(str, "+", 1:1)
                            substr = split(str[2:end])
                            for i = 1:length(substr)
                                qadd[i] = Meta.parse(substr[i])
                            end
                        else
                            substr = split(str)
                            rdat = zeros(nvalue)
                            for i = 1:length(substr)
                                rdat[i] = Meta.parse(substr[i])
                            end
                            clbvalin = rdat.*qfac+qadd
                            for i = 2:length(bvar)
                                for var in keys(beamvarval)
                                    if bvar[i] == var
                                        beamvars[beamvarval[var]] = vcat(beamvars[beamvarval[var]],
                                            [clbvalin[1] clbvalin[i]])
                                    end
                                end
                            end
                            iline = iline+1
                        end
                    end
                end
            end
            spanwise_variables = SpanwiseVariables{R}(beamvars...)
            beam = Beam{R}(beam_number, physical_index, beam_name, spanwise_variables)
            push!(beams, beam)
        else
            error("Undefined input on line $line:\n
            $str")
        end
    end
    close(f)
    asw = Configuration{R}(aswname, si_units, constants, refvals, weights,
        sensors, engines, struts, joints, jangles, grounds, beams)
    return asw
end

"""
    parseblock(f, line, nvalue, skip)

    parses common input blocks, used in readasw
"""
function parseblock(f, line, nvalue, skip)
    qadd = zeros(nvalue)
    qfac = ones(nvalue)
    fdata = zeros(1,nvalue)
    iline = 1
    while !eof(f)
        str = readline(f)
        idx = findfirst(isequal('!'), str[2:end])
        if idx != nothing
            str = str[1:idx]
        end
        line = line + 1
        if isinline(lowercase(str),"end",1:3)
            break
        elseif isinline(str,"!",1:1) || isinline(str,"#",1:1) || isinline(str,"%",1:1)
        elseif strip(str) == ""
        elseif isinline(str,"*",1:1)
            substr = split(str[2:end])
            for i = 1:length(substr)
                qfac[i] = Meta.parse(substr[i])
            end
        elseif isinline(str,"+",1:1)
            substr = split(str[2:end])
            for i = 1:length(substr)
                qadd[i] = Meta.parse(substr[i])
            end
        else
            substr = split(str)
            rdat = zeros(nvalue)
            for i = 1:length(substr)
                rdat[i] = Meta.parse(substr[i])
            end
            flinedata = zeros(1,nvalue)
            for i = 1:length(substr)
                if i<=skip
                    flinedata[1,i] = rdat[i]
                else
                    flinedata[1,i] = rdat[i]*qfac[i-skip]+qadd[i-skip]
                end
            end
            if iline == 1
                fdata = flinedata
            else
                fdata = vcat(fdata,flinedata)
            end
            iline = iline+1
        end
    end
    return fdata,line
end

"""
    isinline(line, stringtofind, range)

    checks for string within line, used in readasw
"""
function isinline(line, stringtofind, range)
    if length(line) < maximum(range)
        return false
    else
        return occursin(line[range], stringtofind)
    end
end

"""
    readpnt(filename::String)

    Reads a .pnt file
"""
function readpnt(filename::AbstractString)
    local constraints, mach_from_airspeed, ground_image, parameters
    open(filename) do f
        oper_pnts = Array{OperatingPoint,1}(undef, 0)
        str = readline(f)
        substr = split(str)
        iptot0 = Meta.parse(substr[1])
        kptot0 = Meta.parse(substr[2])
        kpfree0 = Meta.parse(substr[3])

        if (iptot0 != IPTOT) || (kptot0 != KPTOT) || (kpfree0 != KPFREE)
            error("Dimensions have changed. Point file not read.")
        end

        point = 0
        while !eof(f)
            point = point+1

            # Read constraints
            str = readline(f)
            substr = split(str)

            idx = 0

            linear_acceleration = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            angular_acceleration = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            velocity = Meta.parse(substr[idx+1])
            idx += 1

            beta = Meta.parse(substr[idx+1])
            idx += 1

            alpha = Meta.parse(substr[idx+1])
            idx += 1

            rotation_rate = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            position = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            phi = Meta.parse(substr[idx+1])
            idx += 1

            theta = Meta.parse(substr[idx+1])
            idx += 1

            psi = Meta.parse(substr[idx+1])
            idx += 1

            flap_defl_ctrl_var = zeros(Int, NFLPX)
            for i = 1:NFLPX
                flap_defl_ctrl_var[i] = Meta.parse(substr[idx+i])
            end
            idx += NFLPX

            eng_pwr_ctrl_var = zeros(Int, NENGX)
            for i = 1:NENGX
                eng_pwr_ctrl_var[i] = Meta.parse(substr[idx+i])
            end
            idx += NENGX

            err_int_Vinf = Meta.parse(substr[idx+1])
            idx += 1

            err_int_beta = Meta.parse(substr[idx+1])
            idx += 1

            err_int_alpha = Meta.parse(substr[idx+1])
            idx += 1

            err_int_phi = Meta.parse(substr[idx+1])
            idx += 1

            err_int_theta = Meta.parse(substr[idx+1])
            idx += 1

            err_int_psi = Meta.parse(substr[idx+1])
            idx += 1

            err_int_ROTxyz = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            err_int_VACxyz = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            constraints = Constraints(linear_acceleration,
                angular_acceleration, [velocity, beta, alpha], rotation_rate,
                position, phi, theta, psi, flap_defl_ctrl_var, eng_pwr_ctrl_var,
                err_int_Vinf, err_int_beta, err_int_alpha, err_int_phi,
                err_int_theta, err_int_psi, err_int_ROTxyz, err_int_VACxyz)

            # Read mach_from_airspeed boolean and ground image integer
            str = readline(f)
            substr = split(str)
            if substr[1] == "T"
                mach_from_airspeed = true
            else
                mach_from_airspeed = false
            end
            ground_image = Meta.parse(substr[2])

            # Read specified parameters
            str = readline(f)
            substr = split(str)

            idx = 0

            linear_acceleration = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            angular_acceleration = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            velocity = Meta.parse(substr[idx+1])
            idx += 1

            beta = Meta.parse(substr[idx+1])
            idx += 1

            alpha = Meta.parse(substr[idx+1])
            idx += 1

            rotation_rate = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            position = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            phi = Meta.parse(substr[idx+1])
            idx += 1

            theta = Meta.parse(substr[idx+1])
            idx += 1

            psi = Meta.parse(substr[idx+1])
            idx += 1

            flap_defl_ctrl_var = zeros(Float64, NFLPX)
            for i = 1:NFLPX
                flap_defl_ctrl_var[i] = Meta.parse(substr[idx+i])
            end
            idx += NFLPX

            eng_pwr_ctrl_var = zeros(Float64, NENGX)
            for i = 1:NENGX
                eng_pwr_ctrl_var[i] = Meta.parse(substr[idx+i])
            end
            idx += NENGX

            sum_force = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            sum_mom = [Meta.parse(substr[idx+1]), Meta.parse(substr[idx+2]), Meta.parse(substr[idx+3])]
            idx += 3

            lift = Meta.parse(substr[idx+1])
            idx += 1

            climb_angle = Meta.parse(substr[idx+1])
            idx += 1

            radial_acceleration = Meta.parse(substr[idx+1])
            idx += 1

            usr1 = Meta.parse(substr[idx+1])
            idx += 1

            usr2 = Meta.parse(substr[idx+1])
            idx += 1

            least_squared = Meta.parse(substr[idx+1])
            idx += 1

            parameters = Parameters(
                linear_acceleration = linear_acceleration,
                angular_acceleration = angular_acceleration,
                velocity = velocity,
                beta = beta,
                alpha = alpha,
                rotation_rate = rotation_rate,
                position = position,
                phi = phi,
                theta = theta,
                psi = psi,
                flap_defl_ctrl_var = flap_defl_ctrl_var,
                eng_pwr_ctrl_var = eng_pwr_ctrl_var,
                sum_force = sum_force,
                sum_mom = sum_mom,
                lift = lift,
                climb_angle = climb_angle,
                radial_acceleration = radial_acceleration,
                usr1 = usr1,
                usr2 = usr2,
                least_squared = least_squared)
            # Global parameters set inside ASWING when solved so these are not needed
            str = readline(f)
            # Read flap weights for least squares
            str = readline(f)
            substr = split(str)
            wflap = zeros(Float64, NFLPX)
            for i = 1:min(length(substr), NFLPX)
                wflap[i] = Meta.parse(substr[i])
            end
            # Read engine weights for least squares
            str = readline(f)
            substr = split(str)
            wpeng = zeros(Float64, NENGX)
            for i = 1:min(length(substr), NENGX)
                wpeng[i] = Meta.parse(substr[i])
            end
            oper_pnt = OperatingPoint(constraints, mach_from_airspeed, 0.0,
                ground_image, parameters, wflap, wpeng)
            push!(oper_pnts,oper_pnt)
        end
        return oper_pnts
    end
end


text:
module Dummy

using SCIP

# Define a minimal no-op constraint handler.
# Needs to be mutable for `pointer_from_objref` to work.
mutable struct DummyConsHdlr <: SCIP.AbstractConstraintHandler
    check_called::Int64
    enfo_called::Int64
    lock_called::Int64

    DummyConsHdlr() = new(0, 0, 0)
end

# Implement only the fundamental callbacks:
function SCIP.check(ch::DummyConsHdlr,
                    constraints::Array{Ptr{SCIP.SCIP_CONS}},
                    sol::Ptr{SCIP.SCIP_SOL},
                    checkintegrality::Bool,
                    checklprows::Bool,
                    printreason::Bool,
                    completely::Bool)
    ch.check_called += 1
    return SCIP.SCIP_FEASIBLE
end

function SCIP.enforce_lp_sol(ch::DummyConsHdlr,
                             constraints::Array{Ptr{SCIP.SCIP_CONS}},
                             nusefulconss::Cint,
                             solinfeasible::Bool)
    ch.enfo_called += 1
    return SCIP.SCIP_FEASIBLE
end

function SCIP.enforce_pseudo_sol(ch::DummyConsHdlr,
                                 constraints::Array{Ptr{SCIP.SCIP_CONS}},
                                 nusefulconss::Cint,
                                 solinfeasible::Bool,
                                 objinfeasible::Bool)
    ch.enfo_called += 1
    return SCIP.SCIP_FEASIBLE
end

function SCIP.lock(ch::DummyConsHdlr,
                   constraint::Ptr{SCIP.SCIP_CONS},
                   locktype::SCIP.SCIP_LOCKTYPE,
                   nlockspos::Cint,
                   nlocksneg::Cint)
    ch.lock_called += 1
end

# Corresponding, empty constraint (data) object
mutable struct DummyCons <: SCIP.AbstractConstraint{DummyConsHdlr} end

end # module Dummy

module NeverSatisfied

using SCIP

mutable struct NSCH <: SCIP.AbstractConstraintHandler
    check_called::Int64
    enfo_called::Int64
    lock_called::Int64

    NSCH() = new(0, 0, 0)
end

function SCIP.check(ch::NSCH, constraints, sol, checkintegrality,
                    checklprows, printreason, completely)
    ch.check_called += 1
    return SCIP.SCIP_INFEASIBLE
end

function SCIP.enforce_lp_sol(ch::NSCH, constraints, nusefulconss, solinfeasible)
    ch.enfo_called += 1
    return SCIP.SCIP_INFEASIBLE
end

function SCIP.enforce_pseudo_sol(ch::NSCH, constraints, nusefulconss,
                                 solinfeasible, objinfeasible)
    ch.enfo_called += 1
    return SCIP.SCIP_INFEASIBLE
end

function SCIP.lock(ch::NSCH, constraint, locktype, nlockspos, nlocksneg)
    ch.lock_called += 1
end

# Corresponding, empty constraint (data) object
mutable struct Cons <: SCIP.AbstractConstraint{NSCH} end

end # module AlwaysSatisfied

module NaiveAllDiff

using MathOptInterface
using SCIP

const MOI = MathOptInterface

mutable struct NADCH <: SCIP.AbstractConstraintHandler
    scip::SCIP.Optimizer # for SCIP* and var maps
end

# Constraint data, referencing variables of a single constraint.
mutable struct NADCons <: SCIP.AbstractConstraint{NADCH}
    variables::Array{MOI.VariableIndex}
end

# Helper function used in several callbacks
function anyviolated(ch, constraints, sol=C_NULL)
    for cons_::Ptr{SCIP.SCIP_CONS} in constraints
        cons::NADCons = SCIP.user_constraint(cons_)

        # extract solution values
        values = SCIP.sol_values(ch.scip, cons.variables, sol)

        # check for constraint violation
        if !allunique(values)
            return true
        end
    end
    return false
end

function SCIP.check(ch::NADCH, constraints, sol, checkintegrality,
                    checklprows, printreason, completely)
    if anyviolated(ch, constraints, sol)
        return SCIP.SCIP_INFEASIBLE
    else
        return SCIP.SCIP_FEASIBLE
    end
end

function SCIP.enforce_lp_sol(ch::NADCH, constraints, nusefulconss, solinfeasible)
    if anyviolated(ch, constraints, C_NULL)
        return SCIP.SCIP_INFEASIBLE
    else
        return SCIP.SCIP_FEASIBLE
    end
end

function SCIP.enforce_pseudo_sol(ch::NADCH, constraints, nusefulconss,
                                 solinfeasible, objinfeasible)
    if anyviolated(ch, constraints, C_NULL)
        return SCIP.SCIP_INFEASIBLE
    else
        return SCIP.SCIP_FEASIBLE
    end
end

function SCIP.lock(ch::NADCH, constraint, locktype, nlockspos, nlocksneg)
    cons::NADCons = SCIP.user_constraint(constraint)

    # look all variables in both directions
    for vi in cons.variables
        # TODO: understand why lock is called during SCIPfree, after the
        # constraint should have been deleted already. Does it mean we should
        # implement CONSTRANS?
        var_::Ptr{SCIP.SCIP_VAR} = SCIP.var(ch.scip, vi)
        var_ != C_NULL || continue  # avoid segfault!

        SCIP.@SC SCIP.SCIPaddVarLocksType(
            ch.scip, var_, locktype, nlockspos + nlocksneg, nlockspos + nlocksneg)
    end
end

end # module NaiveAllDiff

module NoGoodCounter

using MathOptInterface
using SCIP

const MOI = MathOptInterface

mutable struct Counter <: SCIP.AbstractConstraintHandler
    scip::SCIP.Optimizer # for SCIP* and var maps
    variables::Array{MOI.VariableIndex}
    solutions::Set{Array{Float64}}

    Counter(scip, variables) = new(scip, variables, Set())
end

function SCIP.check(ch::Counter, constraints, sol, checkintegrality,
                    checklprows, printreason, completely)
    return SCIP.SCIP_INFEASIBLE
end

function enforce(ch::Counter)
    values = SCIP.sol_values(ch.scip, ch.variables)

    # Store solution
    if values in ch.solutions
        # Getting the same solution twice?
        return SCIP.SCIP_INFEASIBLE
    end
    push!(ch.solutions, values)

    # Add no-good constraint: sum_zeros(x) + sum_ones(1-x) >= 1
    zeros = values .≈ 0.0
    ones = values .≈ 1.0
    others = .!zeros .& .!ones
    if any(others)
        error("Found non-binary solution value for $(vars[others])")
    end

    terms = vcat([MOI.ScalarAffineTerm( 1.0, vi) for vi in ch.variables[zeros]],
                 [MOI.ScalarAffineTerm(-1.0, vi) for vi in ch.variables[ones]])
    ci = MOI.add_constraint(ch.scip, MOI.ScalarAffineFunction(terms, 0.0),
                            MOI.GreaterThan(1.0 - sum(ones)))

    return SCIP.SCIP_CONSADDED
end

function SCIP.enforce_lp_sol(ch::Counter, constraints, nusefulconss,
                             solinfeasible)
    @assert length(constraints) == 0
    return enforce(ch)
end

function SCIP.enforce_pseudo_sol(ch::Counter, constraints, nusefulconss,
                                 solinfeasible, objinfeasible)
    @assert length(constraints) == 0
    return enforce(ch)
end

function SCIP.lock(ch::Counter, constraint, locktype, nlockspos, nlocksneg)
end

end # module NoGoodCounter


text:
function matmul()

    # Create data
    n = future(Int32(15e2))
    m = future(Int32(10e3))
    p = future(Int32(5e3))
    A = future() # n x m
    B = future() # m x p

    # Where the data is located
    val(n)
    val(m)
    val(p)
    mem(A, Int64(evaluate(n) * evaluate(m)), Float64)
    mem(B, Int64(evaluate(m) * evaluate(p)), Float64)

    # How the data is partitioned
    pt(n, Div())
    pt(m, Replicate())
    pt(p, Div())
    pt(A, Block(1))
    pt(B, Block(2))
    pc(Cross(A, B))
    pc(Co(A, n))
    pc(Co(B, p))

    mutated(A)
    mutated(B)

    @partitioned A B n m p begin
        A = fill(1, (Int64(n), Int64(m)))
        B = fill(2, (Int64(m), Int64(p)))
    end

    C = future()
    mem(C, Int64(evaluate(n) * evaluate(p)), Float64)

    pt(A, Block(1))
    pt(B, Block(2))
    pt(C, [Block(1), Block(2)])
    pc(Cross(A, B))
    pc(Cross((C, 1), (C, 2))) # Redundant but currently required to assert order of splitting
    pc(Equal((C, 1), (C, 2)))
    pc(Co((C, 1), A))
    pc(Co((C, 2), B))

    mutated(C)

    @partitioned A B C begin
        C = A * B
    end

    # C_new = future()
    # located(C_new, C.location)
    # mutated(C_new)

    # pt(C, [Block(1), Block(2)])
    # pc(Cross((C, 1), (C, 2)))
    # pt(C_new, BlockMulti([1, 2]))

    # @partitioned C C_new begin
    #     C_new = C
    # end
    # C = C_new

    evaluate(C)
end

@testset "Matrix multiplication" begin
    runtest("Matrix Multiplication", j -> begin
        matmul()
    end)
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text: include("avio.jl")


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

mutable struct BinaryOpReq <: ProtoType
    i1::Int64
    i2::Int64
    BinaryOpReq(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type BinaryOpReq
hash(v::BinaryOpReq) = ProtoBuf.protohash(v)
isequal(v1::BinaryOpReq, v2::BinaryOpReq) = ProtoBuf.protoisequal(v1, v2)
==(v1::BinaryOpReq, v2::BinaryOpReq) = ProtoBuf.protoeq(v1, v2)

mutable struct BinaryOpResp <: ProtoType
    result::Int64
    BinaryOpResp(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type BinaryOpResp
hash(v::BinaryOpResp) = ProtoBuf.protohash(v)
isequal(v1::BinaryOpResp, v2::BinaryOpResp) = ProtoBuf.protoisequal(v1, v2)
==(v1::BinaryOpResp, v2::BinaryOpResp) = ProtoBuf.protoeq(v1, v2)

# service methods for TestMath
const _TestMath_methods = MethodDescriptor[
        MethodDescriptor("Mul", 1, BinaryOpReq, BinaryOpResp),
        MethodDescriptor("Add", 2, BinaryOpReq, BinaryOpResp)
    ] # const _TestMath_methods
const _TestMath_desc = ServiceDescriptor("TestMath", 1, _TestMath_methods)

TestMath(impl::Module) = ProtoService(_TestMath_desc, impl)

mutable struct TestMathStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    TestMathStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_TestMath_desc, channel))
end # type TestMathStub

mutable struct TestMathBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    TestMathBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_TestMath_desc, channel))
end # type TestMathBlockingStub

Mul(stub::TestMathStub, controller::ProtoRpcController, inp::BinaryOpReq, done::Function) = call_method(stub.impl, _TestMath_methods[1], controller, inp, done)
Mul(stub::TestMathBlockingStub, controller::ProtoRpcController, inp::BinaryOpReq) = call_method(stub.impl, _TestMath_methods[1], controller, inp)

Add(stub::TestMathStub, controller::ProtoRpcController, inp::BinaryOpReq, done::Function) = call_method(stub.impl, _TestMath_methods[2], controller, inp, done)
Add(stub::TestMathBlockingStub, controller::ProtoRpcController, inp::BinaryOpReq) = call_method(stub.impl, _TestMath_methods[2], controller, inp)

export BinaryOpReq, BinaryOpResp, TestMath, TestMathStub, TestMathBlockingStub, Mul, Add


text:
type Population
    # Manages all the species
    population::Vector{Chromosome}
    popsize::Int
    species::Vector{Species}
    species_history::Vector{Array{Int,1}}
    generation::Int
    avg_fitness::Vector{Float64}
    best_fitness::Vector{Chromosome}
    evaluate::Function # Evaluates population. Override this method in your experiments
    function Population(g::Global, checkpoint::String="")

        if checkpoint != ""
            # start from a previous point: creates an 'empty'
            # population and point its __dict__ to the previous one
            resume_checkpoint(checkpoint)
        else
            population = Chromosome[]
            popsize = g.cf.pop_size # total population size
            for i in 1:popsize

                ch = g.cf.fully_connected? create_fully_connected(g) : create_minimally_connected(g)
                if g.cf.hidden_nodes > 0
                    ch = ch = g.cf.fully_connected? ch: create_unconnected(g)
                    nType = g.cf.feedforward? FeedForward():Recurrent()
                    add_hidden_nodes!(g, ch, g.cf.hidden_nodes,nType)
                end
                push!(population, ch)
            end

            p = new(population, popsize,
                    Species[],Array{Int,1}[], # currently living species and species history
                    -1, # generation
                    Float64[], # avg_fitness
                    Chromosome[]); # best_fitness
            p.evaluate = (f::Function) -> f(p.population); # Evaluates population. Override this method in your experiments
            p
        end
    end
end

function Base.show(io::IO, p::Population)
    @printf(io,"Population size: %3d   Total species: %3d", p.popsize, length(p.species))
end

function remove(p::Population, ch::Chromosome)
    # Removes a chromosome from the population
    deleteat!(p.population,findfirst(p.population,ch))
    return
end

function speciate(g::Global, p::Population, report::Bool)
    # Group chromosomes into species by similarity
        # Speciate the population
    for individual in p.population
        found = false
        for s in p.species
            if distance(individual, s.representant, g.cf) < g.cf.compatibility_threshold
                add(s, individual)
                found = true
                break
            end
        end

        if !found push!(p.species, Species(g, individual)) end
    end

    # eliminate empty species
    keep = map(s->length(s)==0?false:true,p.species)
    if report
        for i = 1:length(keep)
            if !keep[i] println("Removing species $(p.species[i].id) for being empty") end
        end
    end
    p.species = p.species[keep]
    set_compatibility_threshold(g, p)

end

function set_compatibility_threshold(g::Global, p::Population)
    # controls compatibility threshold
    if length(p.species) > g.cf.species_size
        g.cf.compatibility_threshold += g.cf.compatibility_change
    elseif length(p.species) < g.cf.species_size
        if g.cf.compatibility_threshold > g.cf.compatibility_change
            g.cf.compatibility_threshold -= g.cf.compatibility_change
        else
            println("Compatibility threshold cannot be changed (minimum value has been reached)")
        end
    end
end

# Returns the average raw fitness of population
average_fitness(p::Population) = mean([p.population[i].fitness::Float64 for i=1:length(p.population)])

stdeviation(p::Population) = std([p.population[i].fitness::Float64 for i=1:length(p.population)])

function compute_spawn_levels(g::Global, p::Population)
    #  Compute each species' spawn amount (Stanley, p. 40)

    # 1. Boost if young and penalize if old
    # TODO: does it really increase the overall performance?
    species_stats = zeros(length(p.species))
    for i = 1:length(p.species)
        s = p.species[i]
        species_stats[i] = s.age < g.cf.youth_threshold? average_fitness(s) * g.cf.youth_boost:
            s.age > g.cf.old_threshold? average_fitness(s) * g.cf.youth_boost : average_fitness(s)
        if s.age < g.cf.youth_threshold
            species_stats[i] = average_fitness(s) * g.cf.youth_boost
        elseif s.age > g.cf.old_threshold
           species_stats[i] = average_fitness(s) * g.cf.old_penalty
        else
            species_stats[i] = average_fitness(s)
        end
    end

    # 2. Share fitness (only usefull for computing spawn amounts)
    # More info: http://tech.groups.yahoo.com/group/neat/message/2203
    # Sharing the fitness is only meaningful here
    # we don't really have to change each individual's raw fitness
    total_average = sum(species_stats)

     # 3. Compute spawn
    for i= 1:length(p.species)
        s = p.species[i]
        s.spawn_amount = int(round((species_stats[i] * p.popsize / total_average)))
    end
end

function tournamentSelection(p::Population, k=2)
    # Tournament selection with size k (default k=2).
    # randomly select k competitors
    chs = p.population[randperm(length(p.population))[1:k]]
    best = chs[1]
    for ch in chs # choose best among randomly selected
        best = ch.fitness > best.fitness? ch :best
    end
    return best
end

function log_species(p::Population)
    # Logging species data for visualizing speciation
    specById = sort(p.species, by=s->s.id)
    spec_size = zeros(Int,specById[end].id+1)
    map(s->spec_size[s.id]=length(s) , specById)
    push!(p.species_history,spec_size)
end

function population_diversity(p::Population)
    # Calculates the diversity of population: total average weights,
    # number of connections, nodes

    num_nodes = 0
    num_conns = 0
    avg_weights = 0.0

    for ch in p.population
        num_nodes += length(ch.node_genes)
        num_conns += length(ch.conn_genes)
        for cg in ch.conn_genes
            avg_weights += cg.weight
        end
    end

    total = length(p.population)
    return num_nodes/total, num_conns/total, avg_weights/total
end

function epoch(g::Global, p::Population, n::Int, report::Bool=true, save_best::Bool=false,
               checkpoint_interval::Int=15, checkpoint_generation=0)
    #= Runs NEAT's genetic algorithm for n epochs.

        Keyword arguments:
        report -- show stats at each epoch (default true)
        save_best -- save the best chromosome from each epoch (default False)
        checkpoint_interval -- time in minutes between saving checkpoints (default 15 minutes)
        checkpoint_generation -- time in generations between saving checkpoints
            (default 0 -- option disabled)
    =#
    t0 = time() # for saving checkpoints

    for gen in 0:n
        p.generation += 1

        if report println(" ****** Running generation $(p.generation) ******") end

        # Evaluate individuals
        p.evaluate(p.population)
        # Speciates the population
        speciate(g, p, report)

        # Current population's average fitness
        push!(p.avg_fitness, average_fitness(p))

        # Current generation's best chromosome
        bestfit, bestidx = findmax(map(ch-> ch.fitness, p.population))
        best = p.population[bestidx]
        push!(p.best_fitness, best)

        # Which species has the best chromosome?
        for s in p.species
            s.hasBest = false
            if best.species_id == s.id
                s.hasBest = true
            end
        end

#         # saves the best chromo from the current generation
#         if save_best
#             file = open('best_chromo_'+str(self.__generation),'w')
#             pickle.dump(best, file)
#             file.close()
#               end

        #-----------------------------------------
        # Prints chromosome's parents id:  {dad_id, mon_id} -> child_id
#         map(ch-> @printf("{%3d; %3d} -> %3d   Nodes %3d   Connections %3d\n",
#                          ch.parent1_id, ch.parent2_id, ch.id, size(ch)[1], size(ch)[2]), p.population)
        #-----------------------------------------

        # Remove stagnated species and its members (except if it has the best chromosome)
        speciesToKeep = trues(length(p.species))
        deletedSpeciesIds = Int[]
        for i in 1:length(p.species)
            if p.species[i].no_improvement_age > g.cf.max_stagnation
                if !p.species[i].hasBest || p.species[i].no_improvement_age > 2 * g.cf.max_stagnation
                    if report @printf("\n   Species %2d age %2s (with %2d individuals) is stagnated: removing it",
                                      p.species[i].id, p.species[i].age, length(p.species[i])) end
                    speciesToKeep[i] = false
                    push!(deletedSpeciesIds,p.species[i].id)
                end
            end
        end
        p.species = p.species[speciesToKeep] # prune unwanted species

        # remove species' chromosomes from population
        chromosToKeep = trues(length(p.population))
        for i in 1:length(p.population)
            if findfirst(deletedSpeciesIds,p.population[i].species_id) != 0 chromosToKeep[i] = false end
        end
        p.population = p.population[chromosToKeep] # prune unwanted chromosomes

        # Compute spawn levels for each remaining species
        compute_spawn_levels(g, p)

        # Removing species with spawn amount = 0
        speciesToKeep = trues(length(p.species))
        deletedSpeciesIds = Int[]
        for i in 1:length(p.species)

            # This rarely happens
            if p.species[i].spawn_amount == 0
                if report @printf("\n   Species %2d age %2s removed: produced no offspring",p.species[i].id, p.species[i].age) end
                speciesToKeep[i] = false
                push!(deletedSpeciesIds,p.species[i].id)
            end
        end
        p.species = p.species[speciesToKeep] # prune unwanted species

        # remove species' chromosomes from population
        chromosToKeep = trues(length(p.population))
        for i in 1:length(p.population)
            if findfirst(deletedSpeciesIds,p.population[i].species_id) != 0 chromosToKeep[i] = false end
        end
        p.population = p.population[chromosToKeep] # prune unwanted chromosomes

        # Logging speciation stats
        log_species(p)

        if report
            @printf("\nPopulation's average fitness: %3.5f stdev: %3.5f", p.avg_fitness[end], stdeviation(p))
            @printf("\nBest fitness: %2.12s - size: %s - species %s - id %s", best.fitness, size(best), best.species_id, best.id)

            # print some "debugging" information
            @printf("\nSpecies length: %d totalizing %d individuals", length(p.species), sum([length(s) for s in p.species]))
            @printf("\nSpecies ID       : %s",   [s.id for s in p.species])
            @printf("\nEach species size: %s",   [length(s) for s in p.species])
            @printf("\nAmount to spawn  : %s",   [s.spawn_amount for s in p.species])
            @printf("\nSpecies age      : %s",   [s.age for s in p.species])
            @printf("\nSpecies no improv: %s\n", [s.no_improvement_age for s in p.species]) # species no improvement age

            for s in p.species println(s) end
        end

        # Stops the simulation
        if best.fitness > g.cf.max_fitness_threshold
            @printf("Best individual found in epoch %s - complexity: %s\n", p.generation, size(best))
            break
        end

        # -------------------------- Producing new offspring -------------------------- #
        new_population = Chromosome[] # next generation's population

        # Spawning new population
        for s in p.species new_population = vcat(new_population,reproduce(g, s)) end

        # ----------------------------------------------#
        # Controls target population under or overflow  #
        # ----------------------------------------------#
        fill = p.popsize - length(new_population)
        if fill < 0 # overflow
            if report println("\n   Removing $(abs(fill)) excess individual(s) from the new population") end
            # TODO: This is dangerous? I can't remove a species' representant!
            new_population = new_population[1:end+fill] # Removing the last added members
        end

        if fill > 0 # underflow
            if report println("\n   Producing $fill more individual(s) to fill up the new population") end

            # TODO:
            # what about producing new individuals instead of reproducing?
            # increasing diversity from time to time might help
            while fill > 0
                # Selects a random chromosome from population
                parent1 = p.population[rand(1:length(p.population))]
                # Search for a mate within the same species
                found = false
                for parent2 in p.population
                    if parent2.species_id == parent1.species_id && parent2.id != parent1.id
                        child = mutate(crossover(g, parent1, parent2),g)
                        push!(new_population,child)
                        found = true
                        break
                    end
                end
                if !found
                    push!(new_population, mutate(deepcopy(parent1),g)) # will irreversibly mutate parent. ok?
                end # If no mate was found, just mutate it

                fill -= 1
            end
        end

        @assert p.popsize == length(new_population) # Different population sizes!

        # Updates current population
        p.population = new_population

        if time() > t0 + 60 * checkpoint_interval
            create_checkpoint(p,report)
            t0 = time() # updates the counter
        elseif  checkpoint_generation != 0 && p.generation % checkpoint_generation == 0
            create_checkpoint(p,report)
        end
    end
end

function resume_checkpoint(checkpoint)
    # Resumes the simulation from a previous saved point.
    # try:
    #     #file = open(checkpoint)
    #     file = gzip.open(checkpoint)
    # except IOError:
    #     raise
    # print 'Resuming from a previous point: %s' %checkpoint
    # # when unpickling __init__ is not called again
    # previous_pop = pickle.load(file)
    # self.__dict__ = previous_pop.__dict__

    # print 'Loading random state'
    # rstate = pickle.load(file)
    # random.setstate(rstate)
    # #random.jumpahead(1)
    # file.close()
end

function create_checkpoint(p::Population, report)
        # Saves the current simulation state.
        # from time import strftime
        # get current time
        # date = strftime("%Y_%m_%d_%Hh%Mm%Ss")
        # if report print 'Creating checkpoint file at generation: %d' %self.__generation end

        # # dumps 'self'
        # #file = open('checkpoint_'+str(self.__generation), 'w')
        # file = gzip.open('checkpoint_'+str(self.__generation), 'w', compresslevel = 5)
        # # dumps the population
        # pickle.dump(self, file, protocol=2)
        # # dumps the current random state
        # pickle.dump(random.getstate(), file, protocol=2)
        # file.close()
end

# if __name__ ==  '__main__' :

#     # sample fitness function
#     def eval_fitness(population):
#         for individual in population:
#             individual.fitness = 1.0

#     # set fitness function
#     Population.evaluate = eval_fitness

#     # creates the population
#     pop = Population()
#     # runs the simulation for 250 epochs
#     pop.epoch(250)


text:
include("DeclarativePackages.jl")

Pkg.init()
Pkg.add("Compat")
using DeclarativePackages, Compat

if VERSION < v"0.5.0"
    readstring = readall
end

if !haskey(ENV, "DECLARE_VERBOSITY")
    ENV["DECLARE_VERBOSITY"] = 1
end

function installpackages()
    lines = readfile()
    init(lines)
    packages = parselines(lines)
    needbuilding = install(packages)
    resolve(packages, needbuilding)
    finish()
end

function readfile()
    log(1, "Parsing $(ENV["DECLARE"]) ... ")
    lines = split(readstring(ENV["DECLARE"]), '\n')
    lines = map(x->replace(x, r"#.*", ""), lines)
    lines = filter(x->!isempty(x), lines)
    return lines
end

pkgpath(basepath, pkg) = normpath(basepath*"/v$(VERSION.major).$(VERSION.minor)/$pkg/")
markreadonly(path) = run(`chmod a-w $path`)
stepout(path, n=1) = normpath(path*"/"*repeat("../",n))

function hardlinkdirs(existingpath, path)
    log(3, "hardlinking: existingpath: $existingpath\npath: $path")
    assert(existingpath[end]=='/')
    assert(path[end]=='/')
    mkpath(path)
    readdirabs(path) = map(x->(x, path*x), readdir(path))
    items = readdirabs(existingpath)
    for dir in filter(x->isdir(x[2]), items)
        hardlinkdirs(dir[2]*"/", path*dir[1]*"/")
    end
    for file in filter(x->!isdir(x[2]), items)
        hardlinkfile(file[2] , path*file[1])
    end
end

function hardlinkfile(from, to)
    @static if is_apple() ccall((:link, "libc"), Int, (Ptr{UInt8}, Ptr{UInt8}), from, to) end
    @static if is_linux() ccall((:link, "libc.so.6"), Int, (Ptr{UInt8}, Ptr{UInt8}), from, to) end
end

gitcmd(path, cmd) = `git --git-dir=$path.git --work-tree=$path $(split(cmd))`
function gitcommitof(path)
    log(2, "gitcommitof $path")
    cmd = gitcmd(path, "log -n 1 --format=%H")
    log(2, "gitcommitof cmd $cmd")
    r = try
        strip(readstring(cmd))
    catch
        ""
    end
    log(2, "gitcommitof result $r")
    r
end

function gitcommitoftag(path, tag)
    contains(path, "METADATA") && return ""
    length(tag) > 1 && tag[1] != 'v' && return ""
    cmd = gitcmd(path, "rev-list -n 1 $tag")
    strip(readstring(cmd))
end

function gitclone(name, url, path, commit="")
    log(2, "gitclone: name: $name url: $url path: $path commit: $commit")
    run(`git clone $url $path`)
    if isempty(commit)
        commit = gitcommitof(path)
    else
        # check if the repo knows this commit. if not, check in METADATA
        isknown = ismatch(Regex(commit), readstring(gitcmd(path, "tag")))
        if !isknown
            filename = Pkg.dir("METADATA/$name/versions/$(commit[2:end])/sha1")
            if exists(filename)
                commit = strip(readstring(filename))
            else
                if commit[1] == 'v'
                    error("gitclone: Could not find a commit hash for version $commit for package $name ($url)")
                end
            end
        end
    end

    run(gitcmd(path, "checkout --force -b pinned.$commit.tmp $commit"))
end

function existscheckout(pkg, commit)
    log(2, "######## existscheckout $pkg $commit")
    basepath = stepout(Pkg.dir(), 2)
    dirs = readdir(basepath)
    nontmp = filter(x->length(x)>3 && x[1:4]!="tmp_", dirs)
    log(2, "  nontmp dirs: $nontmp")
    for dir in nontmp
        path = pkgpath(basepath*dir, pkg)
        !exists(path) && continue
        existingcommit = gitcommitof(path)
        existingtagcommit = try gitcommitoftag(path, commit) catch "" end
        log(2, "  existinging commit / existingtagcommit / wanted commit:  $existingcommit / $existingtagcommit / $commit")
        if exists(path) && (existingcommit == commit || existingcommit == existingtagcommit)
            log(2, "existscheckout: found $path for $pkg@$commit")
            return path
        end
    end
    return ""
end

function init(lines)
    tmpdir = "tmp_$(randstring(32))"
    ENV["JULIA_PKGDIR"] = normpath(Pkg.dir()*"/../../$tmpdir")
    metadata = filter(x->ismatch(r"METADATA.jl", x), lines)
    commit = ""
    if length(metadata)>0
        assert(length(metadata)==1)
        m = split(metadata[1])
        url = split(metadata[1])[1]
        length(m) > 1 ? commit = m[2] : ""
        log(2, "Found URL $url$(isempty(commit) ? "" : "@$commit") for METADATA")
    else
        url = "https://github.com/JuliaLang/METADATA.jl.git"
    end
    mkpath(Pkg.dir())
    path = Pkg.dir("METADATA/")
    installorlink("METADATA", url, path, commit)
    markreadonly(Pkg.dir("METADATA"))
end

parselines(lines) = filter(x->isa(x,Package), map(parseline, lines))
function parseline(a)
    parts = split(strip(a))

    if parts[1][1] == '@'
        os = parts[1]
        shift!(parts)
    else
        os = ""
    end

    nameorurl = parts[1]
    if contains(nameorurl, "/")
        url = nameorurl
        name = replace(replace(split(url, "/")[end], ".git", ""), ".jl", "")
        isregistered = false
    else
        name = nameorurl
        url = strip(readstring("$(Pkg.dir())/METADATA/$name/url"))
        isregistered = true
    end
    if name=="METADATA"
        return []
    end

    commit = length(parts)>1 ? parts[2] : (isregistered ? "METADATA" : "")
    if length(split(commit,"."))==3
        commit = "v"*commit
    end
    return Package(os, name, url, commit, isregistered)
end

type Package
    os
    name
    url
    commit
    isregistered
end

function install(packages::Array)
    osx = filter(x->x.os=="@osx", packages)
    unix = filter(x->x.os=="@unix", packages)
    linux = filter(x->x.os=="@linux", packages)
    windows = filter(x->x.os=="@windows", packages)
    everywhere = filter(x->x.os=="", packages)
    @static if is_apple() map(install, osx) end
    @static if is_unix() map(install, unix) end
    @static if is_linux() map(install, linux) end
    @static if is_windows() map(install, windows) end
    needbuilding = filter(x->x!=nothing, map(install, everywhere))
end

function installorlink(name, url, path, commit)
    log(2, "Installorlink: $name $url $commit $path")
    existingpath = existscheckout(name, commit)
    if isempty(existingpath)
        gitclone(name, url, path, commit)
        return name
    else
        log(1, "Linking $(name) ...")
        hardlinkdirs(existingpath, path)

        # link the compiled module, too
        # @show name url path commit existingpath path
        # if name != "METADATA"
            # v = "v$(VERSION.major).$(VERSION.minor)"
            # from = joinpath("/",split(existingpath,"/")[1:end-3]..., "lib", v, name*".ji")
            # todir = joinpath("/",split(path        ,"/")[1:end-3]..., "lib", v)
            # mkpath(todir)
            # to =  joinpath(todir, name*".ji")
            # @show from to
            # hardlinkfile(from, to)
        # end
        return
    end
end

function install(a::Package)
    path = Pkg.dir(a.name*"/")

    version(a) = VersionNumber(map(x->parse(Int,x), split(a, "."))...)
    function latest()
        versionsdir = Pkg.dir("METADATA/$(a.name)/versions/")
        if exists(versionsdir)
            "v"*string(maximum(map(version, readdir(versionsdir))))
        else
            ""
        end
    end
    metadatacommit(version) = strip(readstring(Pkg.dir("METADATA/$(a.name)/versions/$(version[2:end])/sha1")))

    commit = a.commit == "METADATA" ? latest() : a.commit
    installorlink(a.name, a.url, path, commit)
end

function resolve(packages, needbuilding)
    requirename = Pkg.dir()*"/REQUIRE"
    log(3, requirename)
    open(requirename,"w") do io
        for pkg in packages
            if !isempty(pkg.commit) && pkg.commit[1]=='v'
                m,n,o = map(x->parse(Int,x), split(split(pkg.commit[2:end],'~')[1], '.'))
                versions = "$m.$n.$o $m.$n.$(o+1)-"
            else
                versions = ""
            end
            log(3, "writing REQUIRE: $(pkg.os) $(pkg.name) $versions\n")
            write(io, "$(pkg.os) $(pkg.name) $versions\n")

            # add test dependencies
            if haskey(ENV, "DECLARE_INCLUDETEST") && ENV["DECLARE_INCLUDETEST"]=="true"
                testrequire = Pkg.dir(pkg.name*"/test/REQUIRE")
                if exists(testrequire)
                    write(io, readstring(testrequire))
                end
            end
        end
    end
    log(1, "Invoking Pkg.resolve() ...")
    Pkg.resolve()
    map(buildifnecessary, needbuilding)
end

function buildifnecessary(x)
    depsdir = Pkg.dir(x,"deps")
    buildscript = joinpath(depsdir, "build.jl")
    declarebuilt = joinpath(depsdir, "built.by.declarepackages.jl")
    exists(buildscript) || return
    exists(declarebuilt) && return
    Pkg.build(x)
    touch(declarebuilt)
end

function finish()
    exportDECLARE(ENV["DECLARE"])

    @static if is_apple() md5 = strip(readstring(`md5 -q $(ENV["DECLARE"])`)) end
    @static if is_linux() md5 = strip(readstring(`md5sum $(ENV["DECLARE"])`)) end
    md5 = split(md5)[1]
    if haskey(ENV, "DECLARE_INCLUDETEST") && ENV["DECLARE_INCLUDETEST"]=="true"
        md5 = md5*"withtest"
    end
    dir = normpath(Pkg.dir()*"/../../"*md5*"-$(VERSION.major).$(VERSION.minor)")

    if exists(dir)
        run(`chmod -R a+w $dir`)
        rm(dir; recursive=true)
    end
    mv(stepout(Pkg.dir(),1), dir)
    symlink(dir, stepout(Pkg.dir())[1:end-1])
    ENV["JULIA_PKGDIR"] = dir

    cp(ENV["DECLARE"], joinpath(dir,"DECLARE"))

    log(1, "Marking $dir read-only ...")
    run(pipeline(`find $dir -maxdepth 1`,`xargs chmod 555 `))
    run(`chmod 755 $dir`)

    log(1, "Finished installing packages for $(ENV["DECLARE"]).")
end

installpackages()


text:
## Day 18: Snailfish ##
#######################

# This is a rewrite, this time using a better data structure & functional programming

abstract type Tree end

struct Pair <: Tree
    left::Tree
    right::Tree
end

struct Leaf <: Tree
    value::Int
end

Pair(l::Int, r::Int) = Pair(Leaf(l), Leaf(r))

parse(data::Int)::Tree = Leaf(data)
parse(data::Array)::Tree = Pair(parse(data[1]), parse(data[2]))

# For debugging
s(pair::Pair)::String = "{$(s(pair.left)),$(s(pair.right))}"
s(leaf::Leaf) = "$(leaf.value)"

function split_one(leaf::Leaf)::Tuple{Tree,Symbol}
    if leaf.value >= 10
        (Pair(Int(floor(leaf.value / 2)), Int(ceil(leaf.value / 2))), :split)
    else
        (leaf, :no_split)
    end
end

function split_one(pair::Pair)::Tuple{Tree,Symbol}
    new_l, action = split_one(pair.left)
    if action == :split
        (Pair(new_l, pair.right), :split)
    else
        new_r, action = split_one(pair.right)
        (Pair(pair.left, new_r), action)
    end
end

@assert (parse([[[[0, 7], 4], [[7, 8], [0, 13]]], [1, 1]]), :split) == split_one(parse([[[[0, 7], 4], [15, [0, 13]]], [1, 1]]))
@assert (parse([[[[0, 7], 4], [[7, 8], [0, [6, 7]]]], [1, 1]]), :split) == split_one(parse([[[[0, 7], 4], [[7, 8], [0, 13]]], [1, 1]]))
@assert (parse([[[[0, 7], 4], [[7, 8], [0, [6, 7]]]], [1, 1]]), :no_split) == split_one(parse([[[[0, 7], 4], [[7, 8], [0, [6, 7]]]], [1, 1]]))

OInt = Union{Int,Nothing}

left_add(tree::Tree, _::Nothing) = tree
left_add(pair::Pair, value::Int) = Pair(left_add(pair.left, value), pair.right)
left_add(leaf::Leaf, value::Int) = Leaf(leaf.value + value)

right_add(tree::Tree, _::Nothing) = tree
right_add(pair::Pair, value::Int) = Pair(pair.left, right_add(pair.right, value))
right_add(leaf::Leaf, value::Int) = Leaf(leaf.value + value)

@assert parse([[6, 2], 6]) == left_add(parse([[1, 2], 6]), 5)
@assert parse([[1, 2], 11]) == right_add(parse([[1, 2], 6]), 5)

ExplodeData = Tuple{OInt,Tree,OInt,Symbol}

function explode_one(tree::Tree)::Tuple{Tree,Symbol}
    _, new_tree, _, action = explode_one(tree, 1)
    (new_tree, action)
end
explode_one(leaf::Leaf, _::Int)::ExplodeData = (nothing, leaf, nothing, :no_explode)
function explode_one(pair::Pair, layer::Int)::ExplodeData
    if layer > 4
        # Found a deeply nested pair, explode it
        # First, assert that there is nothing on layer 6
        @assert typeof(pair.left) == Leaf
        @assert typeof(pair.left) == Leaf
        (pair.left.value, Leaf(0), pair.right.value, :explode)
    else
        # Check if any of the children explodes. If they do, carry up the OInt
        l_int, new_l, r_int, action = explode_one(pair.left, layer + 1)
        if action == :explode
            (l_int, Pair(new_l, left_add(pair.right, r_int)), nothing, :explode)
        else
            l_int, new_r, r_int, action = explode_one(pair.right, layer + 1)
            (nothing, Pair(right_add(pair.left, l_int), new_r), r_int, action)
        end
    end
end

@assert (parse([[[[0, 9], 2], 3], 4]), :explode) == explode_one(parse([[[[[9, 8], 1], 2], 3], 4]))
@assert (parse([7, [6, [5, [7, 0]]]]), :explode) == explode_one(parse([7, [6, [5, [4, [3, 2]]]]]))
@assert (parse([[6, [5, [7, 0]]], 3]), :explode) == explode_one(parse([[6, [5, [4, [3, 2]]]], 1]))
@assert (parse([[3, [2, [8, 0]]], [9, [5, [4, [3, 2]]]]]), :explode) == explode_one(parse([[3, [2, [1, [7, 3]]]], [6, [5, [4, [3, 2]]]]]))
@assert (parse([[3, [2, [8, 0]]], [9, [5, [7, 0]]]]), :explode) == explode_one(parse([[3, [2, [8, 0]]], [9, [5, [4, [3, 2]]]]]))

function normalize(tree::Tree)::Tree
    new_tree, action = explode_one(tree)
    if action == :explode
        normalize(new_tree)
    else
        new_tree, action = split_one(new_tree)
        if action == :split
            normalize(new_tree)
        else
            new_tree
        end
    end
end

@assert parse([[[[0, 7], 4], [[7, 8], [6, 0]]], [8, 1]]) == normalize(parse([[[[[4, 3], 4], 4], [7, [[8, 4], 9]]], [1, 1]]))

add(t1::Tree, t2::Tree)::Tree = normalize(Pair(t1, t2))
sum(trees::Vector)::Tree = foldl(add, trees)

@assert sum(parse.([[[[0, [4, 5]], [0, 0]], [[[4, 5], [2, 6]], [9, 5]]],
    [7, [[[3, 7], [4, 3]], [[6, 3], [8, 8]]]],
    [[2, [[0, 8], [3, 4]]], [[[6, 7], 1], [7, [1, 6]]]],
    [[[[2, 4], 7], [6, [0, 5]]], [[[6, 8], [2, 8]], [[2, 1], [4, 5]]]],
    [7, [5, [[3, 8], [1, 4]]]],
    [[2, [2, 2]], [8, [8, 1]]],
    [2, 9],
    [1, [[[9, 3], 9], [[9, 0], [0, 7]]]],
    [[[5, [7, 4]], 7], 1],
    [[[[4, 2], 2], 6], [8, 7]]])) == parse([[[[8, 7], [7, 7]], [[8, 6], [7, 7]]], [[[0, 7], [6, 6]], [8, 7]]])

include("test-18.jl")
include("input-18.jl")

magnitude(leaf::Leaf)::Int = leaf.value
magnitude(pair::Pair)::Int = 3 * magnitude(pair.left) + 2 * magnitude(pair.right)

@assert 4140 == magnitude(sum(parse.(day18_test)))
@assert 4469 == magnitude(sum(parse.(day18_input)))

function max_magnitude(list_of_trees::Vector)::Int
    best = 0
    for i = 1:length(list_of_trees), j = 1:length(list_of_trees)
        if i == j
            continue
        end
        best = max(best, magnitude(add(list_of_trees[i], list_of_trees[j])))
    end
    best
end

@assert 3993 == max_magnitude(parse.(day18_test))
@assert 4770 == max_magnitude(parse.(day18_input))


text:
##############

@views function ThermalAssembly( Cp, rho, kx, ky, dx, dy, dt, steady, Th_BC, ncx, ncy )

    Num       = reshape(1:ncx*ncy,ncx,ncy)
    iC        = Num
    iW        =  ones(size(Num)); iW[2:end-0,:] = Num[1:end-1,:]
    iE        =  ones(size(Num)); iE[1:end-1,:] = Num[2:end-0,:]
    iS        =  ones(size(Num)); iS[:,2:end-0] = Num[:,1:end-1]
    iN        =  ones(size(Num)); iN[:,1:end-1] = Num[:,2:end-0]
    kW        = zeros(size(Num)); kW[:,:] = kx[1:end-1,:]
    kE        = zeros(size(Num)); kE[:,:] = kx[2:end-0,:]
    kS        = zeros(size(Num)); kS[:,:] = ky[:,1:end-1]
    kN        = zeros(size(Num)); kN[:,:] = ky[:,2:end-0]
    BCW       = zeros(size(Num)); BCW[1,:]   .= Th_BC.type_W
    BCE       = zeros(size(Num)); BCE[end,:] .= Th_BC.type_E
    BCS       = zeros(size(Num)); BCS[:,1]   .= Th_BC.type_S
    BCN       = zeros(size(Num)); BCN[:,end] .= Th_BC.type_S
    # Periodic stencil
    if Th_BC.type_W==0 && Th_BC.type_E==0
        iW[  1,:] = Num[end,:]
        iE[end,:] = Num[  1,:]
    end
    # Matrix coeffcients
    fW     = 1.0 .+ 1.0.*(BCW.==1) .- 1.0.*(BCW.==2)
    fE     = 1.0 .+ 1.0.*(BCE.==1) .- 1.0.*(BCE.==2)
    fS     = 1.0 .+ 1.0.*(BCS.==1) .- 1.0.*(BCS.==2)
    fN     = 1.0 .+ 1.0.*(BCN.==1) .- 1.0.*(BCN.==2)
    cC     =  (steady!=1) .* Cp.*rho./dt .+ fN.*kN./dy^2 .+ fS.*kS./dy^2 .+ fE.*kE./dx^2 .+ fW.*kW./dx^2
    cW     = -kW./dx^2; cW[BCW.!=0] .= 0.0
    cE     = -kE./dx^2; cE[BCE.!=0] .= 0.0
    cS     = -kS./dy^2; cS[BCS.!=0] .= 0.0
    cN     = -kN./dy^2; cN[BCN.!=0] .= 0.0

    # # Matrix coeffcients
    # cC        = zeros(size(Num))
    # cW        = zeros(size(Num))
    # cE        = zeros(size(Num))
    # cS        = zeros(size(Num))
    # cN        = zeros(size(Num))
    # @threads for i = 1:ncx # vectorisation does not work
    #     for j = 1:ncy
    #         # Setup BCs
    #         DirW    = (i==1) & (Th_BC.type_W==1)
    #         NeuW    = (i==1) & (Th_BC.type_W==2)
    #         fW      = 1.0 + 1.0*DirW - 1.0*NeuW
    #         DirE    = (i==ncx) & (Th_BC.type_E==1)
    #         NeuE    = (i==ncx) & (Th_BC.type_E==2)
    #         fE      = 1.0 + 1.0*DirE - 1.0*NeuE
    #         DirS    = (j==1) & (Th_BC.type_S==1)
    #         NeuS    = (j==1) & (Th_BC.type_S==2)
    #         fS      = 1.0 + 1.0*DirS - 1.0*NeuS
    #         DirN    = (j==ncy) & (Th_BC.type_N==1)
    #         NeuN    = (j==ncy) & (Th_BC.type_N==2)
    #         fN      = 1.0 + 1.0*DirN - 1.0*NeuN
    #         # Coefficients
    #         cC[i,j] =  (steady!=1) * Cp[i,j]*rho[i,j]/dt + fN*kN[i,j]/dy^2 + fS*kS[i,j]/dy^2 + fE*kE[i,j]/dx^2 + fW*kW[i,j]/dx^2
    #         cW[i,j] = -(i>1  ) * kW[i,j]/dx^2
    #         cE[i,j] = -(i<ncx) * kE[i,j]/dx^2
    #         cS[i,j] = -(j>  1) * kS[i,j]/dy^2
    #         cN[i,j] = -(j<ncy) * kN[i,j]/dy^2
    #     end
    # end

    # Sparse matrix
    I = [iC[:]; iC[:]; iC[:]; iC[:]; iC[:]]
    J = [iC[:]; iW[:]; iE[:]; iS[:]; iN[:]]
    V = [cC[:]; cW[:]; cE[:]; cS[:]; cN[:]]
    K = sparse( I, J, V)
    droptol!(K, 1e-8)
    # display(UnicodePlots.spy(K))
    return K
end
export ThermalAssembly

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

@views function ThermalSolve!( T, K, F, ncx, ncy )
    # Solve
    b     = zeros(ncx*ncy,1); b[:,1] = F[:]
    Kchol = cholesky(K)
    dT    = Kchol\b
    T   .-= reshape(dT,ncx,ncy)
end
export ThermalSolve!

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

function ThermalResidual(T, T0, rho, Cp, H, kx, ky, dx, dy, dt, steady, Th_BC, ncx, ncy)
    """ Thermal residual """
    Tex = zeros(ncx+2, ncy+2)
    Tex[2:end-1,2:end-1] .= T
    Tex[  1,2:end-1]     .= (Th_BC.type_W==0) * T[end,:] .+ (Th_BC.type_W==2) * T[  1,:] .+ (Th_BC.type_W==1) * (2.0*Th_BC.Dir_W[:] - T[  1,:])
    Tex[end,2:end-1]     .= (Th_BC.type_E==0) * T[  1,:] .+ (Th_BC.type_E==2) * T[end,:] .+ (Th_BC.type_E==1) * (2.0*Th_BC.Dir_E[:] - T[end,:])
    Tex[2:end-1,  1]     .=                                 (Th_BC.type_S==2) * T[:,  1] .+ (Th_BC.type_S==1) * (2.0*Th_BC.Dir_S[:] - T[:,  1])
    Tex[2:end-1,end]     .=                                 (Th_BC.type_N==2) * T[:,end] .+ (Th_BC.type_N==1) * (2.0*Th_BC.Dir_N[:] - T[:,end])

    # Compute residual
    F = zeros( ncx, ncy )
    @tturbo for i = 1:ncx
        for j = 1:ncy
            qW = -kx[i  ,j]*(Tex[i+1,j+1] - Tex[i+0,j+1])/dx
            qE = -kx[i+1,j]*(Tex[i+2,j+1] - Tex[i+1,j+1])/dx
            qS = -ky[i,j  ]*(Tex[i+1,j+1] - Tex[i+1,j+0])/dy
            qN = -ky[i,j+1]*(Tex[i+1,j+2] - Tex[i+1,j+1])/dy
            F[i,j] = (steady!=1) * rho[i,j]*Cp[i,j] * (T[i,j] - T0[i,j])/dt + (qE-qW)/dx + (qN-qS)/dy + H[i,j]
        end
    end
return F
end
export ThermalResidual

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