text:
# Setup everything using the common setup.
include("CloudGraphSetup.jl")

# And check that if we encode and decode this type, it's exactly the same.
# Make a packed data test structure.
fullType = DataTest(rand(10,10), "This is a test string", rand(Int32,10,10));
typePackedRegName = string(PackedDataTest);
typeOriginalRegName = string(DataTest);
# Now lets encode and decode to see.
println("Encoding...")
testPackedType = convert(PackedDataTest, fullType);
# testPackedType = cloudGraph.packedOriginalDataTypes[typeOriginalRegName].encodingFunction(PackedDataTest, fullType);

println("Decoding...")
testFullType = convert(DataTest, testPackedType);
# testFullType = cloudGraph.packedPackedDataTypes[typePackedRegName].decodingFunction(DataTest, testPackedType);

@test json(testFullType) == json(fullType)
println("Success!")

# Creating a local test graph.
print("[TEST] Creating a CloudVertex from an ExVertex...")
localGraph = graph(ExVertex[], ExEdge{ExVertex}[]);
#Make an ExVertex that may be encoded
v = make_vertex(localGraph, "TestVertex");
vertex = Graphs.add_vertex!(localGraph, v);
vertex.attributes["data"] = fullType;
vertex.attributes["age"] = 64;
vertex.attributes["latestEstimate"] = [0.0,0.0,0.0];
bigData = CloudGraphs.BigData();
testElementLegacy = CloudGraphs.BigDataElement("TestLegacy", "Performance test dataset legacy.", rand(UInt8,100), -1); #Data element
testElementDict = CloudGraphs.BigDataElement("TestDictSet", "Performance test dataset new dict type.", Dict{String, Any}("testString"=>"Test String", "randUint8"=>rand(UInt8,100)), -1); #Data element
# TODO: Check BigData scenarios.
# append!(bigData.dataElements, [testElementLegacy, testElementDict]);
vertex.attributes["bigData"] = bigData;
# Now encoding the structure to CloudGraphs vertex
cloudVertex = CloudGraphs.exVertex2CloudVertex(vertex);
println("Success!");

print("[TEST] Adding a vertex...")
CloudGraphs.add_vertex!(cloudGraph, cloudVertex);
println("Success!")

# Get the node from Neo4j.
print("[TEST] Retrieving a node from CloudGraph...")
cloudVertexRet = CloudGraphs.get_vertex(cloudGraph, cloudVertex.neo4jNode.id, false) # fullType not required
# Check that all the important bits match using string comparisons of the JSON form of the structures
@test json(cloudVertex.packed) == json(cloudVertexRet.packed)
#@test json(cloudVertex.bigData) == json(cloudVertexRet.bigData)
@show "Expected = ", json(cloudVertex.properties)
@show "Received = ", json(cloudVertexRet.properties)
@test json(cloudVertex.properties) == json(cloudVertexRet.properties)
@test cloudVertex.neo4jNodeId == cloudVertexRet.neo4jNodeId
@test cloudVertexRet.neo4jNode != Nothing
println("Success!")

print("[TEST] Testing the update method...")
cloudVertex.properties["age"] = 100;
cloudVertex.properties["latestEstimate"] = [5.0, 5.0, 5.0];
CloudGraphs.update_vertex!(cloudGraph, cloudVertex, false);
# Let's retrieve it and see if it is updated.
cloudVertexRet = CloudGraphs.get_vertex(cloudGraph, cloudVertex.neo4jNode.id, false)
# And check that it matches
@show json(cloudVertexRet.properties);
@test json(cloudVertex.properties) == json(cloudVertexRet.properties);
println("Success!")

# Label testing
println("[Test] Creating a vertex with labels...")
cloudVertexWithLabels = deepcopy(cloudVertex);
cloudVertexWithLabels.labels = ["LABEL1", "LABEL2"];
CloudGraphs.add_vertex!(cloudGraph, cloudVertexWithLabels);
println("Success!")
println("[Test] Getting the vertex and checking labels exist...")
cloudVertexWithLabelsRet = CloudGraphs.get_vertex(cloudGraph, cloudVertexWithLabels.neo4jNode.id, false)
@test length(setdiff(cloudVertexWithLabels.labels, cloudVertexWithLabelsRet.labels)) == 0
println("[Test] Adding a label...")
push!(cloudVertexWithLabels.labels, "AnotherLabel")
CloudGraphs.update_vertex!(cloudGraph, cloudVertexWithLabels, false);
cloudVertexWithLabelsRet = CloudGraphs.get_vertex(cloudGraph, cloudVertexWithLabels.neo4jNode.id, false)
@test length(setdiff(cloudVertexWithLabels.labels, cloudVertexWithLabelsRet.labels)) == 0
println("Success!")
# Now clear out all the labels
println("[Test] Clearing all labels...")
cloudVertexWithLabels.labels = Vector{AbstractString}()
CloudGraphs.update_vertex!(cloudGraph, cloudVertexWithLabels, false);
cloudVertexWithLabelsRet = CloudGraphs.get_vertex(cloudGraph, cloudVertexWithLabels.neo4jNode.id, false)
@test length(setdiff(cloudVertexWithLabels.labels, cloudVertexWithLabelsRet.labels)) == 0
println("Success!")

print("[TEST] Deleting a CloudGraph vertex...")
CloudGraphs.delete_vertex!(cloudGraph, cloudVertex);
CloudGraphs.delete_vertex!(cloudGraph, cloudVertexWithLabels);
println("Success!")

print("[TEST] Negative testing for double deletions...")
# Testing a double-delete
try
  CloudGraphs.delete_vertex!(cloudGraph, cloudVertex);
catch
  print("Success!")
end
# Testing the deletion of an apparently existing node
@test_throws ErrorException CloudGraphs.delete_vertex!(cloudGraph, cloudVertexRet);
print("Success!")

print("[TEST] Making an edge...")
# Create two vertices
cloudVert1 = deepcopy(cloudVertex);
cloudVert1.properties["label"] = "Sam's Vertex 1";
cloudVert2 = deepcopy(cloudVertex);
cloudVert2.properties["label"] = "Sam's Vertex 2";
cloudVert3 = deepcopy(cloudVertex);
cloudVert3.properties["label"] = "Sam's Vertex 3";
CloudGraphs.add_vertex!(cloudGraph, cloudVert1);
CloudGraphs.add_vertex!(cloudGraph, cloudVert2);
CloudGraphs.add_vertex!(cloudGraph, cloudVert3);

# Create an edge and add it to the graph.
# Test props
props = Dict{AbstractString, Any}(string("Test") => 8); #UTF8String   utf8(..)
edge12 = CloudGraphs.CloudEdge(cloudVert1, cloudVert2, "DEPENDENCE");
CloudGraphs.add_edge!(cloudGraph, edge12);
edge23 = CloudGraphs.CloudEdge(cloudVert2, cloudVert3, "DEPENDENCE");
CloudGraphs.add_edge!(cloudGraph, edge23);
edge31 = CloudGraphs.CloudEdge(cloudVert3, cloudVert1, "DEPENDENCE");
CloudGraphs.add_edge!(cloudGraph, edge31);
@test edge12.neo4jEdge != nothing
@test edge12.neo4jEdgeId != -1
@test edge23.neo4jEdge != nothing
@test edge23.neo4jEdgeId != -1
@test edge31.neo4jEdge != nothing
@test edge31.neo4jEdgeId != -1
println("Success!")

print("[TEST] Get edge from graph...")
gotedge = CloudGraphs.get_edge(cloudGraph, edge12.neo4jEdgeId)
@test typeof(gotedge) == CloudGraphs.CloudEdge
@test edge12.neo4jEdgeId == gotedge.neo4jEdgeId
@test gotedge.neo4jEdge != edge12.neo4jEdge
@test edge12.edgeType == gotedge.edgeType
@test edge12.neo4jSourceVertexId == gotedge.neo4jSourceVertexId
@test edge12.neo4jDestVertexId == gotedge.neo4jDestVertexId
@test edge12.properties == gotedge.properties

function testCloudGraphsNodeCompares(a::Neo4j.Node, b::Neo4j.Node)
  @test a.paged_traverse == b.paged_traverse
  @test a.labels == b.labels
  @test a.outgoing_relationships == b.outgoing_relationships
  @test a.traverse == b.traverse
  @test a.all_typed_relationships == b.all_typed_relationships
  @test a.all_relationships == b.all_relationships
  @test a.property == b.property
  @test a.self == b.self
  @test a.outgoing_typed_relationships == b.outgoing_typed_relationships
  @test a.properties == b.properties
  @test a.incoming_relationships == b.incoming_relationships
  @test a.incoming_typed_relationships == b.incoming_typed_relationships
  @test a.id == b.id
  # ignore packed data
  # @test a.data == b.data
  nothing
end

testCloudGraphsNodeCompares(edge12.DestVertex.neo4jNode, gotedge.DestVertex.neo4jNode)
testCloudGraphsNodeCompares(edge12.SourceVertex.neo4jNode, gotedge.SourceVertex.neo4jNode)

@test json(edge12.SourceVertex.packed) == json(gotedge.SourceVertex.packed)
@test edge12.SourceVertex.properties == gotedge.SourceVertex.properties
  #@test  edge.SourceVertex.bigData == gotedge.SourceVertex.bigData
@test   edge12.SourceVertex.neo4jNodeId == gotedge.SourceVertex.neo4jNodeId
# @test   edge12.SourceVertex.neo4jNode == gotedge.SourceVertex.neo4jNode
@test   edge12.SourceVertex.isValidNeoNodeId == gotedge.SourceVertex.isValidNeoNodeId
# @test   edge12.SourceVertex.exVertexId == gotedge.SourceVertex.exVertexId
@test   edge12.SourceVertex.isValidExVertex == gotedge.SourceVertex.isValidExVertex
println("Success!")

print("[TEST] Finding all neighbors of a vertex...")
neighs1 = CloudGraphs.get_neighbors(cloudGraph, cloudVert1)
neighs1in = CloudGraphs.get_neighbors(cloudGraph, cloudVert1, incoming=true, outgoing=false)
neighs1out = CloudGraphs.get_neighbors(cloudGraph, cloudVert1, incoming=false, outgoing=true)

@show neighs1
@test length(neighs1) == 2
@test length(neighs1in) == 1
@test length(neighs1out) == 1
@test neighs1[1].neo4jNodeId == cloudVert3.neo4jNodeId || neighs1[2].neo4jNodeId == cloudVert3.neo4jNodeId
@test neighs1[1].neo4jNodeId == cloudVert2.neo4jNodeId || neighs1[2].neo4jNodeId == cloudVert2.neo4jNodeId
@test neighs1in[1].neo4jNodeId == cloudVert3.neo4jNodeId
@test neighs1out[1].neo4jNodeId == cloudVert2.neo4jNodeId
println("Success!")

print("[TEST] Deleting an edge...")
CloudGraphs.delete_edge!(cloudGraph, edge12)
CloudGraphs.delete_edge!(cloudGraph, edge23)
CloudGraphs.delete_edge!(cloudGraph, edge31)
neighs = CloudGraphs.get_neighbors(cloudGraph, cloudVert1)
@test length(neighs) == 0
@test edge12.neo4jEdgeId == -1
@test edge12.neo4jEdge == nothing
println("Success!")

print("[TEST] Deleting nodes...")
CloudGraphs.delete_vertex!(cloudGraph, cloudVert1)
@test cloudVert1.neo4jNode == nothing
@test cloudVert1.neo4jNodeId == -1
CloudGraphs.delete_vertex!(cloudGraph, cloudVert2)
@test cloudVert2.neo4jNode == nothing
@test cloudVert2.neo4jNodeId == -1
CloudGraphs.delete_vertex!(cloudGraph, cloudVert3)
@test cloudVert3.neo4jNode == nothing
@test cloudVert3.neo4jNodeId == -1
println("Success!")


text:
function nlp_cvx_205_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - linear objective
    # - intersection convex quadratic constraints
    # - exponential cones
    # Variants
    #   010 - intersection set

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y >= 0, start = 0.1)
    @variable(model, z)

    @objective(model, Max, y)
    @NLconstraint(model, y * exp(x / y) <= z)
    @NLconstraint(model, y * exp(-x / y) <= z)
    @NLconstraint(model, x^2 + y^2 <= -z + 5)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 1.7912878443121907, tol = objective_tol)
    return check_solution(
        [x, y, z],
        [0.0, 1.7912878443121907, 1.7912878443121907],
        tol = primal_tol,
    )
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

struct GPUmeta
  dev
  ctx
  queue
  mult_kernel
end

function gpu_setup()
  @assert cl ≠ nothing "OpenCL.jl not installed, cannot use GPU"
  @assert clfft ≠ nothing "CLFFT.jl not installed, cannot use GPU"

  devs = cl.devices(:gpu)
  if isempty(devs)
    @warn "GPU not found, falling back to other OpenCL devices"
    devs = cl.devices()
  end
  @assert !isempty(devs) "OpenCL device not found, make sure drivers are installed"

  dev = []
  devnames = map(d -> d[:platform][:name], devs)
  for vendor in ["NVIDIA","AMD","Intel"], (idx,name) in enumerate(devnames)
    if contains(name, vendor)
      dev = devs[idx]
      break
    end
  end

  devtype = uppercase(string(dev[:device_type]))
  devname = dev[:name]

  info("using $devtype $devname")

  ctx = cl.Context(dev)
  queue = cl.CmdQueue(ctx)
  mult_kernel = basic_kernels(ctx)

  GPUmeta(dev, ctx, queue, mult_kernel)
end

function basic_kernels(ctx)
  const mult_kernel = "
    __kernel void mult(__global const float2 *a,
                       __global const float2 *b,
                       __global float2 *c)
    {
      int gid = get_global_id(0);
      c[gid].x = a[gid].x*b[gid].x - a[gid].y*b[gid].y;
      c[gid].y = a[gid].x*b[gid].y + a[gid].y*b[gid].x;
    }
  "

  # build OpenCL program
  prog = cl.Program(ctx, source=mult_kernel) |> cl.build!

  cl.Kernel(prog, "mult")
end

function clfftpad(A::AbstractArray)
  # clFFT releases support powers of 2, 3, 5, ...
  radices = [2,3,5]
  v = clfft.version()
  v ≥ v"2.8.0"  && push!(radices, 7)
  v ≥ v"2.12.0" && push!(radices, 11, 13)

  result = Int[]
  for s in size(A)
    fs = keys(factor(s))
    if fs ⊆ radices
      push!(result, 0)
    else
      # Try a closer number that has prime factors of 2 and 3.
      # Use the next power of 2 (say N) to get multiple new
      # candidates.
      N = nextpow2(s)

      # fractions of N: 100%, 93%, 84%, 75%, 56%
      candidates = [N, 15(N÷16), 27(N÷32), 3(N÷4), 9(N÷16)]
      candidates = candidates[candidates .> s]
      n = minimum(candidates)

      push!(result, n - s)
    end
  end

  result
end


text:
using Polyester

function flatten(dataset::AbstractArray, batch_dim::Int64)
    dims = ndims(dataset)
    len = 1
    for i in 1:dims
        if i!=batch_dim
            len *= size(dataset, i)
        end
    end
    input_data = zeros(Float32, len, size(dataset, batch_dim))
    @batch for i in axes(dataset, batch_dim)
        input_data[:,i] = Array{Float32}(reshape(selectdim(dataset, batch_dim, i), (len,1)))
    end
    return input_data
end


text:
@testset "write fixed width" begin
    res1 = """# Reduced           Mapped To
# === Header Start ===
    # of Elements
30

# === Header End ===
-4                            Int64[]
-3                            [3]
-2                            [2]
-1                            [1]
0                             Int64[]
"""
    io = IOBuffer()
    rm = find_classes(Mapping(x -> [-x]), -4:3)
    em = ExpandMap(rm, collect(-4:3))
    write_fixed_width(io, em, sorted=true)
    out = String(take!(io))
    @test out == res1
end

@testset "JLD" begin
    io = IOBuffer()
    rem = find_classes(Mapping(x -> [-x]), -4:3)
    em = ExpandMap(rem, collect(-4:3))
    write_JLD("tmp_save.jld", rem, em)
    d = load("tmp_save.jld")
    @test em == d["ExpandMap"]
    rm("tmp_save.jld")
end


text:
module Frtb

import ..Clarus

function ima(;params...)
  return Clarus.api_request("FRTB","IMA";params...)
end

function imaimpact(;params...)
  return Clarus.api_request("FRTB","IMAImpact";params...)
end

function modellablerf(;params...)
  return Clarus.api_request("FRTB","ModellableRF";params...)
end

function modellablerftrades(;params...)
  return Clarus.api_request("FRTB","ModellableRFTrades";params...)
end

function plvectors(;params...)
  return Clarus.api_request("FRTB","PLVectors";params...)
end

function sa(;params...)
  return Clarus.api_request("FRTB","SA";params...)
end

function saimpact(;params...)
  return Clarus.api_request("FRTB","SAImpact";params...)
end

end #MODULE-END


text:
@testset "preprocess" begin

@testset "kekulize" begin
    null = parse(SMILES, "")
    kekulize!(null)
    @test nodecount(null) == 0

    furan = parse(SMILES, "o1cccc1")
    kekulize!(furan)
    @test edgeattr(furan, 2).order == 2
    @test edgeattr(furan, 4).order == 2

    pyrrole = parse(SMILES, "[nH]1cccc1")
    kekulize!(pyrrole)
    @test edgeattr(pyrrole, 3).order == 2
    @test edgeattr(pyrrole, 5).order == 2
    wrongpyrrole = parse(SMILES, "n1cccc1")
    @test_throws ErrorException kekulize!(wrongpyrrole)

    pyrene = parse(SMILES, "c1cc2cccc3ccc4cccc1c4c32")
    pyrene = kekulize(pyrene)
    @test sum(bondorder(pyrene) .== 1) == 11
    @test sum(bondorder(pyrene) .== 2) == 8

    pyridineoxide = parse(SMILES, "[n+]1([O-])ccccc1")
    pyridineoxide = kekulize(pyridineoxide)
    @test sum(bondorder(pyridineoxide) .== 2) == 3

    c60 = parse(SMILES, "c12c3c4c5c1c6c7c8c2c9c1c3c2c3c4c4c%10c5c5c6c6c7c7c%11c8c9c8c9c1c2c1c2c3c4c3c4c%10c5c5c6c6c7c7c%11c8c8c9c1c1c2c3c2c4c5c6c3c7c8c1c23")
    c60 = kekulize(c60)
    @test sum(bondorder(c60) .== 2) == 30
end

@testset "allhydrogens" begin
    ethanol = parse(SMILES, "[H]C([H])([H])C([H])([H])O")
    @test issetequal(allhydrogens(ethanol), [1, 3, 4, 6, 7])
    removed = removehydrogens(ethanol)
    @test nodecount(removed) == 3
    @test edgecount(removed) == 2
end

@testset "trivialhydrogens" begin
    ethanold3 = parse(SMILES, "[2H]C([2H])([2H])C([H])([H])O")
    @test issetequal(trivialhydrogens(ethanold3), [6, 7])
    removed = removehydrogens(ethanold3, all=false)
    @test nodecount(removed) == 6
    @test edgecount(removed) == 5
    LAla = parse(SMILES, "[H]N([H])[C@][H](C)C(=O)O")
    @test issetequal(trivialhydrogens(LAla), [1, 3])
end

@testset "addhydrogens" begin
    neop = parse(SMILES, "CC(C)(C)CO")
    neop = addhydrogens(neop)
    @test nodecount(neop) == 18
    @test edgecount(neop) == 17
end

@testset "largestcomponent" begin
    thiols = parse(SMILES, "CCCCS.SCCCCC")
    @test issetequal(largestcomponentnodes(thiols), 6:11)
    lc = extractlargestcomponent(thiols)
    @test nodecount(lc) == 6
    @test edgecount(lc) == 5
end

@testset "protonateacids" begin
    AcOH = parse(SMILES, "CC(=O)[O-]")
    AcOH = protonateacids(AcOH)
    @test nodeattr(AcOH, 4).charge == 0

    thiol = parse(SMILES, "CC[S-]")
    thiol = protonateacids(thiol)
    @test nodeattr(thiol, 3).charge == 0

    noxide = parse(SMILES, "[n+]1([O-])ccccc1")
    noxide = protonateacids(noxide)
    @test nodeattr(noxide, 1).charge == 1
end

@testset "deprotonateoniums" begin
    pyrrolidinium = parse(SMILES, "C1[N+]CCC1")
    pyrrolidinium = deprotonateoniums(pyrrolidinium)
    @test nodeattr(pyrrolidinium, 2).charge == 0

    oxonium = parse(SMILES, "[H][O+]([H])[H]")
    oxonium = deprotonateoniums(oxonium)
    @test nodeattr(oxonium, 2).charge == 0

    ammonium = parse(SMILES, "C[N+](C)(C)C")
    ammonium = deprotonateoniums(ammonium)
    @test nodeattr(ammonium, 2).charge == 1
end

@testset "polarize" begin
    acetone = parse(SMILES, "C[C+]([O-])C")
    acetone = depolarize(acetone)
    @test nodeattr(acetone, 2).charge == 0
    @test nodeattr(acetone, 3).charge == 0
    @test edgeattr(acetone, 2).order == 2
    acetone = polarize(acetone)
    @test nodeattr(acetone, 2).charge == 0
    @test nodeattr(acetone, 3).charge == 0
    @test edgeattr(acetone, 2).order == 2

    phosphorylCl = parse(SMILES, "[O-][P+](Cl)(Cl)Cl")
    phosphorylCl = depolarize(phosphorylCl)
    @test nodeattr(phosphorylCl, 1).charge == 0
    @test nodeattr(phosphorylCl, 1).charge == 0
    @test edgeattr(phosphorylCl, 1).order == 2
    phosphorylCl = polarize(phosphorylCl)
    @test nodeattr(phosphorylCl, 1).charge == 0
    @test nodeattr(phosphorylCl, 1).charge == 0
    @test edgeattr(phosphorylCl, 1).order == 2

    noxide = parse(SMILES, "[n]1(=O)ccccc1")
    noxide = polarize(noxide)
    @test nodeattr(noxide, 1).charge == 1
    @test nodeattr(noxide, 2).charge == -1
    @test edgeattr(noxide, 1).order == 1
    noxide = depolarize(noxide)
    @test nodeattr(noxide, 1).charge == 1
    @test nodeattr(noxide, 2).charge == -1
    @test edgeattr(noxide, 1).order == 1

    dmso = parse(SMILES, "CS(=O)C")
    dmso = polarize(dmso)
    @test nodeattr(dmso, 2).charge == 1
    @test nodeattr(dmso, 3).charge == -1
    @test edgeattr(dmso, 2).order == 1
    dmso = depolarize(dmso)
    @test nodeattr(dmso, 2).charge == 1
    @test nodeattr(dmso, 3).charge == -1
    @test edgeattr(dmso, 2).order == 1
    dmso = depolarize(dmso, positive=[:S])
    @test nodeattr(dmso, 2).charge == 0
    @test nodeattr(dmso, 3).charge == 0
    @test edgeattr(dmso, 2).order == 2
end

@testset "13dipole" begin
    me_azide = parse(SMILES, "C[N-][N+]#N")
    me_azide = toallenelike(me_azide)
    @test nodeattr(me_azide, 2).charge == 0
    @test nodeattr(me_azide, 4).charge == -1
    @test edgeattr(me_azide, 2).order == 2
    @test edgeattr(me_azide, 3).order == 2
    me_azide = totriplebond(me_azide)
    @test nodeattr(me_azide, 2).charge == -1
    @test nodeattr(me_azide, 4).charge == 0
    @test edgeattr(me_azide, 2).order == 1
    @test edgeattr(me_azide, 3).order == 3
end

end # preprocess


text:
module PackageCompilerX

using Base: active_project
using Libdl: Libdl
using Pkg: Pkg
using UUIDs: UUID

export create_sysimage, create_app, audit_app, restore_default_sysimage

include("juliaconfig.jl")

const NATIVE_CPU_TARGET = "native"
const APP_CPU_TARGET = "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)"

current_process_sysimage_path() = unsafe_string(Base.JLOptions().image_file)

all_stdlibs() = readdir(Sys.STDLIB)

yesno(b::Bool) = b ? "yes" : "no"

function get_compiler()
    cc = get(ENV, "JULIA_CC", nothing)
    if cc !== nothing
        return `$cc`
    end
    if Sys.which("gcc") !== nothing
        return `gcc`
    elseif Sys.which("clang") !== nothing
        return `clang`
    end
    if Sys.iswindows()
        if Sys.which("x86_64-w64-mingw32-gcc") !== nothing
            return `x86_64-w64-mingw32-gcc`
        end
    end
    error("could not find a compiler, looked for `gcc` and `clang`")
end

function get_julia_cmd()
    julia_path = joinpath(Sys.BINDIR, Base.julia_exename())
    cmd = `$julia_path --color=yes --startup-file=no`
end

function rewrite_sysimg_jl_only_needed_stdlibs(stdlibs::Vector{String})
    sysimg_source_path = Base.find_source_file("sysimg.jl")
    sysimg_content = read(sysimg_source_path, String)
    # replaces the hardcoded list of stdlibs in sysimg.jl with
    # the stdlibs that is given as argument
    return replace(sysimg_content, r"stdlibs = \[(.*?)\]"s => string("stdlibs = [", join(":" .* string.(stdlibs), ",\n"), "]"))
end

function create_fresh_base_sysimage(stdlibs::Vector{String}; cpu_target::String)
    tmp = mktempdir(cleanup=false)
    sysimg_source_path = Base.find_source_file("sysimg.jl")
    base_dir = dirname(sysimg_source_path)
    tmp_corecompiler_ji = joinpath(tmp, "corecompiler.ji")
    tmp_sys_ji = joinpath(tmp, "sys.ji")
    compiler_source_path = joinpath(base_dir, "compiler", "compiler.jl")

    @info "PackageCompilerX: creating base system image (incremental=false)..."
    cd(base_dir) do
        # Create corecompiler.ji
        cmd = `$(get_julia_cmd()) --cpu-target $cpu_target --output-ji $tmp_corecompiler_ji
                                  -g0 -O0 $compiler_source_path`
        @debug "running $cmd"
        read(cmd)

        # Use that to create sys.ji
        new_sysimage_content = rewrite_sysimg_jl_only_needed_stdlibs(stdlibs)
        new_sysimage_source_path = joinpath(base_dir, "sysimage_packagecompiler_x.jl")
        write(new_sysimage_source_path, new_sysimage_content)
        try
            cmd = `$(get_julia_cmd()) --cpu-target $cpu_target
                                      --sysimage=$tmp_corecompiler_ji
                                      -g1 -O0 --output-ji=$tmp_sys_ji $new_sysimage_source_path`
            @debug "running $cmd"
            read(cmd)
        finally
            rm(new_sysimage_source_path; force=true)
        end
    end

    return tmp_sys_ji
end

function run_precompilation_script(project::String, precompile_file::Union{String, Nothing})
    tracefile = tempname()
    if precompile_file == nothing
        arg = `-e ''`
    else
        arg = `$precompile_file`
    end
    touch(tracefile)
    cmd = `$(get_julia_cmd()) --sysimage=$(current_process_sysimage_path()) --project=$project
            --compile=all --trace-compile=$tracefile $arg`
    @debug "run_precompilation_script: running $cmd"
    read(cmd)
    return tracefile
end

function create_sysimg_object_file(object_file::String, packages::Vector{Symbol};
                            project::String,
                            base_sysimage::String,
                            precompile_execution_file::Vector{String},
                            precompile_statements_file::Vector{String},
                            cpu_target::String,
                            compiled_modules::Bool)
    # include all packages into the sysimg
    julia_code = """
        if !isdefined(Base, :uv_eventloop)
            Base.reinit_stdio()
        end
        Base.__init__();
        """
    for package in packages
        julia_code *= "using $package\n"
    end

    # handle precompilation
    precompile_statements = ""
    @debug "running precompilation execution script..."
    tracefiles = String[]
    for file in (isempty(precompile_execution_file) ? (nothing,) : precompile_execution_file)
        tracefile = run_precompilation_script(project, file)
        precompile_statements *= "    append!(precompile_statements, readlines($(repr(tracefile))))\n"
    end
    for file in precompile_statements_file
        precompile_statements *=
            "    append!(precompile_statements, readlines($(repr(file))))\n"
    end

    precompile_code = """
        # This @eval prevents symbols from being put into Main
        @eval Module() begin
            PrecompileStagingArea = Module()
            for (_pkgid, _mod) in Base.loaded_modules
                if !(_pkgid.name in ("Main", "Core", "Base"))
                    eval(PrecompileStagingArea, :(const \$(Symbol(_mod)) = \$_mod))
                end
            end
            precompile_statements = String[]
            $precompile_statements
            for statement in sort(precompile_statements)
                # println(statement)
                try
                    Base.include_string(PrecompileStagingArea, statement)
                catch
                    # See julia issue #28808
                    @error "failed to execute \$statement"
                end
            end
        end # module
        """
    julia_code *= precompile_code

    # finally, make julia output the resulting object file
    @debug "creating object file at $object_file"
    @info "PackageCompilerX: creating system image object file, this might take a while..."

    cmd = `$(get_julia_cmd()) --compiled-modules=$(yesno(compiled_modules)) --cpu-target=$cpu_target
                              --sysimage=$base_sysimage --project=$project --output-o=$(object_file) -e $julia_code`
    @debug "running $cmd"
    run(cmd)
end

default_sysimg_path() = joinpath(julia_private_libdir(), "sys." * Libdl.dlext)
default_sysimg_name() = basename(default_sysimg_path())
backup_default_sysimg_path() = default_sysimg_path() * ".backup"
backup_default_sysimg_name() = basename(backup_default_sysimg_path())

# TODO: Also check UUIDs for stdlibs, not only names
function gather_stdlibs_project(project::String)
    project_toml_path = abspath(Pkg.Types.projectfile_path(project; strict=true))
    ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
    @assert ctx.env.manifest !== nothing
    stdlibs = all_stdlibs()
    stdlibs_project = String[]
    for (uuid, pkg) in ctx.env.manifest
        if pkg.name in stdlibs
            push!(stdlibs_project, pkg.name)
        end
    end
    return stdlibs_project
end

"""
    create_sysimage(packages::Union{Symbol, Vector{Symbol}}; kwargs...)

Create a system image that includes the package(s) in `packages`.  An attempt
to automatically find a compiler will be done but can also be given explicitly
by setting the envirnment variable `JULIA_CC` to a path to a compiler

### Keyword arguments:

- `sysimage_path::Union{String,Nothing}`: The path to where
   the resulting sysimage should be saved. If set to `nothing` the keyword argument
   `replace_defalt` needs to be set to `true`.

- `project::String`: The project that should be active when the sysmage is created,
   defaults to the current active project.

- `precompile_execution_file::Union{String, Vector{String}}`: A file or list of
   files that contain code which precompilation statements should be recorded from.

- `precompile_statements_file::Union{String, Vector{String}}`: A file or list of
   files that contains precompilation statements that should be included in the sysimage.

- `incremental::Bool`: If `true`, build the new sysimage on top of the sysimage
   of the current process otherwise build a new sysimage from scratch. Defaults to `true`.

- `filter_stdlibs::Bool`: If `true`, only include stdlibs that are in the project file.
   Defaults to `false`, only set to `true` if you know the potential pitfalls.

- `replace_default::Bool`: If `true`, replaces the default system image which is automatically
   used when Julia starts. To replace with the one Julia ships with, use [`restore_default_sysimage()`](@ref)
"""
function create_sysimage(packages::Union{Symbol, Vector{Symbol}};
                         sysimage_path::Union{String,Nothing}=nothing,
                         project::String=active_project(),
                         precompile_execution_file::Union{String, Vector{String}}=String[],
                         precompile_statements_file::Union{String, Vector{String}}=String[],
                         incremental::Bool=true,
                         filter_stdlibs=false,
                         replace_default::Bool=false,
                         cpu_target::String=NATIVE_CPU_TARGET,
                         base_sysimage::Union{Nothing, String}=nothing,
                         compiled_modules=true)
    if replace_default==true
        if sysimage_path !== nothing
            error("cannot specify `sysimage_path` when `replace_default` is `true`")
        end
    end
    if sysimage_path === nothing
        if replace_default == false
            error("`sysimage_path` cannot be `nothing` if `replace_default` is `false`")
        end
        # We will replace the default sysimage so just put it somewhere for now
        tmp = mktempdir()
        sysimage_path = joinpath(tmp, string("sys.", Libdl.dlext))
    end
    if filter_stdlibs && incremental
        error("must use `incremental=false` to use `filter_stdlibs=true`")
    end

    # Functions lower down handles `packages` and precompilation file as arrays so convert here
    packages = vcat(packages)
    precompile_execution_file  = vcat(precompile_execution_file)
    precompile_statements_file = vcat(precompile_statements_file)

    if !incremental
        if base_sysimage !== nothing
            error("cannot specify `base_sysimage`  when `incremental=false`")
        end
        if filter_stdlibs
            stdlibs = gather_stdlibs_project(project)
        else
            stdlibs= all_stdlibs()
        end
        base_sysimage = create_fresh_base_sysimage(stdlibs; cpu_target=cpu_target)
    else
        if base_sysimage == nothing
            base_sysimage = current_process_sysimage_path()
        end
    end

    object_file = tempname() * ".o"
    create_sysimg_object_file(object_file, packages;
                              project=project,
                              base_sysimage=base_sysimage,
                              precompile_execution_file=precompile_execution_file,
                              precompile_statements_file=precompile_statements_file,
                              cpu_target=cpu_target,
                              compiled_modules=compiled_modules)
    create_sysimg_from_object_file(object_file, sysimage_path)
    if replace_default
        if !isfile(backup_default_sysimg_path())
            @debug "making a backup of default sysimg"
            cp(default_sysimg_path(), backup_default_sysimg_path())
        end
        @info "PackageCompilerX: default sysimg replaced, restart Julia for the new sysimg to be in effect"
        mv(sysimage_path, default_sysimg_path(); force=true)
    end
    # TODO: Remove object file
end

function create_sysimg_from_object_file(input_object::String, sysimage_path::String)
    julia_libdir = dirname(Libdl.dlpath("libjulia"))

    # Prevent compiler from stripping all symbols from the shared lib.
    # TODO: On clang on windows this is called something else
    if Sys.isapple()
        o_file = `-Wl,-all_load $input_object`
    else
        o_file = `-Wl,--whole-archive $input_object -Wl,--no-whole-archive`
    end
    extra = Sys.iswindows() ? `-Wl,--export-all-symbols` : ``
    cmd = `$(get_compiler()) -shared -L$(julia_libdir) -o $sysimage_path $o_file -ljulia $extra`
    @debug "running $cmd"
    run(cmd)
    return nothing
end

"""
    restore_default_sysimage()

Restores the default system image to the one that Julia shipped with.
Useful after running [`create_sysimage`](@ref) with `replace_default=true`.
"""
function restore_default_sysimage()
    if !isfile(backup_default_sysimg_path())
        error("did not find a backup sysimg")
    end
    cp(backup_default_sysimg_path(), default_sysimg_path(); force=true)
    rm(backup_default_sysimg_path())
    @info "PackageCompilerX: default sysimg restored, restart Julia for the new sysimg to be in effect"
    return nothing
end

const REQUIRES = "Requires" => UUID("ae029012-a4dd-5104-9daa-d747884805df")

"""
    audit_app(app_dir::String)

Check for possible problems with regards to relocatability for
the project at `app_dir`.

!!! warning
    This cannot guarantee that the project is free of relocatability problems,
    it can only detect some known bad cases and warn about those.
"""
function audit_app(app_dir::String)
    project_toml_path = abspath(Pkg.Types.projectfile_path(app_dir; strict=true))
    ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
    return audit_app(ctx)
end
function audit_app(ctx::Pkg.Types.Context)
    # Check for Requires.jl usage
    if REQUIRES in ctx.env.project.deps
        @warn "Project has a dependency on Requires.jl, code in `@require` will not be run"
    end
    for (uuid, pkg) in ctx.env.manifest
        if REQUIRES in pkg.deps
            @warn "$(pkg.name) has a dependency on Requires.jl, code in `@require` will not be run"
        end
    end

    # Check for build script usage
    if isfile(joinpath(dirname(ctx.env.project_file), "deps", "build.jl"))
        @warn "Project has a build script, this might indicate that it is not relocatable"
    end
    pkgs = Pkg.Types.PackageSpec[]
    Pkg.Operations.load_all_deps!(ctx, pkgs)
    for pkg in pkgs
        pkg_source = Pkg.Operations.source_path(pkg)
        pkg_source === nothing && continue
        if isfile(joinpath(pkg_source, "deps", "build.jl"))
            @warn "Package $(pkg.name) has a build script, this might indicate that it is not relocatable"
        end
    end
    return
end

"""
    create_app(app_source::String, compiled_app::String)

Compile an app with the source in `app_source` to the folder `compiled_app`.
The folder `app_source` needs to contain a package where the package include a
function with the signature

```
Base.@ccallable julia_main()::Cint
    # Perhaps do something based on ARGS
    ...
end
```

The executable will be placed in a folder called `bin` in `compiled_app` and
when the executabl run the `julia_main` function is called.

An attempt to automatically find a compiler will be done but can also be given
explicitly by setting the envirnment variable `JULIA_CC` to a path to a
compiler.

### Keyword arguments:

- `precompile_execution_file::Union{String, Vector{String}}`: A file or list of
   files that contain code which precompilation statements should be recorded from.

- `precompile_statements_file::Union{String, Vector{String}}`: A file or list of
   files that contains precompilation statements that should be included in the sysimage
   for the app.

- `incremental::Bool`: If `true`, build the new sysimage on top of the sysimage
   of the current process otherwise build a new sysimage from scratch. Defaults to `false`.

- `filter_stdlibs::Bool`: If `true`, only include stdlibs that are in the project file.
   Defaults to `false`, only set to `true` if you know the potential pitfalls.

- `audit::Bool`: Warn about eventual relocatability problems with the app, defaults
   to `true`.

- `force::Bool`: Remove the folder `compiled_app` if it exists before creating the app.
"""
function create_app(package_dir::String,
                    app_dir::String;
                    precompile_execution_file::Union{String, Vector{String}}=String[],
                    precompile_statements_file::Union{String, Vector{String}}=String[],
                    incremental=false,
                    filter_stdlibs=false,
                    audit=true,
                    force=false)
    project_toml_path = abspath(Pkg.Types.projectfile_path(package_dir; strict=true))
    if project_toml_path === nothing
        error("no project found in $(repr(package_dir))")
    end
    manifest_toml_path = abspath(Pkg.Types.manifestfile_path(package_dir))
    if manifest_toml_path === nothing
        @warn "it is not recommended to create an app without a preexisting manifest"
    end
    project_toml = Pkg.TOML.parsefile(project_toml_path)
    project_path = abspath(package_dir)
    app_name = get(project_toml, "name") do
        error("expected package to have a `name`-entry")
    end
    sysimg_file = app_name * "." * Libdl.dlext

    ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
    @debug "instantiating project at \"$project_toml_path\""
    Pkg.instantiate(ctx)

    if isdir(app_dir)
        if !force
            error("directory $(repr(app_dir)) already exists, use `force=true` to overwrite (will completely",
                  " remove the directory")
        end
        rm(app_dir; force=true, recursive=true)
    end

    audit && audit_app(ctx)

    mkpath(app_dir)

    bundle_julia_libraries(app_dir)
    bundle_artifacts(ctx, app_dir)

    # TODO: Create in a temp dir and then move it into place?
    binpath = joinpath(app_dir, "bin")
    mkpath(binpath)
    cd(binpath) do
        if !incremental
            tmp = mktempdir()
            # Use workaround at https://github.com/JuliaLang/julia/issues/34064#issuecomment-563950633
            # by first creating a normal "empty" sysimage and then use that to finally create the one
            # with the @ccallable function
            tmp_base_sysimage = joinpath(tmp, "tmp_sys.so")
            create_sysimage(Symbol[]; sysimage_path=tmp_base_sysimage, project=project_path,
                            incremental=false, filter_stdlibs=filter_stdlibs,
                            cpu_target=APP_CPU_TARGET)

            create_sysimage(Symbol(app_name); sysimage_path=sysimg_file, project=project_path,
                            incremental=true,
                            precompile_execution_file=precompile_execution_file,
                            precompile_statements_file=precompile_statements_file,
                            cpu_target=APP_CPU_TARGET,
                            base_sysimage=tmp_base_sysimage,
                            compiled_modules=false #= workaround julia#34076=#)
        else
            create_sysimage(Symbol(app_name); sysimage_path=sysimg_file, project=project_path,
                                              incremental=incremental, filter_stdlibs=filter_stdlibs,
                                              precompile_execution_file=precompile_execution_file,
                                              precompile_statements_file=precompile_statements_file,
                                              cpu_target=APP_CPU_TARGET,
                                              compiled_modules=false #= workaround julia#34076=#)
        end
        create_executable_from_sysimg(; sysimage_path=sysimg_file, executable_path=app_name)
        if Sys.isapple()
            cmd = `install_name_tool -change $sysimg_file @rpath/$sysimg_file $app_name`
            @debug "running $cmd"
            run(cmd)
        end
    end
    return
end

# This requires that the sysimg have been built so that there is a ccallable `julia_main`
# in Main.
function create_executable_from_sysimg(;sysimage_path::String,
                                        executable_path::String)
    flags = join((cflags(), ldflags(), ldlibs()), " ")
    flags = Base.shell_split(flags)
    wrapper = joinpath(@__DIR__, "embedding_wrapper.c")
    if Sys.iswindows()
        rpath = ``
    elseif Sys.isapple()
        rpath = `-Wl,-rpath,'@executable_path' -Wl,-rpath,'@executable_path/../lib'`
    else
        rpath = `-Wl,-rpath,\$ORIGIN:\$ORIGIN/../lib`
    end
    cmd = `$(get_compiler()) -DJULIAC_PROGRAM_LIBNAME=$(repr(sysimage_path)) -o $(executable_path) $(wrapper) $(sysimage_path) -O2 $rpath $flags`
    @debug "running $cmd"
    run(cmd)
    return nothing
end

function bundle_julia_libraries(app_dir)
    app_libdir = joinpath(app_dir, Sys.isunix() ? "lib" : "bin")
    cp(julia_libdir(), app_libdir; force=true)
    # We do not want to bundle the sysimg (nor the backup):
    rm(joinpath(app_libdir, "julia", default_sysimg_name()); force=true)
    rm(joinpath(app_libdir, "julia", backup_default_sysimg_name()); force=true)
    return
end

function bundle_artifacts(ctx, app_dir)
    @debug "bundling artifacts..."

    pkgs = Pkg.Types.PackageSpec[]
    Pkg.Operations.load_all_deps!(ctx, pkgs)

    # Also want artifacts for the project itself
    if ctx.env.pkg !== nothing
        # This is kinda ugly...
        ctx.env.pkg.path = dirname(ctx.env.project_file)
        push!(pkgs, ctx.env.pkg)
    end

    # Collect all artifacts needed for the project
    artifact_paths = String[]
    for pkg in pkgs
        pkg_source_path = Pkg.Operations.source_path(pkg)
        pkg_source_path === nothing && continue
        # Check to see if this package has an (Julia)Artifacts.toml
        for f in Pkg.Artifacts.artifact_names
            artifacts_toml_path = joinpath(pkg_source_path, f)
            if isfile(artifacts_toml_path)
                @debug "bundling artifacts for $(pkg.name)"
                artifact_dict = Pkg.Artifacts.load_artifacts_toml(artifacts_toml_path)
                for name in keys(artifact_dict)
                    @debug "  \"$name\""
                    push!(artifact_paths, Pkg.Artifacts.ensure_artifact_installed(name, artifacts_toml_path))
                end
                break
            end
        end
    end

    # Copy the artifacts needed to the app directory
    artifact_app_path = joinpath(app_dir, "artifacts")
    if !isempty(artifact_paths)
        mkpath(artifact_app_path)
    end
    for artifact_path in artifact_paths
        artifact_name = basename(artifact_path)
        # force=true?
        cp(artifact_path, joinpath(artifact_app_path, artifact_name))
    end
    return
end

end # module


text:
export to_long_mi, logsumexp_cuda,
    pop_cuda!, push_cuda!, all_empty, length_cuda,
    generate_all, generate_data_all, kfold

using DataFrames
using CUDA: CUDA

###################
# Misc.
####################

function to_long_mi(m::Matrix{Float64}, min_int, max_int)::Matrix{Int64}
    δmi = maximum(m) - minimum(m)
    δint = max_int - min_int
    return @. round(Int64, m * δint / δmi + min_int)
end

# TODO: get rid of all copies
@inline function logsumexp_cuda(x,y)
    Δ = ifelse(x == y, zero(x), CUDA.abs(x - y))
    max(x, y) + CUDA.log1p(CUDA.exp(-Δ))
end

###################
# Rudimentary CUDA-compatible stack data structure
####################

# sadly making `i` varargs doesn't work; kernel won't compile

function pop_cuda!(stack, i)
    if @inbounds stack[i,1] == zero(eltype(stack))
        return zero(eltype(stack))
    else
        @inbounds stack[i,1] -= one(eltype(stack))
        @inbounds return stack[i,stack[i,1]+2]
    end
end

function pop_cuda!(stack, i, j)
    if @inbounds stack[i,j,1] == zero(eltype(stack))
        return zero(eltype(stack))
    else
        @inbounds stack[i,j,1] -= one(eltype(stack))
        @inbounds return stack[i,j,stack[i,j,1]+2]
    end
end

function push_cuda!(stack, v, i)
    @inbounds stack[i,1] += one(eltype(stack))
    @inbounds CUDA.@cuassert 1+stack[i,1] <= size(stack, ndims(stack)) "CUDA stack overflow"
    @inbounds stack[i,1+stack[i,1]] = v
    return nothing
end

function push_cuda!(stack, v, i, j)
    @inbounds stack[i, j,1] += one(eltype(stack))
    @inbounds CUDA.@cuassert 1+stack[i, j,1] <= size(stack, ndims(stack)) "CUDA stack overflow"
    @inbounds stack[i, j,1+stack[i, j,1]] = v
    return nothing
end

all_empty(stack::AbstractArray{T,2}) where T =
    all(x -> iszero(x), stack[:,1])

all_empty(stack::AbstractArray{T,3}) where T =
    all(x -> iszero(x), stack[:,:,1])

length_cuda(stack, i...) = stack[i...,1]

###################
# One-Hot Encoding
####################

"""
One-hot encode data (2-D Array) based on categories (1-D Array)
Each row of the return value is a concatenation of one-hot encoding of elements of the same row in data
Assumption: both input arrays have elements of same type
"""
function one_hot_encode(X::Array{T, 2}, categories::Array{T,1}) where {T<:Any}
    X_dash = zeros(Bool, size(X)[1], length(categories)*size(X)[2])
    for i = 1:size(X)[1], j = 1:size(X)[2]
            X_dash[i, (j-1)*length(categories) + findfirst(==(X[i,j]), categories)] = 1
    end
    X_dash
end

###################
# Testing Utils
####################

"""
Given some missing values generates all possible fillings
"""
function generate_all(row::Vector)
    miss_count = count(ismissing, row)
    lits = length(row)
    result = Bool.(zeros(1 << miss_count, lits))

    if miss_count == 0
        result[1, :] = copy(row)
    else
        for mask = 0: (1<<miss_count) - 1
            cur = missings(Bool, lits)
            cur .= row
            cur[ismissing.(row)] = transpose(parse.(Bool, split(bitstring(mask)[end-miss_count+1:end], "")))
            result[mask+1,:] = cur
        end
    end
    DataFrame(result, :auto)
end

"""
Generates all possible binary configurations of size N
"""
function generate_data_all(N::Int)
    data_all = transpose(parse.(Bool, split(bitstring(0)[end-N+1:end], "")));
    for mask = 1: (1<<N) - 1
        data_all = vcat(data_all,
            transpose(parse.(Bool, split(bitstring(mask)[end-N+1:end], "")))
        );
    end
    DataFrame(data_all, :auto)
end

#####################
# K-fold partitioning
#####################

"Returns a(n index) partitioning a la k-fold."
function kfold(n::Int, p::Int)::Vector{Tuple{UnitRange, Vector{Int}}}
    F = Vector{Tuple{UnitRange, Vector{Int}}}(undef, p)
    j = s = 1
    k = n÷p
    for i ∈ 1:n%p
        if s > 1
            I = collect(1:s-1)
            if s+k < n append!(I, s+k+1:n) end
        else I = collect(s+k+1:n) end
        F[j] = (s:s+k, I)
        s += k+1
        j += 1
    end
    k = n÷p-1
    for i ∈ 1:p-n%p
        if s > 1
            I = collect(1:s-1)
            if s+k < n append!(I, s+k+1:n) end
        else I = collect(s+k+1:n) end
        F[j] = (s:s+k, I)
        s += k+1
        j += 1
    end
    return F
end


text: @warn "TODO: Implement NONCVXU2"


text:
function linearizepermutation(p1::NTuple{N₁,Int}, p2::NTuple{N₂}, n₁::Int, n₂::Int) where {N₁,N₂}
    p1′ = ntuple(StaticLength(N₁)) do n
        p1[n] > n₁ ? n₂+2n₁+1-p1[n] : p1[n]
    end
    p2′ = ntuple(StaticLength(N₂)) do n
        p2[N₂+1-n] > n₁ ? n₂+2n₁+1-p2[N₂+1-n] : p2[N₂+1-n]
    end
    return (p1′..., p2′...)
end

function permutation2swaps(perm)
    p = collect(perm)
    @assert isperm(p)
    swaps = Vector{Int}()
    N = length(p)
    for k = 1:N-1
        append!(swaps, p[k]-1:-1:k)
        for l = k+1:N
            if p[l] < p[k]
                p[l] += 1
            end
        end
        p[k] = k
    end
    return swaps
end

function _kron(A, B)
    sA = size(A)
    sB = size(B)
    s = map(*, sA, sB)
    C = similar(A, promote_type(eltype(A),eltype(B)), s)
    for IA in eachindex(IndexCartesian(), A)
        for IB in eachindex(IndexCartesian(), B)
            I = CartesianIndex(IB.I .+ (IA.I .- 1) .* sB)
            C[I] = A[IA]*B[IB]
        end
    end
    return C
end


text:
r = RandomSequences(10, 10)
@test r isa RandomSequences


text:
import Base: push!, length, show, getindex, setindex!, eachindex, isassigned,
             isapprox, zero, one, lastindex

export Leaf, Tape, Node, Branch, ∇

""" Basic unit on the computational graph."""
abstract type Node{T} end

""" A topologically ordered collection of Nodes. """
struct Tape
    tape::Vector{Any}
    Tape() = new(Vector{Any}())
    Tape(N::Int) = new(Vector{Any}(undef, N))
end
function show(io::IO, t::Tape)
    n = length(t)
    print(io, "Tape with ", n, " element", n == 1 ? "" : "s", n > 0 ? ":" : "")
    for i in eachindex(t)
        print(io, "\n  [", i, "]: ")
        if isassigned(tape(t), i)
            show(io, t[i])
        else
            print(io, "#undef")
        end
    end
end
@inline getindex(t::Tape, n::Int) = unthunk(getindex(tape(t), n))
@inline getindex(t::Tape, node::Node) = getindex(t, pos(node))
@inline lastindex(t::Tape) = length(t)
@inline setindex!(t::Tape, x, n::Int) = (tape(t)[n] = x; t)
@inline eachindex(t::Tape) = eachindex(tape(t))
@inline length(t::Tape) = length(tape(t))
@inline push!(t::Tape, node::Node) = (push!(tape(t), node); t)
@inline isassigned(t::Tape, n::Int) = isassigned(tape(t), n)
@inline isassigned(t::Tape, node::Node) = isassigned(t, pos(node))

# Make `Tape`s broadcast as scalars without a warning on 0.7
Base.Broadcast.broadcastable(tape::Tape) = Ref(tape)

"""
An element at the 'bottom' of the computational graph.

Fields:
val - the value of the node.
tape - The Tape to which this Leaf is assigned.
pos - the location of this Leaf in the tape to which it is assigned.
"""
struct Leaf{T} <: Node{T}
    val::T
    tape::Tape
    pos::Int
end
function Leaf(tape::Tape, val)
    leaf = Leaf(val, tape, length(tape) + 1)
    push!(tape, leaf)
    return leaf
end
show(io::IO, tape::Leaf{T}) where T = print(io, "Leaf{$T} $(unbox(tape))")
show(io::IO, tape::Leaf{T}) where T<:AbstractArray = print(io, "Leaf{$T} $(size(unbox(tape)))")

"""
A `Branch` is a Node with parents (args).

Fields:
val::T - the value of this node produced in the forward pass.
f - the function used to generate this Node.
args - Values indicating which elements in the tape will require updating by this node.
tape - The Tape to which this Branch is assigned.
pos - the location of this Branch in the tape to which it is assigned.
pullback::B - if there is a custom primative rule (a `ChainRulesCore.rrule`) then this holds
    the pullback to propagate gradients back through the operation. If there is not a rule
    then this is set to `nothing`.
    It may also be set to `nothing` by legacy Nabla rules that have not moved to ChainRules.
"""
struct Branch{T, B} <: Node{T}
    val::T
    f
    args::Tuple
    kwargs::NamedTuple
    tape::Tape
    pos::Int
    pullback::B
end
function Branch(f, args::Tuple, tape::Tape; kwargs...)
    unboxed = unbox.(args)

    # We could check for an `rrule` here if we wanted but we don't,
    # because we should never reach this point if we have an rrule
    primal_val = f(unboxed...; kwargs...)
    pullback = nothing

    branch = Branch(primal_val, f, args, kwargs.data, tape, length(tape) + 1, pullback)
    push!(tape, branch)
    return branch
end
show(io::IO, branch::Branch{T}) where T =
    print(io, "Branch{$T} $(unbox(branch)) f=$(getfield(branch, :f))")
show(io::IO, branch::Branch{T}) where T<:AbstractArray =
    print(io, "Branch{$T} $(size(unbox(branch))) f=$(getfield(branch, :f))")

"""
    tape(x::Node)
    tape(x::Tape)

Retrieve the `Tape` in a `Node`, or the underyling vector in a `Tape`.
"""
tape(x::Union{Node,Tape}) = getfield(x, :tape)

"""
    pos(x::Node)
    pos(x)

Location of Node on tape. -1 if not a Node object.
"""
pos(x::Node) = getfield(x, :pos)
pos(x) = -1

"""
    unbox(x::Node)
    unbox(x)

Get `.val` if `x` is a Node, otherwise is equivalent to `identity`.
"""
unbox(x::Node) = getfield(x, :val)
unbox(x) = x

isapprox(n::Node, f) = unbox(n) ≈ f
isapprox(f, n::Node) = n ≈ f
isapprox(n::Node, f::Node) = unbox(n) ≈ unbox(f)

zero(n::Node) = zero(unbox(n))
one(n::Node) = one(unbox(n))

# Let the user get the `size` and `length` of `Node`s.
Base.size(x::Node, dims...) = size(unbox(x), dims...)
Base.length(x::Node) = length(unbox(x))

# Leafs do nothing, Branches compute their own sensitivities and update others.
@inline propagate(y::Leaf, rvs_tape::Tape) = nothing
function propagate(y::Branch, rvs_tape::Tape)
    ȳ = rvs_tape[y]  # the gradient we are going to propagate through the operation in y
    d_tape = Nabla.tape(rvs_tape)  # strips off the Tape abstration leaving a plain Vector
    f = getfield(y, :f)
    args = getfield(y, :args)
    kwargs = getfield(y, :kwargs)
    xs = map(unbox, args)
    xids = map(pos, args)
    p = preprocess(f, y, ȳ, args...)  # inlining CSE will avoid unboxing twice.
    for j in eachindex(xs)
        x, xid = xs[j], xids[j]
        if xid > 0
            d_tape[xid] = isassigned(d_tape, xid) ?
                ∇(d_tape[xid], f, Arg{j}, p, unbox(y), ȳ, xs...; kwargs...) :  # maybe-inplace version
                ∇(f, Arg{j}, p, unbox(y), ȳ, xs...; kwargs...)
        end
    end
    return nothing
end

function propagate(fwd_tape::Tape, rvs_tape::Tape)
    for n in eachindex(rvs_tape)
        δ = length(rvs_tape) - n + 1
        isassigned(tape(rvs_tape), δ) && propagate(fwd_tape[δ], rvs_tape)
    end
    return rvs_tape
end

""" Initialise a Tape appropriately for being used as a reverse-tape. """
function reverse_tape(y::Node, ȳ)
    tape = Tape(pos(y))
    tape[end] = ȳ
    return tape
end

""" Used to flag which argument is being specified in x̄. """
struct Arg{N} end

"""
    ∇(y::Node{<:∇Scalar})
    ∇(y::Node{T}, ȳ::T) where T

Return a `Tape` object which can be indexed using `Node`s, each element of which contains
the result of multiplying `ȳ` by the transpose of the Jacobian of the function specified by
the `Tape` object in `y`. If `y` is a scalar and `ȳ = 1` then this is equivalent to
computing the gradient of `y` w.r.t. each of the elements in the `Tape`.

    ∇(f::Function, ::Type{Arg{N}}, p, y, ȳ, x...)

To implement a new reverse-mode sensitivity for the `N^{th}` argument of function `f`. p
is the output of `preprocess`. `x1`, `x2`,... are the inputs to the function, `y` is its
output and `ȳ` the reverse-mode sensitivity of `y`.

    ∇(x̄, f::Function, ::Type{Arg{N}}, p, y, ȳ, x...)

This is the optional in-place version of `∇` that should, if implemented, mutate
x̄ to have the gradient added to it.
"""
∇(y::Node, ȳ) = propagate(tape(y), reverse_tape(y, ȳ))
@inline ∇(y::Node{<:∇Scalar}) = ∇(y, one(unbox(y)))

@inline function ∇(x̄, f, ::Type{Arg{N}}, args...; kwargs...) where N
    return ChainRulesCore.add!!(x̄, ∇(f, Arg{N}, args...; kwargs...))
end

"""
    ∇(f; get_output::Bool=false)

Returns a function which, when evaluated with arguments that are accepted by `f`, will
return the gradient w.r.t. each of the arguments. If `get_output` is `true`, the result
of calling `f` on the given arguments is also returned.
"""
function ∇(f; get_output::Bool=false)
    return function(args...; kwargs...)
        args_ = Leaf.(Tape(), args)
        y = f(args_...; kwargs...)
        if y isa Node
            ∇f = ∇(y)
            ∇args = map(args_, args) do arg_, arg
                isassigned(∇f, arg_) ? ∇f[arg_] : zero(arg)
            end
        else
            ∇args = zero.(args)
        end
        ∇args_public = map(unthunk, ∇args)
        return get_output ? (y, ∇args_public) : ∇args_public
    end
end

# """
#     ∇(f::Function)

# Returns a function which, when evaluated with arguments that are accepted by `f` (`x`),
# will return a Tuple, the first element of which is the output of the function `f` and then
# second element of which is (yet another) function `g`. `g` can either be evaluated with no
# arguments, in which case it will return the gradient of `f` evaluated at `x`.
# Alternatively, it can be evaluated with arguments of the same type and shape as the output
# of `f(x)`, in which case it is equivalent to multiplying them 'from the left' by the
# Jacobian ∂(f(x)) / ∂x.
# """
# function ∇(f::Function)
#     return function(args...)
#         args_ = Leaf.(Tape(), args)
#         y = f(args_...)
#         ∇fx = (ȳ)->∇

#     end
# end

# A collection of methods for initialising nested indexable containers to zero.
for (f_name, scalar_init, array_init) in
    zip((:zerod_container, :oned_container, :randned_container),
        (:zero, :one, nothing),
        (:zeros, :ones, nothing))
    if scalar_init !== nothing
        @eval @inline $f_name(x::Number) = $scalar_init(x)
    end
    if array_init !== nothing
        @eval @inline $f_name(x::AbstractArray{<:Real}) = $array_init(eltype(x), size(x))
    end
    eval(quote
        @inline $f_name(x::Tuple) = map($f_name, x)
        @inline function $f_name(x)
            y = Base.copy(x)
            for n in eachindex(y)
                @inbounds y[n] = $f_name(y[n])
            end
            return y
        end
        $f_name(x::Ref) = Ref($f_name(x[]))
    end)
end
@inline randned_container(x::Number) = randn(typeof(x))
@inline randned_container(x::AbstractArray{<:Real}) = randn(eltype(x), size(x)...)
for T in (:Diagonal, :UpperTriangular, :LowerTriangular)
    @eval @inline randned_container(x::$T{<:Real}) = $T(randn(eltype(x), size(x)...))
end

# Bare-bones FMAD implementation based on internals of ForwardDiff.
# Accepts a Tuple of args and returns a Tuple of gradients.
# Currently scales almost exactly linearly with the number of inputs.
# The coefficient of this scaling could be improved by fully utilizing ForwardDiff
# and computing from multiple seeds at the same time.
function dual_call_expr(f, x::Type{<:Tuple}, ::Type{Type{Val{n}}}) where n
    dual_call = Expr(:call, :f)
    for m in 1:Base.length(x.parameters)
        push!(dual_call.args, n == m ? :(ForwardDiff.Dual(x[$m], 1)) : :(x[$m]))
    end
    return :(first(ForwardDiff.partials($dual_call)))
end
@generated fmad(f, x, n) = dual_call_expr(f, x, n)

function fmad_expr(f, x::Type{<:Tuple})
    body = Expr(:tuple)
    for n in 1:Base.length(x.parameters)
        push!(body.args, dual_call_expr(f, x, Type{Val{n}}))
    end
    return body
end
@generated fmad(f, x) = fmad_expr(f, x)


text:
__precompile__()

module DiffEqDSL

# package code goes here

end # module


text:
function testaddr(expr, expectedwithlines)
    Stochy.resetaddressingcounter()
    actual = Stochy.at(expr) |> Stochy.striplineinfo
    expected = expectedwithlines |> Stochy.striplineinfo
    @test actual == expected
    # if actual == expected
    #     println("*")
    # else
    #     println("Failure")
    #     println("Expected:")
    #     println(expected)
    #     println("Got:")
    #     println(actual)
    #     println()
    # end
end

# Values.

testaddr(:(0), :(0))
testaddr(:(x), :(x))
testaddr(:(false), :(false))
testaddr(:(:x), :(:x))
testaddr(Expr(:quote, :x), Expr(:quote, :x))
testaddr(:("foo"), :("foo"))

# Blocks.

testaddr(:(begin; end), :(begin; end))
testaddr(:(begin; x; end), :(begin; x; end))
testaddr(:(begin; x; y; end), :(begin; x; y; end))
testaddr(:(begin; x; begin; y; end; end), :(begin; x; begin; y; end; end))

# Function definitions.
# TODO: Lose these extra blocks.

testaddr(:(function f(); 0; end), :(function f($(Stochy.addressarg)); begin; 0; end; end))
testaddr(:(function f(x); 0; end), :(function f($(Stochy.addressarg), x); begin; 0; end; end))
testaddr(:(function f(x, y); 0; end), :(function f($(Stochy.addressarg), x, y); begin; 0; end; end))

testaddr(:(() -> false), :(($(Stochy.addressarg),) -> begin; false; end))
testaddr(:(x -> false), :(($(Stochy.addressarg), x) -> begin; false; end))
testaddr(:((x) -> false), :(($(Stochy.addressarg), x) -> begin; false; end))
testaddr(:((x,y) -> false), :(($(Stochy.addressarg), x, y) -> begin; false; end))

testaddr(:((args...) -> false), :(($(Stochy.addressarg), args...) -> begin; false; end))
testaddr(:((x, args...) -> false), :(($(Stochy.addressarg), x, args...) -> begin; false; end))

# Function application.
testaddr(:(f()), :(f(cons(1, $(Stochy.addressarg)))))
testaddr(:(f(x)), :(f(cons(1, $(Stochy.addressarg)), x)))

# Compound examples.

testaddr(:(function f(); g(); end), :(function f($(Stochy.addressarg)); begin; g(cons(1, $(Stochy.addressarg))); end; end))
testaddr(:(() -> f()), :(($(Stochy.addressarg),) -> begin; f(cons(1, $(Stochy.addressarg))); end))
testaddr(:(f(g())), :(f(cons(2, $(Stochy.addressarg)), g(cons(1, $(Stochy.addressarg))))))

# Other syntax.
testaddr(:(begin; local x = f(); end), :(begin; local x = f(cons(1, $(Stochy.addressarg))); end))
testaddr(:(f() ? g(0) : g(1)), :(f(cons(1, $(Stochy.addressarg))) ? g(cons(2, $(Stochy.addressarg)), 0) : g(cons(3, $(Stochy.addressarg)), 1)))
testaddr(:(f() && g()), :(f(cons(1, $(Stochy.addressarg))) && g(cons(2, $(Stochy.addressarg)))))
testaddr(:(f() || g()), :(f(cons(1, $(Stochy.addressarg))) || g(cons(2, $(Stochy.addressarg)))))
testaddr(:(f().g()), :(f(cons(1, $(Stochy.addressarg))).g(cons(2, $(Stochy.addressarg)))))
testaddr(:(f() == g()), :(f(cons(1, $(Stochy.addressarg))) == g(cons(2, $(Stochy.addressarg)))))
testaddr(:([f(), g()]), :([f(cons(1, $(Stochy.addressarg))), g(cons(2, $(Stochy.addressarg)))]))
testaddr(:(f()[g()]), :(f(cons(1, $(Stochy.addressarg)))[g(cons(2, $(Stochy.addressarg)))]))
testaddr(:([f()...]), :([f(cons(1,$(Stochy.addressarg)))...]))


text:
sig = [0.5 0; 0 0.5]
mvnormal = MvNormal([3, 5], sig)
mvnormal2 = MvNormal([2, 3], sig)
likelihood(x) = pdf(mvnormal, [x[1], x[2]]) + pdf(mvnormal2, [x[1], x[2]])

prior(x) = 1
priorsample(n) = rand(Uniform(0.001,10),(n, 2))

propsig = [0.2 0; 0 0.2]
propdist = MvNormal([0,0], propsig)
propdistsample() = rand(propdist)
propdistpdf(x) = pdf(propdist, x)

function term(sset, j)
    rv = true
    if j >= 1 rv = false end
    return rv
 end

@testset "Metropolis-Hasting" begin
    Random.seed!(8128)
    mhsample = mh(likelihood, prior, propdistpdf, propdistsample, [3, 1], 15000, 500)
    Random.seed!()

    @test size(mhsample, 1) == 14500
    @test isapprox(mean(mhsample, dims=1), [2.5 4], atol=0.05)
end

@testset "TMCMC" begin
    Random.seed!(8128)
    tmcmcsample = tmcmc(likelihood, prior, priorsample, 15000, 0.2)
    Random.seed!()

    @test size(tmcmcsample, 1) == 15000
    @test isapprox(mean(tmcmcsample, dims=1), [2.5 4], atol=0.05)
end

@testset "SMC" begin
    Random.seed!(8128)
    sample = smc(likelihood, prior, priorsample, propdistpdf, propdistsample, 15000, term)
    Random.seed!()

    @test size(sample, 1) == 15000
    @test isapprox(mean(sample, dims=1), [2.5 4], atol=0.05)
end

@testset "Gelman-Rubin_convergence" begin
    Random.seed!(8128)
    rhat = grconvergence(likelihood, prior, priorsample, propdistpdf, propdistsample, 15000, 8)
    Random.seed!()

    @test isapprox(rhat, 1, atol=0.01)
end


text:
abstract type AbstractBaseCampComputation end

mutable struct LogReturnComputationModel <: AbstractBaseCampComputation

    # data -
    ticker::String
    data::DataFrame
    map::Pair{Symbol,Symbol}
    from::Union{Date,Nothing}
    to::Union{Date,Nothing}

    # constructor -
    LogReturnComputationModel() = new()
end

mutable struct LinearReturnComputationModel <: AbstractBaseCampComputation

    # data -
    ticker::String
    data::DataFrame
    map::Pair{Symbol,Symbol}
    from::Union{Date,Nothing}
    to::Union{Date,Nothing}

    # constructor -
    LinearReturnComputationModel() = new()
end


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

N = 100
b_AU = 0.5
b_AC = 3
b_UX = -1
b_UB = 2
b_CB = -1.5
b_CY = 1
b_XY = 5

df = DataFrame(
  :a => rand(Normal(), N)
);
df[!, :u] = rand(Normal(1, 2), N) + b_AU * df[:, :a]
df[!, :c] = rand(Normal(0, 1), N) + b_AC * df[:, :a]
df[!, :b] = rand(Normal(1, 1), N) + b_UB * df[:, :u]
df[!, :x] = rand(Normal(-2, 1), N) + b_UX * df[:, :u]
df[!, :y] = rand(Normal(1, 2), N) + b_XY * df[:, :x] + b_CY * df[:, :c]

u = [:u]                    # Latent variable
d = OrderedDict(
  :y => [:c, :x],
  :b => [:c, :u],
  :x => :u,
  :c => :a,
  :u => :a
  )

dag = DAG("sr6_4_2", d, df=df);

fn = joinpath(mktempdir(), "sr6_4_2.dot")
to_graphviz(dag, fn)
#Sys.isapple() && run(`open -a GraphViz.app $(fn)`)

basisset = basis_set(dag)
allpaths  = all_paths(dag, :x, :y)
backdoorpaths = backdoor_paths(dag, allpaths, :x)
adjustmentsets = adjustment_sets(dag, :x, :y)

@testset "sr6_4_2" begin

  @test length(basisset) == 15
  @test basisset[6][1] == :x
  @test basisset[6][2] == :b
  @test basisset[6][3:end] == [:u]
  @test basisset[8][1] == :x
  @test basisset[8][2] == :a
  @test basisset[8][3:end] == [:u]

  @test allpaths[1] == [:x, :u, :b, :c, :y]
  @test allpaths[2] == [:x, :u, :a, :c, :y]

  @test backdoorpaths[1] == [:x, :u, :b, :c, :y]
  @test backdoorpaths[2] == [:x, :u, :a, :c, :y]

  @test adjustmentsets == [[:u], [:c], [:a]]

end


text:
#=

A call to app/render using tiny image.

>julia
julia> include("render.jl")
julia> @time test(image)
=#

# For "load" local  files
using FileIO
# The tested item
include("../../src/apps/render.jl")

# mock test data
# global so it persists and out of the test
image = load("/home/bootch/git/resynthesizerJulia/test/data/in/mediumlighthouse.png")

# selection mask for corpus
# medium
mask = falses(size(image))
mask[15:30,85:100] .= true

function test(image, mask)
    render(image, mask)
end

# Expect an image with texture of the lighthouse, everywhere

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using ITensors,
      Test

import ITensors: In, Out, Neither

@testset "QN Index" begin

  @testset "hasqns function" begin
    i = Index(4,"i")
    @test hasqns(i) == false
    j = Index(QN(0)=>1,QN(1)=>1)
    @test hasqns(j) == true
  end

  @testset "Array of QN Constructor" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    @test hasqns(i)
    @test dim(i) == 3
    @test hastags(i,"i")
  end

  @testset "Vararg Constructor" begin
    i = Index(QN(0)=>1,QN(1)=>2;tags="i")
    @test hasqns(i)
    @test dim(i) == 3
    @test hastags(i,"i")
    @test dir(i) == Out

    j = Index(QN(0)=>1,QN(1)=>2;tags="j",dir=In)
    @test hasqns(j)
    @test dim(j) == 3
    @test hastags(j,"j")
    @test dir(j) == In
  end

  @testset "flux and qn" begin
    i = dag(Index([QN(0)=>2, QN(1)=>2], "i"))

    @test flux(i => 1) == QN(0)
    @test flux(i => 2) == QN(0)
    @test flux(i => 3) == QN(-1)
    @test flux(i => 4) == QN(-1)
    @test flux(i => Block(1)) == QN(0)
    @test flux(i => Block(2)) == QN(-1)

    @test qn(i => 1) == QN(0)
    @test qn(i => 2) == QN(0)
    @test qn(i => 3) == QN(1)
    @test qn(i => 4) == QN(1)
    @test qn(i => Block(1)) == QN(0)
    @test qn(i => Block(2)) == QN(1)

  end
end

nothing


text:
using Quats
using Test

# write your own tests here
include("quaternion-test.jl")


text:
@testset "MPO" begin

    @testset "xxz" begin
        smpo = xxz_u1mpo(Float64, 4, 2, 1.0)
        mpo = xxz_mpo(Float64, 4, 2, 1.0);

        h2sym = contract(smpo.Ws[1], (1,2,-1,5),
                         smpo.Ws[4], (-1,3,4,6))
        h2sym = fuselegs(fuselegs(h2sym, 1, 3), 2, 3, true)
        # h2sym = dropdummyleg(h2sym, 3)
        # h2sym = dropdummyleg(h2sym, 1)

        h4sym = contract(
            contract(
                contract(smpo.Ws[1], (1,2,-1,5),
                         smpo.Ws[2], (-1,3,4,6)), (1,2,3,-1,6,7),
                smpo.Ws[3], (-1,4,5,8)), (1,2,3,4,-1,7,8,9),
            smpo.Ws[4], (-1,5,6,10))

        h4sym = fuselegs(fuselegs(h4sym, 1, 5), 2, 5, true)
        # h4sym = dropdummyleg(h4sym, 3)
        # h4sym = dropdummyleg(h4sym, 1)

        @tensor h2[:] :=
            mpo.Ws[1][-1,-2,1,-5] * mpo.Ws[4][1,-3,-4,-6]
        h2 = reshape(h2,4,4)

        @tensor h3[:] :=
            mpo.Ws[1][-1,-2,1,-6] * mpo.Ws[2][1,-3,2,-7] *
            mpo.Ws[4][2,-4,-5,-8]
        h3 = reshape(h3,8,8)

        @tensor h4[:] :=
            mpo.Ws[1][-1,-2,1,-7] * mpo.Ws[2][1,-3,2,-8] *
            mpo.Ws[3][2,-4,3,-9] * mpo.Ws[4][3,-5,-6,-10]
        h4 = reshape(h4,16,16)

        H2 = Matrix(xxz_hamiltonian(2, 1.0))
        H3 = Matrix(xxz_hamiltonian(3, 1.0))
        H4 = Matrix(xxz_hamiltonian(4, 1.0))

        @test h2 ≈ H2
        @test h3 ≈ H3
        @test h4 ≈ H4

        nums = collect(0:15)
        index4 = count_ones.(nums) .== 4
        index3 = count_ones.(nums) .== 3
        index2 = count_ones.(nums) .== 2
        index1 = count_ones.(nums) .== 1
        index0 = count_ones.(nums) .== 0

        @test h4[index0,index0] == h4sym.blocks[Sector(U1(0), U1(0))]
        @test h4[index1,index1] == h4sym.blocks[Sector(U1(1), U1(1))]
        @test h4[index2,index2] == h4sym.blocks[Sector(U1(2), U1(2))]
        @test h4[index3,index3] == h4sym.blocks[Sector(U1(3), U1(3))]
        @test h4[index4,index4] == h4sym.blocks[Sector(U1(4), U1(4))]

    end

    @testset "xxzlong" begin
        delta, r = 2.0, 0.6

        mpo = xxzlong_mpo(Float64,4,2,delta,r)
        @tensor h2[:] :=
            mpo.Ws[1][-1,-2,1,-5] * mpo.Ws[4][1,-3,-4,-6]
        h2 = reshape(h2,4,4)

        @tensor h3[:] :=
            mpo.Ws[1][-1,-2,1,-6] * mpo.Ws[2][1,-3,2,-7] *
            mpo.Ws[4][2,-4,-5,-8]
        h3 = reshape(h3,8,8)

        @tensor h4[:] :=
            mpo.Ws[1][-1,-2,1,-7] * mpo.Ws[2][1,-3,2,-8] *
            mpo.Ws[3][2,-4,3,-9] * mpo.Ws[4][3,-5,-6,-10]
        h4 = reshape(h4,16,16)

        H2 = xxz_longrange(2,delta,r)
        H3 = xxz_longrange(3,delta,r)
        H4 = xxz_longrange(4,delta,r)

        @test h2 ≈ Matrix(H2)
        @test h3 ≈ Matrix(H3)
        @test h4 ≈ Matrix(H4)
    end
end

@testset "MPOgen" begin

    @testset "j1j2" begin
        n = 4
        lx = 2*n
        j1, j2 = 1.0, 0.3
        model = j1j2model(lx, j1, j2)
        mpo = generatempo(model)
        legacympo = MatrixProductStateTools.j1j2_mpo(lx, j1, j2, 2)
        @test mpo2hamiltonian(legacympo) ≈ mpo2hamiltonian(mpo)

        tmodel = triangularspinmodel((2, n),
                                     j2, j1, j1, 0.0, 0.0, 0.0,
                                     boundary=(:OBC, :OBC))
        tmpo = generatempo(tmodel)
        @test mpo2hamiltonian(tmpo) ≈ mpo2hamiltonian(mpo)
    end

    @testset "j1j2sym" begin
        j1, j2 = 1.0, 0.3
        for lx in [3, 4]
            smodel = j1j2model(lx, j1, j2, symmetry=U1)
            smpo = generatempo(smodel)

            model = j1j2model(lx, j1, j2)
            mpo = generatempo(model)

            sH = mpo2hamiltonian(smpo)
            H = mpo2hamiltonian(mpo)

            nums = collect(0:2^lx-1)
            indexes = [count_ones.(nums) .== m for m in 0:lx]
            for n = 1:lx
                @test H[indexes[n], indexes[n]] ≈ sH.blocks[Sector(U1(n-1), -U1(n-1))]
            end
        end

        for n in [2, 3]
            smodel = j1j2model(2*n, j1, j2, symmetry=U1)
            smpo = generatempo(smodel)

            tmodel = triangularspinmodel((2, n),
                                         0.3, 1.0, 1.0, 0.0, 0.0, 0.0,
                                         boundary=(:OBC, :OBC),
                                         symmetry=U1)
            tmpo = generatempo(tmodel)
            @test mpo2hamiltonian(tmpo) ≈ mpo2hamiltonian(smpo)
        end
    end

    @testset "2legtriangle" begin
        ⊗ = kron
        n = 4
        j1, j2, k1, k2 = 1.0, 0.3, 0.4, 0.7
        k = 0.3
        H = j1j2_explicit(2*n, j1, j2)
        ringop = ringexchangeoperator(4)
        for i=1:2*n-3
            if i % 2 == 1
                H += k1 * I(2^(i-1)) ⊗ permutespins(ringop, [1,3,4,2]) ⊗ I(2^(2*n-i-3))
            else
                H += k2 * I(2^(i-1)) ⊗ permutespins(ringop, [1,3,4,2]) ⊗ I(2^(2*n-i-3))
            end
        end

        model = triangularspinmodel((2, n), j2, j1, j1, k1, k2, 0.,
                                    boundary=(:OBC, :OBC))
        mpo = generatempo(model)
        Hmpo = mpo2hamiltonian(mpo)

        nk1s = div(2*n-3, 2) + (2*n-3) % 2
        nk2s = div(2*n-3, 2)
        @test Hmpo + 0.25*I(2^(2*n))*k1*nk1s + 0.25*I(2^(2*n))*k2*nk2s ≈ H
    end

    @testset "2legtrianglesym" begin
        j1, j2, j3, k1, k2, k3 = rand(6)
        for n in [2, 3]
            smodel = triangularspinmodel((2,n), j1, j2, j3, k1, k2, k3, symmetry=U1)
            smpo = generatempo(smodel)

            model = triangularspinmodel((2,n), j1, j2, j3, k1, k2, k3)

            mpo = generatempo(model)

            sH = mpo2hamiltonian(smpo)
            H = mpo2hamiltonian(mpo)

            lx = 2*n
            nums = collect(0:2^lx-1)
            indexes = [count_ones.(nums) .== m for m in 0:lx]
            for n = 1:lx
                @test H[indexes[n], indexes[n]] ≈ sH.blocks[Sector(U1(n-1), -U1(n-1))]
            end
        end
    end
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
function save_output(env, solution)

    return
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
using Plots
using StatsBase
using Distributions
using StatsPlots

# %%
n = 1_000
p_grid = range(0, 1, length = n)
prior = ones(n)
# prior[p_grid .< 0.5] .= 0.0
likelihood = @. pdf(Binomial(15, p_grid), 8)
posterior = likelihood .* prior
posterior ./= sum(posterior)

plot(p_grid, posterior)

# %%
ns = 10_000
weights = pweights(posterior)
samples = sample(p_grid, weights, ns)

density(samples)
plot!(p_grid, posterior * n)

# %%
using Turing
chn = Chains(samples, ["waters"])
hpd(chn; alpha = 0.1)

mix = MixtureModel(Binomial.(15, samples))
predpos_samples = rand(mix, 10_000)
histogram(predpos_samples, normalize = :probability)
mean(predpos_samples) do s
    s == 8
end

mix = MixtureModel(Binomial.(9, samples))
predpos_samples = rand(mix, 10_000)
histogram(predpos_samples, normalize = :probability)
mean(predpos_samples) do s
    s == 6
end


text:
"""
    energyplot(data; kwargs...)
    energyplot!(data; kwargs...)

Plot energy transition distribution and marginal energy distribution in HMC algorithms.

This may help to diagnose poor exploration by gradient-based algorithms like HMC or NUTS.

# Positional Arguments
- `data`: any object that can be converted to an `ArviZ.Dataset` representing a
    `sample_stats` group with an `energy` variable.

# Keyword Arguments
- `showbfmi=true`: If `true` add to the plot the value of the estimated Bayesian fraction of
    missing information.
- `kwargs`: Additional attributes understood by Plots.jl. In particular, see keyword
    arguments for [`distplot`](@ref).
"""
RecipesBase.@userplot EnergyPlot

RecipesBase.@recipe function f(p::EnergyPlot; showbfmi=true)
    data = only(p.args)
    energy = ArviZ.convert_to_dataset(data; group=:sample_stats).energy.values

    if showbfmi
        bfmi = ArviZ.bfmi(energy)
        nchains = length(bfmi)
        labels = ["Chain $i BFMI = $(@sprintf("%.2f",bfmi[i]))" for i in 1:nchains]
        # hack to just add legend entries
        # TODO: work out how to do this without incrementing color cycle. Ideally would
        # use `primary=false`, but this removes the series from the legend
        RecipesBase.@series begin
            seriestype := :path
            label := permutedims(labels)
            seriesalpha := false
            seriescolor := :white
            [[0] for _ in 1:nchains], [[0] for _ in 1:nchains]
        end
    end

    energy_marginal = vec(energy .- StatsBase.mean(energy))
    energy_trans = vec(diff(energy; dims=2))
    seriestype := distplot
    label := ["Marginal energy" "Energy transition"]
    legend --> :best
    fill --> get(plotattributes, :kind, :default) !== :ecdf
    fillalpha --> [1 0.75]
    # use first 2 colors in color cycle
    seriescolor --> [1 2]
    grid --> false
    ticks --> nothing
    return [energy_marginal, energy_trans]
end


text:
using Test
using CyclicArrays

#=
face: 1:nfaces
direction:X=1 / Y=2 / Z=3,
side: - = 1 / + = 2
destiny:
1  face: 1:nfaces
2  direction:X=1 / Y=2 / Z=3,
3  side: - = 1 / + = 2
4  flip 1 / no flip 0
=#

# 3d array example

nx=4
ny=4
nz=1
nfaces=1;

faces=zeros(nfaces,3,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];
faces[1,2,1,:]=[1,2,2,0];
faces[1,2,2,:]=[1,2,1,0];
faces[1,3,1,:].=-1;
faces[1,3,2,:].=-1;

grid=CyclicArray(faces);

x0=reshape(collect(0:15).%4 .+1,(4,4))

x=reshape(x0,(1,4,4));
y=reshape(x0',(1,4,4));

x1=CyclicArray(x,grid)
y1=CyclicArray(y,grid)

i=5;j=1;k=1;
x1[k,j,i],y1[k,j,i]

x1[k,j,i]=2

@test x1[1,5,1]==2
@test [x1[1,2,2],y1[1,2,2]]==[2,2]
@test [x1[1,3,2],y1[1,3,2]]==[3,2]
@test [x1[1,5,2],y1[1,5,2]]==[1,2]
@test [x1[1,2,5],y1[1,2,5]]==[2,1]

x2=CyclicArray(x1.data,"3D")
@test all(Bool.((x1==x2).data))

println("Pass 3D test")

# 2d array example

nx=5
ny=5
nfaces=1;

faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];
faces[1,2,1,:]=[1,2,2,0];
faces[1,2,2,:]=[1,2,1,0];

grid=CyclicArray(faces);

x=reshape(collect(0:15).%4 .+1,(4,4));
y=x'

x1=CyclicArray(x,grid)
y1=CyclicArray(y,grid)

i=5;j=1;
x1[i,j],y1[i,j]

x1[i,j]=2

@test x1[5,1]==2
@test [x1[2,2],y1[2,2]]==[2,2]
@test [x1[3,2],y1[3,2]]==[3,2]
@test [x1[5,2],y1[5,2]]==[1,2]
@test [x1[2,5],y1[2,5]]==[2,1]
@test x1[1:2,1:2]==x[1:2,1:2]
@test x1[5:6,5:6]==x[1:2,1:2]

z1=CyclicArray(ones(1,4,4),grid)

@test z1.data==(z1*z1.data).data
@test z1.data==(z1+z1.data).data/2
@test z1.data==(z1-z1.data).data.+1
@test z1.data==(z1^z1.data).data
@test z1.data==(z1/z1.data).data
@test z1.data==(z1\z1.data).data

z1[1,1,1]=2

@test stagger(z1,dims=2,frac=0.7)[1,1,1]==1.3
@test shiftc(z1,dims=2,shift=1)[1,2,1]==2

x2=CyclicArray(x1.data,"2D")
@test all(Bool.((x1==x2).data))

println("Pass 2D test")

# 1d array example

faces=zeros(1,1,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];

grid=CyclicArray(faces);
x=[0:3;]
x1=CyclicArray(x,grid)

@test x1[5]==0
@test x1[0]==3

x2=CyclicArray(x,"1D")

@test all(Bool.((x1==x2).data))

println("Pass 1D test")

# 1d array example 2 faces

x=zeros(2,4)
x[1,:]=[0:3;]
x[2,:]=[4:7;]
x1=CyclicArray(x,"1D2F")

@test x1[1,0]==7
@test x1[1,5]==4
@test x1[1,9]==0
@test x1[2,5]==0
@test x1[2,0]==3

println("Pass 1D2F test")

# cubed sphere example

nfaces=6
faces=zeros(Int64,nfaces,3,2,4);

faces[1,1,1,:]=[4,1,2,0]
faces[1,1,2,:]=[2,1,1,0]
faces[1,2,1,:]=[6,2,2,0]
faces[1,2,2,:]=[5,2,1,0]
faces[1,3,1,:].=-1
faces[1,3,2,:].=-1

faces[2,1,1,:]=[1,1,2,0]
faces[2,1,2,:]=[3,1,1,0]
faces[2,2,1,:]=[6,1,2,1]
faces[2,2,2,:]=[5,1,2,0]
faces[2,3,1,:].=-1
faces[2,3,2,:].=-1

faces[3,1,1,:]=[2,1,2,0]
faces[3,1,2,:]=[4,1,1,0]
faces[3,2,1,:]=[6,2,1,1]
faces[3,2,2,:]=[5,2,2,1]
faces[3,3,1,:].=-1
faces[3,3,2,:].=-1

faces[4,1,1,:]=[3,1,2,0]
faces[4,1,2,:]=[1,1,1,0]
faces[4,2,1,:]=[6,1,1,0]
faces[4,2,2,:]=[5,1,1,1]
faces[4,3,1,:].=-1
faces[4,3,2,:].=-1

faces[5,1,1,:]=[4,2,2,1]
faces[5,1,2,:]=[2,2,2,0]
faces[5,2,1,:]=[1,2,2,0]
faces[5,2,2,:]=[3,2,2,0]
faces[5,3,1,:].=-1
faces[5,3,2,:].=-1

faces[6,1,1,:]=[4,2,1,0]
faces[6,1,2,:]=[2,2,1,1]
faces[6,2,1,:]=[3,2,1,1]
faces[6,2,2,:]=[1,2,1,0]
faces[6,3,1,:].=-1
faces[6,3,2,:].=-1

grid=CyclicArray(faces);

x0=repeat(reshape(collect(0:15).%4 .+1,(1,1,4,4)),6,1,1,1)
y0=repeat(permutedims(reshape(collect(0:15).%4 .+1,(1,1,4,4)),(1,2,4,3)),6,1,1,1)
a0=repeat(collect(1:6),1,1,4,4)

x1=CyclicArray(x0,grid)
y1=CyclicArray(y0,grid)
a1=CyclicArray(a0,grid)

f=1;i=5;j=1;k=1;
a1[f,k,j,i],x1[f,k,j,i],y1[f,k,j,i]
@test [a1[f,k,j,i],x1[f,k,j,i],y1[f,k,j,i]]==[2,1,1]
f=1;i=1;j=5;k=1;
@test [a1[f,k,j,i],x1[f,k,j,i],y1[f,k,j,i]]==[5,1,1]
f=1;i=0;j=1;k=1;
@test [a1[f,k,j,i],x1[f,k,j,i],y1[f,k,j,i]]==[4,1,4]
f=1;i=1;j=0;k=1;
@test [a1[f,k,j,i],x1[f,k,j,i],y1[f,k,j,i]]==[6,4,1]

x2=CyclicArray(x0,"cubed")
@test all(Bool.((x1==x2).data))

println("Pass cubed sphere 3D test")

# cubed sphere 2D example

x0=permutedims(reshape((1:1:96)',(4,4,6)),(3,2,1))
x1=CyclicArray(x0,"cubed2D")

@test x1[3,1,0]==x1[2,1,4]
@test x1[1,1,0]==x1[4,1,4]

println("Pass cubed sphere 2D test")

# 2d array example XY exchange

x=reshape(collect(0:15).+1,(4,4));
x2=CyclicArray(x,"2DXY")
@test x2[1,9]==1
@test x2[1,10]==5
@test x2[9,1]==1
@test x2[10,1]==2

println("Pass 2D exchange test")

# 2d array example XY exchange

x=reshape(collect(0:15).+1,(1,4,4));
x2=CyclicArray(x,"3DXY")
@test x2[1,1,9]==1
@test x2[1,1,10]==5
@test x2[1,9,1]==1
@test x2[1,10,1]==2

println("Pass 3D exchange test")

nfaces=1;
faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,2,2,1];
faces[1,1,2,:]=[1,2,1,1];
faces[1,2,1,:]=[1,1,2,0];
faces[1,2,2,:]=[1,1,1,0];

x=reshape(collect(0:15).+1,(4,4));
x2=CyclicArray(x,faces)

@test x2[1,5] == 13
@test x2[1,0] == 16
@test x2[0,1] == 13
@test x2[5,1] == 1

nfaces=1;
faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,1,2,1];
faces[1,1,2,:]=[1,1,1,1];
faces[1,2,1,:]=[1,2,2,1];
faces[1,2,2,:]=[1,2,1,1];

x=reshape(collect(0:15) .+1,(4,4));
x2=CyclicArray(x,faces)

@test x2[1,0]==x2[4,4]
@test x2[1,5]==x2[4,1]
@test x2[0,1]==x2[4,4]
@test x2[5,1]==x2[1,4]

x=reshape(collect(0:15) .+1,(4,4));
nfaces=1;
faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,1,1,0];
faces[1,1,2,:]=[1,1,2,0];
faces[1,2,1,:]=[1,2,1,0];
faces[1,2,2,:]=[1,2,2,0];

x2=CyclicArray(x,faces)
@test x2[1,0]==x2[1,1]
@test x2[1,5]==x2[1,4]
@test x2[0,1]==x2[1,1]
@test x2[5,1]==x2[4,1]

println("Pass 2D flip test")

# tests by sides

x=reshape(collect(0:15) .+1,(4,4));
nfaces=1;
faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,1,1,0];
faces[1,1,2,:]=[1,1,1,0];
faces[1,2,1,:]=[1,1,1,0];
faces[1,2,2,:]=[1,1,1,0];
faces[1,1,1,:]=[1,2,2,0];

faces[1,1,1,:]=[1,1,1,0];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[1,1]
@test x2[1,-1]==x2[1,2]

faces[1,1,1,:]=[1,2,1,0];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[1,1]
@test x2[1,-1]==x2[2,1]

faces[1,1,1,:]=[1,1,2,0];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[1,4]
@test x2[1,-1]==x2[1,3]

faces[1,1,1,:]=[1,2,2,0];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[4,1]
@test x2[1,-1]==x2[3,1]
#
faces[1,1,1,:]=[1,1,1,1];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[4,1]
@test x2[1,-1]==x2[4,2]

faces[1,1,1,:]=[1,2,1,1];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[1,4]
@test x2[1,-1]==x2[2,4]

faces[1,1,1,:]=[1,1,2,1];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[4,4]
@test x2[1,-1]==x2[4,3]

faces[1,1,1,:]=[1,2,2,1];
x2=CyclicArray(x,faces)
@test x2[1,0]==x2[4,4]
@test x2[1,-1]==x2[3,4]

println("Pass tests by sides 2D 1")

faces[1,1,2,:]=[1,1,1,0];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[1,1]
@test x2[1,6]==x2[1,2]

faces[1,1,2,:]=[1,2,1,0];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[1,1]
@test x2[1,6]==x2[2,1]

faces[1,1,2,:]=[1,1,2,0];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[1,4]
@test x2[1,6]==x2[1,3]

faces[1,1,2,:]=[1,2,2,0];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[4,1]
@test x2[1,6]==x2[3,1]
#
faces[1,1,2,:]=[1,1,1,1];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[4,1]
@test x2[1,6]==x2[4,2]

faces[1,1,2,:]=[1,2,1,1];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[1,4]
@test x2[1,6]==x2[2,4]

faces[1,1,2,:]=[1,1,2,1];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[4,4]
@test x2[1,6]==x2[4,3]

faces[1,1,2,:]=[1,2,2,1];
x2=CyclicArray(x,faces)
@test x2[1,5]==x2[4,4]
@test x2[1,6]==x2[3,4]

println("Pass tests by sides 2D 2")

faces[1,2,1,:]=[1,1,1,0];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[1,1]
@test x2[-1,1]==x2[1,2]

faces[1,2,1,:]=[1,2,1,0];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[1,1]
@test x2[-1,1]==x2[2,1]

faces[1,2,1,:]=[1,1,2,0];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[1,4]
@test x2[-1,1]==x2[1,3]

faces[1,2,1,:]=[1,2,2,0];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[4,1]
@test x2[-1,1]==x2[3,1]
#
faces[1,2,1,:]=[1,1,1,1];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[4,1]
@test x2[-1,1]==x2[4,2]

faces[1,2,1,:]=[1,2,1,1];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[1,4]
@test x2[-1,1]==x2[2,4]

faces[1,2,1,:]=[1,1,2,1];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[4,4]
@test x2[-1,1]==x2[4,3]

faces[1,2,1,:]=[1,2,2,1];
x2=CyclicArray(x,faces)
@test x2[0,1]==x2[4,4]
@test x2[-1,1]==x2[3,4]

println("Pass tests by sides 2D 3")

faces[1,2,2,:]=[1,2,1,0];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[1,1]
@test x2[6,1]==x2[2,1]

faces[1,2,2,:]=[1,1,2,0];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[1,4]
@test x2[6,1]==x2[1,3]

faces[1,2,2,:]=[1,2,2,0];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[4,1]
@test x2[6,1]==x2[3,1]
#
faces[1,2,2,:]=[1,1,1,1];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[4,1]
@test x2[6,1]==x2[4,2]

faces[1,2,2,:]=[1,2,1,1];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[1,4]
@test x2[6,1]==x2[2,4]

faces[1,2,2,:]=[1,1,2,1];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[4,4]
@test x2[6,1]==x2[4,3]

faces[1,2,2,:]=[1,2,2,1];
x2=CyclicArray(x,faces)
@test x2[5,1]==x2[4,4]
@test x2[6,1]==x2[3,4]

println("Pass tests by sides 2D 4")

faces=zeros(1,1,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];

x=[0:3;]
x1=CyclicArray(x,faces)

faces[1,1,1,:]=[1,1,1,0];
x1=CyclicArray(x,faces)
x1[0]==x1[1]
x1[-1]==x1[2]

faces[1,1,1,:]=[1,1,2,0];
x1=CyclicArray(x,faces)
x1[0]==x1[4]
x1[-1]==x1[3]

faces[1,1,2,:]=[1,1,1,0];
x1=CyclicArray(x,faces)
x1[5]==x1[1]
x1[6]==x1[2]

faces[1,1,2,:]=[1,1,2,0];
x1=CyclicArray(x,faces)
x1[5]==x1[4]
x1[6]==x1[3]

println("Pass tests by sides 1D")

max_issues_repo_name: udistr/CircularArrays


text:
using Base.Test
using DataFrames

let
	@assert isequal(DataVector[1, 2, NA],
		            DataArray(PooledDataVector[1, 2, NA]))

	# Test vector() and matrix() conversion tools
	dv = dataones(5)
	@assert isa(vector(dv), Vector{Float64})
	@assert isa(convert(Vector{Float64}, dv), Vector{Float64})
	dv[1] = NA
	# Should raise errors:
	# vector(dv)
	# convert(Vector{Float64}, dv)
	@assert isa(vector(dv, Any), Vector{Any})
	@assert isnan(vector(dv, Float64, NaN)[1])

	dm = dataones(3, 3)
	@assert isa(matrix(dm), Matrix{Float64})
	@assert isa(convert(Matrix{Float64}, dm), Matrix{Float64})
	dm[1, 1] = NA
	# Should raise errors:
	# matrix(dm)
	# convert(Matrix{Float64}, dm)
	@assert isa(matrix(dm, Any), Matrix{Any})
	@assert isnan(matrix(dm, Float64, NaN)[1, 1])

	df = DataFrame()
	df["A"] = 1:5
	df["B"] = ["A", "B", "C", "D", "E"]
	@assert isa(matrix(df), Matrix{Any})
	@assert isa(matrix(df, Any), Matrix{Any})

	df = DataFrame()
	df["A"] = 1:5
	df["B"] = 1.0:5.0
	@assert isa(matrix(df), Matrix{Real})
	@assert isa(matrix(df, Any), Matrix{Any})
	@assert isa(matrix(df, Float64), Matrix{Float64})

	df = DataFrame()
	df["A"] = 1.0:5.0
	df["B"] = 1.0:5.0
	@assert isa(matrix(df), Matrix{Float64})
	@assert isa(matrix(df, Any), Matrix{Any})
	@assert isa(matrix(df, Int), Matrix{Int})
end


text:
using StringDistances, Random
using BenchmarkTools

N = if length(ARGS) > 0
    try
        parse(Int, ARGS[1])
    catch _
        100
    end
else
    100 # default value
end

Maxlength = if length(ARGS) > 1
    try
        parse(Int, ARGS[2])
    catch _
        100
    end
else
    100 # default value
end

# If there are strings already cached to disk we start with them and only
# add new ones if needed.
using Serialization
const CacheFile = joinpath(@__DIR__(), "perfteststrings_$(Maxlength).juliabin")
SaveCache = false

S = if isfile(CacheFile)
    try
        res = deserialize(CacheFile)
        println("Read $(length(res)) strings from cache file: $CacheFile")
        res
    catch err
        String[]
    end
else
    println("Creating $N random strings.")
    SaveCache = true
    String[randstring(rand(3:Maxlength)) for _ in 1:N]
    SaveCache = true
    String[randstring(rand(3:Maxlength)) for _ in 1:N]
end

if length(S) < N
    for i in (length(S)+1):N
        push!(S, randstring(rand(3:Maxlength)))
    end
    SaveCache = true
end

if SaveCache
    println("Saving cache file with $(length(S)) strings: $CacheFile")
    serialize(CacheFile, S)
end

println("For ", Threads.nthreads(), " threads and ", N, " strings of max length ", Maxlength, ":")

dist = Cosine(2)
t1 = @belapsed dm1 = pairwise(dist, S; preprocess = false)
t2 = @belapsed dm2 = pairwise(dist, S; preprocess = true)

println("  - time WITHOUT pre-calculation: ", round(t1, digits = 3))
println("  - time WITH    pre-calculation: ", round(t2, digits = 3))
println("  - speedup with pre-calculation: ", round(t1/t2, digits = 1))


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04az(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 50
    LDB = NMAX
    LDC =   NMAX
    LDD = NMAX
    LDFG = NMAX
    LDQ = 2*NMAX
    LDU = NMAX
    LDWORK = 18*NMAX*NMAX + NMAX + max( 2*NMAX, 24 ) + 3
    LDZ = NMAX
    LIWORK = 2*NMAX + 9
    LZWORK = 8*NMAX + 28
    Z = Array{ComplexF64,2}(undef, LDZ,NMAX)
    B = Array{ComplexF64,2}(undef, LDB,NMAX)
    FG = Array{ComplexF64,2}(undef, LDFG,NMAX)
    D = Array{ComplexF64,2}(undef, LDD,NMAX)
    C = Array{ComplexF64,2}(undef, LDC,NMAX)
    Q = Array{ComplexF64,2}(undef, LDQ,2*NMAX)
    U = Array{ComplexF64,2}(undef, LDU,2*NMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX)
    ALPHAI = Array{Float64,1}(undef, NMAX)
    BETA = Array{Float64,1}(undef, NMAX)
    BWORK = Array{BlasBool,1}(undef, NMAX)
    ZWORK = Array{ComplexF64,1}(undef, LZWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    COMPQ =  vs[2][1]
    COMPU =  vs[3][1]
    N = parse(BlasInt, vs[4])
    if ( N<0 || N>NMAX || mod( N, 2 )!=0 )
        @error "Illegal N=$N"
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       Z[i,1:_jsz] .= parsex.(ComplexF64, vs[_i0+1:_i0+_jsz])
    end
    vs = String[]
    _isz,_jsz = (N÷2,N÷2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(ComplexF64, vs[_i0+1:_i0+_jsz])
    end
    vs = String[]
    _isz,_jsz = (N÷2,N÷2+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       FG[i,1:_jsz] .= parsex.(ComplexF64, vs[_i0+1:_i0+_jsz])
    end

# interp call 1
    close(f)

    INFO = SLICOT.mb04az!(JOB, COMPQ, COMPU, N, Z, B, FG, D, C, Q, U, ALPHAR, ALPHAI, BETA, LIWORK, BWORK)
    @test INFO == 0
    INFO == 0 || return

    M = N÷2
if ( LSAME( JOB, 'T' ) )
# interp output 1
    println(io, "Z:")
    _nc = N
    _nr = N
    show(io, "text/plain", Z[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "B:")
    _nc = N
    _nr = N
    show(io, "text/plain", UpperTriangular(B[1:_nr,1:_nc]))
    println(io)

# interp output 3
    println(io, "FG:")
    _nc = N
    _nr = N
    show(io, "text/plain", UpperTriangular(FG[1:_nr,1:_nc]))
    println(io)

# interp output 4
    println(io, "D:")
    _nc = N
    _nr = N
    show(io, "text/plain", D[1:_nr,1:_nc])
    println(io)

# interp output 5
    println(io, "C:")
    _nc = N
    _nr = N
    show(io, "text/plain", C[1:_nr,1:_nc])
    println(io)

end # if
if ( LSAME( COMPQ, 'C' ) )
# interp output 6
    println(io, "Q:")
    _nc = 2*N
    _nr = 2*N
    show(io, "text/plain", Q[1:_nr,1:_nc])
    println(io)

end # if
if ( LSAME( COMPU, 'C' ) )
# interp output 7
    println(io, "U:")
    _nc = 2*N
    _nr = N
    show(io, "text/plain", U[1:_nr,1:_nc])
    println(io)

end # if
# interp output 8
    println(io, "ALPHAR:")
    _nr = N
    show(io, "text/plain", ALPHAR[1:_nr])
    println(io)

# interp output 9
    println(io, "ALPHAI:")
    _nr = N
    show(io, "text/plain", ALPHAI[1:_nr])
    println(io)

# interp output 10
    println(io, "BETA:")
    _nr = N
    show(io, "text/plain", BETA[1:_nr])
    println(io)

end # run_mb04az()


text:
const EM = (mu  = 0.01215, Rp  = 6371E3, Rs  = 1.737E6, r12 = 384403e3, m1  = 5.972e24)
const SE = (mu  = 0.01215, Rp  = 6371E3, Rs  = 1.737E6, r12 = 384403e3, m1  = 5.972e24)
# When a system is added the function "createSystem(sys::System,G=6.673e-11::Float64)"
# has to be updated and a new struct created


text:
using Pkg; "C:/Users/Arman Angaji/OneDrive - Universität zu Köln/Dokumente/Uni-Köln/Masterarbeit/Workspace/Julia_Master/MasterProject_Julia" |> Pkg.activate

using StaticArrays, Distributions, DataFrames, StatsBase, Plots, ProgressMeter, Statistics, Random

##########################################
## "integrated wt size" distribution

function neutral_growth(endsize; b , d, μ, t= 0. )

    n = 1
    ndist = [1]
    times = [0.]

    while 0 < n < endsize

        p = rand()*(b+d)
        t += rand(Exponential(1.))/((b+d)*n)

        if p > d
            if rand() > μ
                n += 1
            end
        else
            n -= 1
        end
        push!(ndist, n)
        push!(times, t)
    end
    return n, ndist, times
end

let reps = 4000
    b, μ = 1., 0.2
    drange = 0.:0.025:0.9
    q = similar(collect(drange))
    for (i,d)=enumerate(drange)
        for _=1:reps
            q[i] += iszero(neutral_growth(500; b = b, d = d, μ = μ)[1])
        end
    end
    q ./= reps
    scatter(drange, q, legend=:none, xlab=:d, ylab=:q)
    plot!(drange, [min(1., d/(b*(1-μ))) for d in drange])
end

# haplotype extinction prob is q = d / (b*(1-μ))

let reps = 1000
    b, d, μ = 1., 0.3, 0.2
    ndist = Array{Array{Int,1}}(undef, reps)
    times = Array{Array{Float64,1}}(undef, reps)
    for i=1:reps
        out = neutral_growth(500; b = b, d = d, μ = μ)
        ndist[i] = out[2]
        times[i] = out[3]
    end

    # 1) wt/haplotype growth curve
    # timesmean = [mean(vcat(times...)[findall(isequal(n), vcat(ndist...))])  for n in 1:500]
    # plot(timesmean, log.(1:500), legend=:none)
    # plot!(timesmean, t -> (b*(1-μ)-d)*t, lw=3.)

    # 2) wt/haplotype sizes at eventtimes
    # histogram(vcat(ndist...))

    # 3) wt/haplotype sizes integrated over time
    dtimes = [ts[2:end].-ts[1:end-1] for ts in times]
    pop!.(ndist)
    hist = fit(Histogram, vcat(ndist...), weights(vcat(dtimes...)), nbins=500 )
    scatter(log.(hist.edges[1][2:end]), log.(hist.weights) .- log(sum(hist.weights)))
    plot!(log.(1:500), logn -> -logn - log(sum(1/n for n=1:500)))
    plot!(legend=:none, xlab=:log_n, ylab=:log_pn, title="wt size over time ~ 1/n" )

end

# 1) haplotypes grow at rate b*(1-μ) - d
# 2) the distribution of wt sizes at events (!) is uniform
# waiting times (between events) on the other hand are ∝ 1/n
# 3) so wt size integrated over time is ~1/n

##########################################
## haplotype-parent size at birth

function neutral_growth(endsize; b , d, μ, t= 0. )

    n = 1
    nmut = 0
    ndist = Int[]

    while 0 < n < endsize

        p = rand()*(b+d)
        t += rand(Exponential(1.))/((b+d)*n)
        if p > d
            if rand() < μ
                nmut += 1
                push!(ndist, n)
            else
                n += 1
            end
        else
            n -= 1
        end
    end
    return n, nmut, ndist
end

let reps = 10000
    b, d, μ = 1., 0.3, 0.2
    ndist = Array{Array{Int,1}}(undef, reps)
    nmut = Array{Int}(undef, reps)
    for i=1:reps
        out = neutral_growth(500; b = b, d = d, μ = μ)
        ndist[i] =  out[3]
        nmut[i] = out[2]
    end
    # histogram(vcat(ndist...))

    histogram(nmut, nbins=100, xlab=:m, lab="")
    vline!([μ*b/(b*(1-μ)-d)*500], lab="expected | surviving")
    vline!([mean(nmut)] , lab="mean")
    vline!([(1-d/(b*(1-μ)))*μ*b/(b*(1-μ)-d)*500  ], lab="expected")
end

# 1) following from wt size distribution: parent (haplotype) size at birth ~ Unif

# 2)
# haplotype growing to size M has μ*b/λ'*M mutants when conditioned on survival
# in general: p(m=0) = (1-q), p(m>0) ~ q*Gaussian(μ*b/λ'*M)
# here λ' is haplotype growth rate b*(1-μ) - d !!

##########################################
## haplotype turnover - estranged

let N = 200
    b , μ = 1., 0.3

    dup = 0.85  # b*(1-μ)
    drange = 0.:0.025:dup
    turnover = similar(collect(drange))

    for (i,d) = enumerate(drange)

        norph = 0
        ntot = 0
        for i=1:10000
            n, nmut, ndist = neutral_growth(N; b = b, d = d, μ = μ)
            norph += iszero(n) ? nmut : 0
            ntot += nmut
        end
        turnover[i] = norph/ntot
        println(d)
    end
    scatter(drange, turnover)

    estim = [ sum( min(1., d/(b*(1-μ)) )^n for n=1:N) / N for d = 0.:0.005:dup]
    plot!( 0.:0.005:dup , estim, legend=:none, xlab=:d, ylab=:turnover)
end

# savefig("singleclone_estranged_turnover_μ0.3.png")

##########################################
## wt background (clone-parent) size at birth

function neutral_growth(endsize; b , d, μ, t= 0. )

    n = 1
    nmut = 0
    ndist = Int[]

    while 0 < n < endsize

        p = rand()*(b+d)
        t += rand(Exponential(1.))/((b+d)*n)
        if p > d
            if rand() < μ
                nmut += 1
                push!(ndist, n)
            end
            n += 1
        else
            n -= 1
        end
    end
    return n, nmut, ndist
end

let reps = 10000
    b, d, μ = 1., 0.3, 0.2
    ndist = Array{Array{Int,1}}(undef, reps)
    nmut = Array{Int}(undef, reps)
    for i=1:reps
        out = neutral_growth(500; b = b, d = d, μ = μ)
        ndist[i] = out[3]
        nmut[i] = out[2]
    end
    histogram(vcat(ndist...), xlab=:N1, lab="")

    # histogram(nmut, nbins=100, xlab=:m, lab="")
    # vline!([μ*b/(b-d)*500], lab="expected | surviving")
    # vline!([mean(nmut)] , lab="mean")
    # vline!([(1-d/b)*μ*b/(b-d)*500  ], lab="expected")
end

# 1) wt (clone) size at birth is uniform N1 ~ unif

# 2)
# clone growing to size M has μ*β*M mutants when conditioned on survival
# in general: p(m=0) = (1-q), p(m>0) ~ q*Gaussian(μ*β*M)

include("neutral_haplotype_growth_v1.jl")
using TumorGrowth: clones_by_mutations, DataFrame

let M=20
    b, μ = 1., 0.1

    drange = collect(0.0:0.05:0.95)
    push!(drange, 0.975)
    reps = 300

    countorph = zeros(Float64, length(drange))

    for (i,dp)=enumerate(drange)
        q = dp/b

        k = 0
        ntot = 0
        for _=1:reps
            out = neutral_growth(10*M ; Nthresh=M, b = b, d = dp, μ = μ)
            tumor = DataFrame(out[1])
            obs = out[2][:, 2]
            isempty(obs) && continue
            mtypes = clones_by_mutations(tumor)[1]
            freqs = sum.(getproperty.(mtypes,:n))
            filter!(!iszero, freqs)
            isempty(freqs) && continue
            k += 1
            ntot += length(freqs)
            countorph[i] += count(iszero, 10*M .- freqs )
        end
        countorph[i] /= ntot

        print("d",dp," ")
        sleep(0.1)
    end

    estim = [ sum( min(1., d/b )^n for n=1:M) / M for d = 0.:0.005:1.]
    plot(0.:0.005:1., estim, lab="M$M", legend=:topleft, xlab=:d, ylab=:turnover)
    scatter!(drange, countorph, lab="M$M μ$μ")
end

############################################
# 1/n size distribution
# this simulation also shows the 1/n distribution of haplotype sizes over time

function subclonegrowth()
    Random.seed!(1234)
    tstep=0.001
    a=1.
    b=.5
    μ = 0.2
    nwtvec=zeros(Int,0)
    for run in 1:10000
        nwt=1
        while 0<nwt<1000
            if rand()<a*tstep*nwt
                if rand() > μ
                    nwt+=1
                end
            end
            if rand()<b*tstep*nwt
                nwt-=1
            end
            if rand()<0.001
                push!(nwtvec,nwt)
            end
        end
    end
    return(nwtvec)
end

ns = subclonegrowth()

hist = fit(Histogram, ns, nbins = maximum(ns))
scatter(log.(hist.edges[1]), log.(hist.weights./sum(hist.weights)) )
plot!(log.(midpoints(hist.edges[1])), logn -> -logn)

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
import MusicAlbums
include(joinpath(dirname(pathof(MusicAlbums)), "../test", "runtests.jl"))