text:
# Define labelled nodes from A to T
a = Node("A")
b = Node("B")
c = Node("C")
d = Node("D")
e = Node("E")
f = Node("F")
g = Node("G")
h = Node("H")
i = Node("I")
j = Node("J")
k = Node("K")
l = Node("L")
m = Node("M")
n = Node("N")
o = Node("O")
p = Node("P")
q = Node("Q")
r = Node("R")
s = Node("S")
t = Node("T")

# Link them
link!(a, b, 4.465957012470789)
link!(b, c, 32.620624765198855)
link!(c, d, 3.3650512429729265)
link!(d, e, 7.263546792740656)
link!(e, f, 10.193537222701895)
link!(f, g, 28.38523219145209)
link!(f, h, 4.956537575274135)
link!(e, i, 11.855818691462021)
link!(i, j, 1.9066195090479867)
link!(j, k, 10.627409065542043)
link!(i, l, 16.331078235033235)
link!(i, m, 14.711619982982215)
link!(c, n, 28.28461325395832)
link!(n, o, 7.41317192325869)
link!(o, p, 11.35796124012394)
link!(o, q, 0.6669348982306939)
link!(n, r, 7.594073912870435)
link!(r, s, 22.386191336343458)
link!(r, t, 0.6063155522048979)

#= The resulting tree should have the following Newick string from A:
    "((((((G,H)F,((K)J,L,M)I)E)D,((P,Q)O,(S,T)R)N)C)B)A;"
=#

tree = Tree(a)
show(a)
show(tree)

@testset "Node properties" begin
    @test label(a) == "A"
    label!(a, "Foo")
    @test label(a) == "Foo"
    label!(a, "A")
    @test getspecies(a) == getspecies(b) == 0
    setspecies!(a, 1)
    @test getspecies(a) == 1
    annotate!(a, "Bootstrap", 0.95)
    @test annotation(a, "Bootstrap") == 0.95
    deannotate!(a, "Bootstrap")
    @test_throws KeyError annotation(a, "Bootstrap")
end

@testset "Branch properties" begin
    show(getbranch(a, b))
    @test brlabel(a, b) == ""
    brlabel!(a, b, "Branchiphoo")
    @test brlabel(a, b) == "Branchiphoo"
    brannotate!(a, b, "Bootstrap", 0.95)
    @test brannotation(a, b, "Bootstrap") == 0.95
    brdeannotate!(a, b, "Bootstrap")
    @test_throws KeyError brannotation(a, b, "Bootstrap")
    @test brlength(a, b) == 4.465957012470789
    brlength!(a, b, nothing)
    @test isnothing(brlength(a, b))
    brlength!(a, b, 4.465957012470789)
    @test_throws Phylodendron2.InvalidTopology brlength(a, e)
end

@testset "Subtree functions" begin
    @test n_node(tree) == 20
    @test n_node(c, n) == 7
    @test n_node(c, b) == 2
    @test n_node(h, f) == 19
    @test n_tip(tree) == 10
    @test n_tip(c, n) == 4
    @test n_tip(c, b) == 1
    @test n_tip(h, f) == 9
    @test n_branch(tree) == 19
    @test n_branch(c, n) == 6
    @test n_branch(c, b) == 1
    @test n_branch(h, f) == 18
    @test nodelabels(tree) == string.('A':'T')
    @test nodelabels(n, o) == ["O", "P", "Q"]
    @test tiplabels(tree) == ["A", "G", "H", "K", "L", "M", "P", "Q", "S", "T"]
    @test find_nonsplitting(tree) == [b, d, j]
    @test tips(tree) == [a, g, h, k, l, m, p, q, s, t]
    dir = SpeciesDirectory(tiplabels(tree))
    tree.dir = dir
    setspecies!(tree)
    @test n_species(tree) == 10
    @test n_species(c, n) == 4
    @test n_species(c, b) == 1
end

@testset "Species directories" begin
    dir = SpeciesDirectory(tiplabels(tree))
    @test dir.list == tiplabels(tree)
    @test length(dir.dict) == length(dir.list)
    tree.dir = dir
    setspecies!(tree)
end

@testset "Tree properties" begin
    @test label(tree) == ""
    label!(tree, "Tree")
    @test label(tree) == "Tree"
    @test origin(tree) == a
    root!(tree, c)
    @test isrooted(tree)
    @test origin(tree) == c
    unroot!(tree)
    @test ! isrooted(tree)
    @test origin(tree) == c
    origin!(tree, a)
    @test origin(tree) == a
    annotate!(a, "Posterior probability", 0.95)
    @test annotation(a, "Posterior probability") == 0.95
    deannotate!(a, "Posterior probability")
    @test_throws KeyError annotation(a, "Posterior probability")
end

@testset "Tree cloning" begin
    tree2 = clone(tree)
	@test tree ≢ tree2
	@test newick_string(tree) == newick_string(tree2)
	@test tree.origin ≢ tree2.origin
	@test tree.rooted == tree2.rooted
end

@testset "Peek topology" begin
    @test n_neighbour(Node()) == 0
    @test neighbours(Node()) == Node[]
    @test n_neighbour(a) == 1
    @test neighbours(a) == [b]
    @test istip(a)
    @test areneighbours(a, b)
    @test n_neighbour(b) == 2
    @test neighbours(b) == [a, c]
    @test ! istip(b)
    @test areneighbours(b, a)
    @test areneighbours(b, c)
    @test ! areneighbours(b, d)
    @test n_neighbour(n) == 3
    @test neighbours(n) == [c, o, r]
    @test n_neighbour(i) == 4
    @test neighbours(i) == [e, j, l, m]
end

@testset "Traversals" begin
    @testset "Preorder" begin
        @test preorder_vector(tree) ==
            [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t]
        @test preorder(tree) ==
            [
                (a, a),(a, b),(b, c),(c, d),(d, e),(e, f),(f, g),(f, h),(e, i),(i, j),(j, k),(i, l),(i, m),(c, n),(n, o),(o, p),(o, q),(n, r),(r, s),(r, t)
            ]
        @test preorder_vector(e, d) == [d, c, b, a, n, o, p , q, r, s, t]
        @test preorder_vector(d, e) == [e, f, g, h, i, j, k, l, m]
        @test preorder_vector(f) ==
            [f, e, d, c, b, a, n, o, p, q, r, s, t, i, j, k, l, m, g, h]
        @test_throws Phylodendron2.InvalidTopology preorder(a, f)
    end
    @testset "Postorder" begin
        @test postorder_vector(tree) ==
            [t, s, r, q, p, o, n, m, l, k, j, i, h, g, f, e, d, c, b, a]
        @test postorder(tree) ==
            [
                (r, t), (r, s), (n, r), (o, q), (o, p), (n, o), (c, n), (i, m), (i, l), (j, k), (i, j), (e, i), (f, h), (f, g), (e, f), (d, e), (c, d), (b, c), (a, b), (a, a)
            ]
        @test postorder_vector(e, d) == [t, s, r, q, p, o, n, a, b, c, d]
        @test postorder_vector(d, e) == [m, l, k, j, i, h, g, f, e]
        @test_throws Phylodendron2.InvalidTopology postorder(a, f)
    end
end

@testset "Find nodes" begin
    origin!(tree, a)
    @test findfirst(x -> label(x) == "F", tree) == f
    @test findfirst(x -> label(x) == "F", tree.origin) == f
    @test isnothing(findfirst(x -> label(x) == "W", tree))
    @test findfirst(istip, d, e) == g
    @test findall(x -> label(x) == "F", tree) == [f]
    @test findall(x -> label(x) == "F", tree.origin) == [f]
    @test isempty(findall(x -> label(x) == "W", tree))
    @test findall(istip, d, e) == [g, h, k, l, m]
end

@testset "Grafting and plucking" begin
    x = Node()
    brlen0 = brlength(a, b)
    graft!(x, a, b, 0.4)
    @test neighbours(x) == [a, b]
    @test neighbours(a) == [x]
    @test neighbours(b) == [c, x]
    @test brlength(a, x) == 0.4 * brlen0
    @test brlength(x, b) == 0.6 * brlen0
    pluck!(x, a, b)
    @test neighbours(x) == Node[]
    @test neighbours(a) == [b]
    @test neighbours(b) == [c, a]
    @test brlength(a, b) == brlen0
    brax = Branch(nothing)
    brxb = Branch(0.1)
    graft!(x, a, b, brax, brxb)
    @test neighbours(x) == [a, b]
    @test neighbours(a) == [x]
    @test neighbours(b) == [c, x]
    @test isnothing(brlength(a, x))
    @test brlength(x, b) == 0.1
    pluck!(x, a, b)
    @test isnothing(brlength(a, b))
    brlength!(a, b, brlen0)
end

@testset "More topological manipulation" begin
    tree2 = deepcopy(tree)
    origin2 = origin(tree2)
    pluck_nonsplitting!(tree2)
    @test find_nonsplitting(tree2) == []
    @test origin(tree2) == origin2
    tree2 = deepcopy(tree)
    origin2 = origin(tree2)
    neworigin = find_nonsplitting(tree2)[1]
    lastorigin = neighbours(neworigin)[1]
    origin!(tree2, neworigin)
    pluck_nonsplitting!(tree2)
    @test find_nonsplitting(tree2) == []
    @test origin(tree2) == lastorigin
end

@testset "Node paths and distances" begin
    @test nodepath(s, f) == [s, r, n, c, d, e, f]
    @test nodepath(f, s) == [f, e, d, c, n, r, s]
    @test nodepath(f, f) == Node[]
    @test nodepath(f, e) == [f, e]
    @test nodepath(l, k) == [l, i, j, k]
    @test_throws Phylodendron2.InvalidTopology nodepath(a, Node())
    @test nodedist(s, f) == 5
    @test nodedist(f, s) == 5
    @test nodedist(f, f) == 0
    @test nodedist(f, e) == 0
    @test nodedist(l, k) == 2
    @test_throws Phylodendron2.InvalidTopology nodedist(a, Node())
    @test brdist(s, f) ==
        brlength(s, r) + brlength(r, n) + brlength(n, c) + brlength(c, d) + brlength(d, e) + brlength(e, f)
    @test brdist(s, f) == brdist(f, s)
    @test brdist(f, f) == 0.0
    @test brdist(f, e) == brlength(f, e)
    @test brdist(l, k) == brlength(l, i) + brlength(i, j) + brlength(j, k)
    @test_throws Phylodendron2.InvalidTopology brdist(a, Node())
end

@testset "Bipartitions" begin
    origin!(tree, a)
    @test ! isinformative(Bipartition(BitArray([0, 0, 0, 0])))
    @test isinformative(Bipartition(BitArray([0, 0, 0, 1])))
    @test istrivial(Bipartition(BitArray([0, 1, 1, 1])))
    @test istrivial(Bipartition(BitArray([0, 0, 0, 1])))
    @test ! istrivial(Bipartition(BitArray([0, 0, 0, 0])))
    @test are_compatible(
        Bipartition(BitArray([0, 1, 1, 1])),
        Bipartition(BitArray([0, 1, 1, 0]))
        )
    @test are_conflicting(
        Bipartition(BitArray([0, 1, 1, 0])),
        Bipartition(BitArray([0, 0, 1, 1]))
        )
    dir = SpeciesDirectory(nodelabels(tree))
    tree.dir = dir
    setspecies!(tree)
    compute_bipartitions!(tree)
    v = falses(20)
    v[14:20] .= true
    @test getbranch(c, n).bipart.v == v
    v = falses(20)
    v[5:13] .= true
    @test getbranch(d, e).bipart.v == v
    set1 = Set(bipartitions(tree))
    for p in preorder_vector(tree)[2:end]
        origin!(tree, p)
        compute_bipartitions!(tree, dir)
        set2 = Set(bipartitions(tree))
        @test all([x in set2 for x in set2])
        set1 = set2
    end
end


text:
# See https://github.com/JuliaIO/JLD.jl/blob/master/doc/jld.md#custom-serialization
__precompile__()
module PyCallJLD

using PyCall, JLD

const dumps = PyNULL()
const loads = PyNULL()

function __init__()
    pickle = pyimport(PyCall.pyversion.major ≥ 3 ? "pickle" : "cPickle")
    copy!(dumps, pickle["dumps"])
    copy!(loads, pickle["loads"])
end

struct PyObjectSerialization
    repr::Vector{UInt8}
end

function JLD.writeas(pyo::PyObject)
    b = PyCall.PyBuffer(pycall(dumps, PyObject, pyo))
    # We need a `copy` here because the PyBuffer might be GC'ed after we've
    # left this scope, but see
    # https://github.com/JuliaPy/PyCallJLD.jl/pull/3/files/17b052d018f79905baf855b40e440d2cacc171ae#r115525173
    PyObjectSerialization(copy(unsafe_wrap(Array, Ptr{UInt8}(pointer(b)), sizeof(b))))
end

function JLD.readas(pyo_ser::PyObjectSerialization)
    pycall(loads, PyObject,
           PyObject(PyCall.@pycheckn ccall(@pysym(PyCall.PyString_FromStringAndSize),
                                           PyPtr, (Ptr{UInt8}, Int),
                                           pyo_ser.repr, sizeof(pyo_ser.repr))))
end

end # module

ext: jl


text:
using Mocking

cancel_instance_refresh_1_patch = @patch function cancel_instance_refresh(a...; b...)
    return Dict{String, Any}("InstanceRefreshId" => "08b91cf7-8fa6-48af-b6a6-d227f40f1b9b")
end

describe_account_limits_1_patch = @patch function describe_account_limits(a...; b...)
    return Dict{String, Any}("MaxNumberOfLaunchConfigurations" => 100, "MaxNumberOfAutoScalingGroups" => 20, "NumberOfAutoScalingGroups" => 3, "NumberOfLaunchConfigurations" => 5)
end

describe_adjustment_types_1_patch = @patch function describe_adjustment_types(a...; b...)
    return Dict{String, Any}("AdjustmentTypes" => Any[Dict{String, Any}("AdjustmentType" => "ChangeInCapacity"), Dict{String, Any}("AdjustmentType" => "ExactCapcity"), Dict{String, Any}("AdjustmentType" => "PercentChangeInCapacity")])
end

describe_auto_scaling_groups_1_patch = @patch function describe_auto_scaling_groups(a...; b...)
    return Dict{String, Any}("AutoScalingGroups" => Any[Dict{String, Any}("DefaultCooldown" => 300, "EnabledMetrics" => Any[], "MaxSize" => 1, "AvailabilityZones" => Any["us-west-2c"], "CreatedTime" => "2013-08-19T20:53:25.584Z", "DesiredCapacity" => 1, "Tags" => Any[], "AutoScalingGroupName" => "my-auto-scaling-group", "SuspendedProcesses" => Any[], "HealthCheckType" => "EC2", "LaunchConfigurationName" => "my-launch-config", "LoadBalancerNames" => Any[], "AutoScalingGroupARN" => "arn:aws:autoscaling:us-west-2:123456789012:autoScalingGroup:930d940e-891e-4781-a11a-7b0acd480f03:autoScalingGroupName/my-auto-scaling-group", "NewInstancesProtectedFromScaleIn" => false, "HealthCheckGracePeriod" => 300, "TerminationPolicies" => Any["Default"], "VPCZoneIdentifier" => "subnet-12345678", "Instances" => Any[Dict{String, Any}("LifecycleState" => "InService", "AvailabilityZone" => "us-west-2c", "HealthStatus" => "Healthy", "InstanceId" => "i-4ba0837f", "ProtectedFromScaleIn" => false, "LaunchConfigurationName" => "my-launch-config")], "MinSize" => 0)])
end

describe_auto_scaling_instances_1_patch = @patch function describe_auto_scaling_instances(a...; b...)
    return Dict{String, Any}("AutoScalingInstances" => Any[Dict{String, Any}("LifecycleState" => "InService", "AvailabilityZone" => "us-west-2c", "HealthStatus" => "HEALTHY", "InstanceId" => "i-4ba0837f", "AutoScalingGroupName" => "my-auto-scaling-group", "ProtectedFromScaleIn" => false, "LaunchConfigurationName" => "my-launch-config")])
end

describe_auto_scaling_notification_types_1_patch = @patch function describe_auto_scaling_notification_types(a...; b...)
    return Dict{String, Any}("AutoScalingNotificationTypes" => Any["autoscaling:EC2_INSTANCE_LAUNCH", "autoscaling:EC2_INSTANCE_LAUNCH_ERROR", "autoscaling:EC2_INSTANCE_TERMINATE", "autoscaling:EC2_INSTANCE_TERMINATE_ERROR", "autoscaling:TEST_NOTIFICATION"])
end

describe_instance_refreshes_1_patch = @patch function describe_instance_refreshes(a...; b...)
    return Dict{String, Any}("InstanceRefreshes" => Any[Dict{String, Any}("StartTime" => "2020-06-02T18:11:27Z", "InstanceRefreshId" => "08b91cf7-8fa6-48af-b6a6-d227f40f1b9b", "Status" => "InProgress", "AutoScalingGroupName" => "my-auto-scaling-group", "InstancesToUpdate" => 5, "PercentageComplete" => 0), Dict{String, Any}("StartTime" => "2020-06-02T16:43:19Z", "InstanceRefreshId" => "dd7728d0-5bc4-4575-96a3-1b2c52bf8bb1", "EndTime" => "2020-06-02T16:53:37Z", "Status" => "Successful", "AutoScalingGroupName" => "my-auto-scaling-group", "InstancesToUpdate" => 0, "PercentageComplete" => 100)])
end

describe_launch_configurations_1_patch = @patch function describe_launch_configurations(a...; b...)
    return Dict{String, Any}("LaunchConfigurations" => Any[Dict{String, Any}("BlockDeviceMappings" => Any[], "InstanceMonitoring" => Dict{String, Any}("Enabled" => true), "LaunchConfigurationName" => "my-launch-config", "InstanceType" => "t1.micro", "AssociatePublicIpAddress" => true, "SecurityGroups" => Any["sg-67ef0308"], "ImageId" => "ami-043a5034", "EbsOptimized" => false, "LaunchConfigurationARN" => "arn:aws:autoscaling:us-west-2:123456789012:launchConfiguration:98d3b196-4cf9-4e88-8ca1-8547c24ced8b:launchConfigurationName/my-launch-config", "CreatedTime" => "2014-05-07T17:39:28.599Z")])
end

describe_lifecycle_hook_types_1_patch = @patch function describe_lifecycle_hook_types(a...; b...)
    return Dict{String, Any}("LifecycleHookTypes" => Any["autoscaling:EC2_INSTANCE_LAUNCHING", "autoscaling:EC2_INSTANCE_TERMINATING"])
end

describe_lifecycle_hooks_1_patch = @patch function describe_lifecycle_hooks(a...; b...)
    return Dict{String, Any}("LifecycleHooks" => Any[Dict{String, Any}("LifecycleHookName" => "my-lifecycle-hook", "NotificationTargetARN" => "arn:aws:sns:us-west-2:123456789012:my-sns-topic", "DefaultResult" => "ABANDON", "GlobalTimeout" => 172800, "HeartbeatTimeout" => 3600, "AutoScalingGroupName" => "my-auto-scaling-group", "RoleARN" => "arn:aws:iam::123456789012:role/my-auto-scaling-role", "LifecycleTransition" => "autoscaling:EC2_INSTANCE_LAUNCHING")])
end

describe_load_balancer_target_groups_1_patch = @patch function describe_load_balancer_target_groups(a...; b...)
    return Dict{String, Any}("LoadBalancerTargetGroups" => Any[Dict{String, Any}("LoadBalancerTargetGroupARN" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067", "State" => "Added")])
end

describe_load_balancers_1_patch = @patch function describe_load_balancers(a...; b...)
    return Dict{String, Any}("LoadBalancers" => Any[Dict{String, Any}("LoadBalancerName" => "my-load-balancer", "State" => "Added")])
end

describe_metric_collection_types_1_patch = @patch function describe_metric_collection_types(a...; b...)
    return Dict{String, Any}("Granularities" => Any[Dict{String, Any}("Granularity" => "1Minute")], "Metrics" => Any[Dict{String, Any}("Metric" => "GroupMinSize"), Dict{String, Any}("Metric" => "GroupMaxSize"), Dict{String, Any}("Metric" => "GroupDesiredCapacity"), Dict{String, Any}("Metric" => "GroupInServiceInstances"), Dict{String, Any}("Metric" => "GroupPendingInstances"), Dict{String, Any}("Metric" => "GroupTerminatingInstances"), Dict{String, Any}("Metric" => "GroupStandbyInstances"), Dict{String, Any}("Metric" => "GroupTotalInstances")])
end

describe_notification_configurations_1_patch = @patch function describe_notification_configurations(a...; b...)
    return Dict{String, Any}("NotificationConfigurations" => Any[Dict{String, Any}("TopicARN" => "arn:aws:sns:us-west-2:123456789012:my-sns-topic-2", "NotificationType" => "autoscaling:TEST_NOTIFICATION", "AutoScalingGroupName" => "my-auto-scaling-group"), Dict{String, Any}("TopicARN" => "arn:aws:sns:us-west-2:123456789012:my-sns-topic", "NotificationType" => "autoscaling:TEST_NOTIFICATION", "AutoScalingGroupName" => "my-auto-scaling-group")])
end

describe_policies_1_patch = @patch function describe_policies(a...; b...)
    return Dict{String, Any}("ScalingPolicies" => Any[Dict{String, Any}("ScalingAdjustment" => -1, "PolicyName" => "ScaleIn", "Alarms" => Any[], "AutoScalingGroupName" => "my-auto-scaling-group", "PolicyARN" => "arn:aws:autoscaling:us-west-2:123456789012:scalingPolicy:2233f3d7-6290-403b-b632-93c553560106:autoScalingGroupName/my-auto-scaling-group:policyName/ScaleIn", "AdjustmentType" => "ChangeInCapacity"), Dict{String, Any}("MinAdjustmentStep" => 2, "PolicyName" => "ScalePercentChange", "ScalingAdjustment" => 25, "Cooldown" => 60, "Alarms" => Any[], "AutoScalingGroupName" => "my-auto-scaling-group", "PolicyARN" => "arn:aws:autoscaling:us-west-2:123456789012:scalingPolicy:2b435159-cf77-4e89-8c0e-d63b497baad7:autoScalingGroupName/my-auto-scaling-group:policyName/ScalePercentChange", "AdjustmentType" => "PercentChangeInCapacity")])
end

describe_scaling_activities_1_patch = @patch function describe_scaling_activities(a...; b...)
    return Dict{String, Any}("Activities" => Any[Dict{String, Any}("StartTime" => "2013-08-19T20:53:29.930Z", "Description" => "Launching a new EC2 instance: i-4ba0837f", "EndTime" => "2013-08-19T20:54:02Z", "Progress" => 100, "ActivityId" => "f9f2d65b-f1f2-43e7-b46d-d86756459699", "StatusCode" => "Successful", "AutoScalingGroupName" => "my-auto-scaling-group", "Cause" => "At 2013-08-19T20:53:25Z a user request created an AutoScalingGroup changing the desired capacity from 0 to 1.  At 2013-08-19T20:53:29Z an instance was started in response to a difference between desired and actual capacity, increasing the capacity from 0 to 1.", "Details" => "details")])
end

describe_scaling_process_types_1_patch = @patch function describe_scaling_process_types(a...; b...)
    return Dict{String, Any}("Processes" => Any[Dict{String, Any}("ProcessName" => "AZRebalance"), Dict{String, Any}("ProcessName" => "AddToLoadBalancer"), Dict{String, Any}("ProcessName" => "AlarmNotification"), Dict{String, Any}("ProcessName" => "HealthCheck"), Dict{String, Any}("ProcessName" => "Launch"), Dict{String, Any}("ProcessName" => "ReplaceUnhealthy"), Dict{String, Any}("ProcessName" => "ScheduledActions"), Dict{String, Any}("ProcessName" => "Terminate")])
end

describe_scheduled_actions_1_patch = @patch function describe_scheduled_actions(a...; b...)
    return Dict{String, Any}("ScheduledUpdateGroupActions" => Any[Dict{String, Any}("ScheduledActionName" => "my-scheduled-action", "ScheduledActionARN" => "arn:aws:autoscaling:us-west-2:123456789012:scheduledUpdateGroupAction:8e86b655-b2e6-4410-8f29-b4f094d6871c:autoScalingGroupName/my-auto-scaling-group:scheduledActionName/my-scheduled-action", "StartTime" => "2016-12-01T00:30:00Z", "DesiredCapacity" => 4, "Recurrence" => "30 0 1 12 0", "AutoScalingGroupName" => "my-auto-scaling-group", "MaxSize" => 6, "Time" => "2016-12-01T00:30:00Z", "MinSize" => 2)])
end

describe_tags_1_patch = @patch function describe_tags(a...; b...)
    return Dict{String, Any}("Tags" => Any[Dict{String, Any}("Value" => "Research", "ResourceType" => "auto-scaling-group", "ResourceId" => "my-auto-scaling-group", "Key" => "Dept", "PropagateAtLaunch" => true), Dict{String, Any}("Value" => "WebServer", "ResourceType" => "auto-scaling-group", "ResourceId" => "my-auto-scaling-group", "Key" => "Role", "PropagateAtLaunch" => true)])
end

describe_termination_policy_types_1_patch = @patch function describe_termination_policy_types(a...; b...)
    return Dict{String, Any}("TerminationPolicyTypes" => Any["ClosestToNextInstanceHour", "Default", "NewestInstance", "OldestInstance", "OldestLaunchConfiguration"])
end

detach_instances_1_patch = @patch function detach_instances(a...; b...)
    return Dict{String, Any}("Activities" => Any[Dict{String, Any}("StartTime" => "2015-04-12T15:02:16.179Z", "Description" => "Detaching EC2 instance: i-93633f9b", "Progress" => 50, "ActivityId" => "5091cb52-547a-47ce-a236-c9ccbc2cb2c9", "StatusCode" => "InProgress", "AutoScalingGroupName" => "my-auto-scaling-group", "Cause" => "At 2015-04-12T15:02:16Z instance i-93633f9b was detached in response to a user request, shrinking the capacity from 2 to 1.", "Details" => "details")])
end

enter_standby_1_patch = @patch function enter_standby(a...; b...)
    return Dict{String, Any}("Activities" => Any[Dict{String, Any}("StartTime" => "2015-04-12T15:10:23.640Z", "Description" => "Moving EC2 instance to Standby: i-93633f9b", "Progress" => 50, "ActivityId" => "ffa056b4-6ed3-41ba-ae7c-249dfae6eba1", "StatusCode" => "InProgress", "AutoScalingGroupName" => "my-auto-scaling-group", "Cause" => "At 2015-04-12T15:10:23Z instance i-93633f9b was moved to standby in response to a user request, shrinking the capacity from 2 to 1.", "Details" => "details")])
end

exit_standby_1_patch = @patch function exit_standby(a...; b...)
    return Dict{String, Any}("Activities" => Any[Dict{String, Any}("StartTime" => "2015-04-12T15:14:29.886Z", "Description" => "Moving EC2 instance out of Standby: i-93633f9b", "Progress" => 30, "ActivityId" => "142928e1-a2dc-453a-9b24-b85ad6735928", "StatusCode" => "PreInService", "AutoScalingGroupName" => "my-auto-scaling-group", "Cause" => "At 2015-04-12T15:14:29Z instance i-93633f9b was moved out of standby in response to a user request, increasing the capacity from 1 to 2.", "Details" => "details")])
end

put_scaling_policy_1_patch = @patch function put_scaling_policy(a...; b...)
    return Dict{String, Any}("Alarms" => Any[Dict{String, Any}("AlarmName" => "TargetTracking-my-asg-AlarmHigh-fc0e4183-23ac-497e-9992-691c9980c38e", "AlarmARN" => "arn:aws:cloudwatch:us-west-2:123456789012:alarm:TargetTracking-my-asg-AlarmHigh-fc0e4183-23ac-497e-9992-691c9980c38e"), Dict{String, Any}("AlarmName" => "TargetTracking-my-asg-AlarmLow-61a39305-ed0c-47af-bd9e-471a352ee1a2", "AlarmARN" => "arn:aws:cloudwatch:us-west-2:123456789012:alarm:TargetTracking-my-asg-AlarmLow-61a39305-ed0c-47af-bd9e-471a352ee1a2")], "PolicyARN" => "arn:aws:autoscaling:us-west-2:123456789012:scalingPolicy:228f02c2-c665-4bfd-aaac-8b04080bea3c:autoScalingGroupName/my-auto-scaling-group:policyName/alb1000-target-tracking-scaling-policy")
end

start_instance_refresh_1_patch = @patch function start_instance_refresh(a...; b...)
    return Dict{String, Any}("InstanceRefreshId" => "08b91cf7-8fa6-48af-b6a6-d227f40f1b9b")
end


text:
function f_value(F, nEq)

    f = zeros(nEq)
    f[:] = [F]

    return f

end

function f_value(F::Union{Vector{<:Real}, Matrix{<:Real}}, nEq)

    f = zeros(nEq)
    f[:] = F

    return f

end

## Gamma function for vectors ##
Γ(b) = map(gamma, b)

##
# Based on Roberto Garrappa's codes
function DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn_pred, zn_corr, N, METH, METH_fft, Probl)

    # I defined them bc they couldn't sum themselves in a loop (e.g. nxi += r)
    nxi = 0
    nxf = 0
    nyi = 0
    nyf = 0

    nxi = nx0
    nxf = nx0 + L * r - 1
    nyi = ny0
    nyf = ny0 + L * r - 1
    is = 1

    s_nxi = Int64.(zeros(L))
    s_nxf = Int64.(zeros(L))
    s_nyi = Int64.(zeros(L))
    s_nyf = Int64.(zeros(L))

    s_nxi[is] = nxi
    s_nxf[is] = nxf
    s_nyi[is] = nyi
    s_nyf[is] = nyf

    i_triangolo = 0
    stop = 0
    stop = false # We have to define stop (e.g. by zeros) before the loop to keep changes inside the loop!

    while !stop

        stop = ((nxi + r - 1 == nx0 + L * r - 1) || (nxi + r - 1 >= Nr - 1)) # It stops when current triangle ends at the end of the square

        zn_pred, zn_corr = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, METH_fft, Probl)

        y, fy = Triangolo(nxi, nxi + r - 1, t, y, fy, zn_pred, zn_corr, N, METH, Probl)
        i_triangolo += 1

        if !stop

            if nxi + r - 1 == nxf   # The triangle ends where the square ends -> you level down

                i_Delta = Int64.(ff[i_triangolo])
                Delta = i_Delta * r

                nxi = s_nxf[is] + 1
                nxf = s_nxf[is]  + Delta
                nyi = s_nxf[is] - Delta + 1
                nyf = s_nxf[is]

                s_nxi[is] = nxi
                s_nxf[is] = nxf
                s_nyi[is] = nyi
                s_nyf[is] = nyf

            else # The triangle ends before the square -> a square is made next to it

                nxi += r
                nxf = nxi + r - 1
                nyi = nyf + 1
                nyf += r

                is += 1

                s_nxi[is] = nxi
                s_nxf[is] = nxf
                s_nyi[is] = nyi
                s_nyf[is] = nyf

            end

        end

    end

    return y, fy

end

##
# Based on Roberto Garrappa's codes
function Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, METH_fft, Probl)

    coef_end = nxf - nyi + 1
    i_fft = Int64.(log2(coef_end / METH.r))

    funz_beg = nyi + 1
    funz_end = nyf + 1

    Nnxf = min(N, nxf)

 # Evaluation convolution segment for the predictor
    vett_funz = zeros(Probl.problem_size, coef_end)
    vett_funz[:, 1:funz_end - funz_beg + 1] = fy[:, funz_beg:funz_end]
    vett_funz_fft = fft(vett_funz, 2)
    zzn_pred = zeros(Probl.problem_size, coef_end)

    for i in 1:Probl.problem_size

        i_β = min(Probl.β_length, i)
        Z = METH_fft.bn_fft[i_β,METH_fft.index_fft[1, i_fft]:METH_fft.index_fft[2, i_fft]] .* vett_funz_fft[i, :]
        zzn_pred[i, :] = real(ifft(Z))

    end

    zzn_pred = zzn_pred[:, nxf - nyf:end - 1]
    zn_pred[:, nxi + 1:Nnxf + 1] += zzn_pred[:, 1:Nnxf - nxi + 1]

    # Evaluation convolution segment for the corrector
    if METH.μ > 0

        if nyi == 0 # Evaluation of the lowest square

            vett_funz = zeros(Probl.problem_size, coef_end)
            vett_funz[:, 1:funz_end - funz_beg + 1] = [zeros(Probl.problem_size, 1) fy[:, funz_beg + 1:funz_end]]
            vett_funz_fft = fft(vett_funz, 2)

        end

        zzn_corr = zeros(Probl.problem_size, coef_end)

        for i in 1:Probl.problem_size

            i_β = min(Probl.β_length, i)
            Z = METH_fft.an_fft[i_β, METH_fft.index_fft[1, i_fft]:METH_fft.index_fft[2, i_fft]] .* vett_funz_fft[i, :]
            zzn_corr[i, :] = real(ifft(Z))

        end

        zzn_corr = zzn_corr[:, nxf - nyf + 1:end]
        zn_corr[:, nxi + 1:Nnxf + 1] += zzn_corr[:, 1:Nnxf - nxi + 1]

    else

        zn_corr = 0

    end

    return zn_pred, zn_corr

end

##
# Based on Roberto Garrappa's codes
function Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, Probl)

    for n in nxi:min(N, nxf)

        # Evaluation of the predictor
        Φ = zeros(Probl.problem_size, 1)

        j_beg = ifelse(nxi == 1, 0, nxi) # if nxi == 1 -> first triangle, else -> any other triangle

        for j in j_beg:n - 1

            Φ += METH.bn[1:Probl.β_length, n - j] .* fy[:, j + 1]

        end

        St = taylor_expansion(t[n + 1], Probl.ic)
        y_pred = St .+ METH.hα1 .* (zn_pred[:, n + 1] .+ Φ)
        f_pred = f_value(Probl.f_fun(t[n+1], y_pred, Probl.param...), Probl.problem_size)

        # Evaluation of the corrector
        if METH.μ == 0

            y[:, n + 1] = y_pred
            fy[:, n + 1] = f_pred

        else

            j_beg = nxi
            Φ = zeros(Probl.problem_size, 1)

            for j in j_beg:n - 1

                Φ += METH.an[1:Probl.β_length, n - j + 1] .* fy[:, j + 1]

            end

            Φ_n = St + METH.hα2 .* (METH.a0[1:Probl.β_length, n + 1] .* fy[:, 1] .+ zn_corr[:, n + 1] .+ Φ)

            yn0 = y_pred
            fn0 = f_pred

            stop = 0
            stop = false # it is defined for counting in the following loop while
            mu_it = 0

            yn1 = zeros(Probl.problem_size, 1)
            fn1 = zeros(Probl.problem_size, 1)

            while !stop

                yn1 = Φ_n + METH.hα2 .* fn0
                mu_it += 1

                if METH.StopIt == "Convergence"

                    stop = norm(yn1 - yn0, Inf) < METH.μTol

                    if (mu_it > METH.itmax && !stop)

                        stop = true

                    end

                else

                    stop = (mu_it == METH.μ)

                end

                fn1 = f_value(Probl.f_fun(t[n+1], yn1, Probl.param...),Probl.problem_size)
                yn0 = yn1
                fn0 = fn1

            end

            y[:, n + 1] = yn1
            fy[:, n + 1] = fn1

        end

    end

    return y, fy

end

##
# Based on Roberto Garrappa's code (if we remove our previous code without FFT, then we can changed this name to taylor_expansion)
function taylor_expansion(t, ic)

    ys = zeros(size(ic.y0, 1), 1)

    for k in 1:maximum(ic.m_β)

        if length(ic.m_β) == 1

            ys = ys .+ (t - ic.t0).^(k - 1) ./ ic.m_β_factorial[k] .* ic.y0[:, k]

        else

            i_β = findall(k .<= ic.m_β)
            ys[i_β, 1] += (t - ic.t0)^(k - 1) * ic.y0[i_β,k] ./ ic.m_β_factorial[i_β, k]

        end

    end

    return ys

end


text:
#FT_DEVICE (DWORD)
if ~isdefined(:FT_DEVICE_232BM)
  export FT_DEVICE_232BM, FT_DEVICE_232AM, FT_DEVICE_100AX
  export FT_DEVICE_UNKNOWN, FT_DEVICE_2232C, FT_DEVICE_232R
  export FT_DEVICE_2232H, FT_DEVICE_4232H, FT_DEVICE_232H
  export FT_DEVICE_X_SERIES
  const FT_DEVICE_232BM = 0 # PFT_EEPROM_232B
  const FT_DEVICE_232AM = 1 # discontinued device
  const FT_DEVICE_100AX = 2 # discontinued device
  const FT_DEVICE_UNKNOWN = 3
  const FT_DEVICE_2232C = 4 # PFT_EEPROM_2232
  const FT_DEVICE_232R = 5 # PFT_EEPROM_232R
  const FT_DEVICE_2232H = 6 # PFT_EEPROM_2232H
  const FT_DEVICE_4232H = 7 # PFT_EEPROM_4232H
  const FT_DEVICE_232H = 8 # PFT_EEPROM_232H
  const FT_DEVICE_X_SERIES = 9 # PFT_EEPROM_X_SERIES
  const FT_DEVICE_4222H_0 = 10 #
  const FT_DEVICE_4222H_1_2 = 11 #
  const FT_DEVICE_4222H_3 = 12 #
  const FT_DEVICE_4222_PROG =13 #

  export ftdevicetypedict
  ftdevicetypedict = Dict(
    FT_DEVICE_232BM => "232BM",
    FT_DEVICE_232AM => "232AM",
    FT_DEVICE_100AX => "100AX",
    FT_DEVICE_UNKNOWN => "UNKNOWN",
    FT_DEVICE_2232C => "2232C",
    FT_DEVICE_232R => "232R",
    FT_DEVICE_2232H => "2232H",
    FT_DEVICE_4232H => "4232H",
    FT_DEVICE_232H => "232H",
    FT_DEVICE_X_SERIES => "X_SERIES",
    FT_DEVICE_4222H_0 => "FT_DEVICE_4222H_0",
    FT_DEVICE_4222H_1_2 => "FT_DEVICE_4222H_1_2",
    FT_DEVICE_4222H_3 => "FT_DEVICE_4222H_3",
    FT_DEVICE_4222_PROG => "FT_DEVICE_4222_PROG")
end


text:
using FactorGraph, Random

# generate a random with 5 nodes (factors) and 10 edges (variables)
Random.seed!(2); rg = random_fg(5, 10)

# remove the first two edges from this incidence list, notice these labels (numbers) of edges are dynamicaly assigned.
rg = rg |> rem_edge(1) |> rem_edge(1)
# remove the first vertex from this incidence list
rg = rg |> rem_vertex(2)

# eliminate a variable (sum over an edge)
rg = rg |> contract(3) |> contract(1)

# the resulting tensors
rg.tensors .|> size

# the topology of the incidence list of the resulting V-3, E-4 graph.
rg.legmap


text: abstract type AbstractImmutableSimplex{D, T} <: AbstractSimplex{D, T} end


text:
function k_tr{T}(x::T)
    if(isnan(x) || abs(x)<= one(T))
        return one(Float64)
    else
        return zero(Float64)
    end
end

function k_bt{T}(x::T)
    if isnan(x)
        one(Float64)
    end
    float(max(one(T)-abs(x), zero(T)))
end

function k_pr{T}(x::T)
    if isnan(x)
        one(Float64)
    end
    ax = abs(x)
    if(ax > one(T))
        zero(Float64)
    elseif ax <= .5
        float(1 - 6 * x^2 + 6 * ax^3)
    else
        float(2 * (1-ax)^3)
    end
end

function k_qs{T <: Number}(x::T)
    if isnan(x)
        one(Float64)
    end
    if(isequal(x, zero(eltype(x))))
        one(Float64)
    else
        return (25/(12*π²*x^2))*(sin(sixπ*x/5)/(sixπ*x/5)-cos(sixπ*x/5))
    end
end

function k_th{T <: Number}(x::T)
    if isnan(x)
        one(Float64)
    end
    ax = abs(x)
    if(ax < one(T))
        (1 + cos(π*x))/2
    else
        zero(Float64)
    end
end

##############################################################################
##
## Optimal band-width
##
##############################################################################
type TruncatedKernel <: HAC
  kernel::Function
  bw::Function
end

type BartlettKernel <: HAC
  kernel::Function
  bw::Function
end

type ParzenKernel <: HAC
  kernel::Function
  bw::Function
end

type TukeyHanningKernel <: HAC
  kernel::Function
  bw::Function
end

type QuadraticSpectralKernel <: HAC
  kernel::Function
  bw::Function
end

type VARHAC <: HAC
    imax::Int64
    ilag::Int64
    imodel::Int64
end

typealias TRK TruncatedKernel
typealias BTK BartlettKernel
typealias PRK ParzenKernel
typealias THK TukeyHanningKernel
typealias QSK QuadraticSpectralKernel

TruncatedKernel()                   = TRK(k_tr, optimalbw_ar_one)
BartlettKernel()                    = BTK(k_bt, optimalbw_ar_one)
ParzenKernel()                      = PRK(k_pr, optimalbw_ar_one)
TukeyHanningKernel()                = THK(k_th, optimalbw_ar_one)
QuadraticSpectralKernel()           = QSK(k_qs, optimalbw_ar_one)

TruncatedKernel(bw::Number)         = TRK(k_tr, (x, k) -> float(bw))
BartlettKernel(bw::Number)          = BTK(k_bt, (x, k) -> float(bw))
ParzenKernel(bw::Number)            = PRK(k_pr, (x, k) -> float(bw))
TukeyHanningKernel(bw::Number)      = THK(k_th, (x, k) -> float(bw))
QuadraticSpectralKernel(bw::Number) = QSK(k_qs, (x, k) -> float(bw))

VARHAC()                            = VARHAC(2, 2, 1)
VARHAC(imax::Int64)                 = VARHAC(imax, 2, 1)

function bandwidth(k::HAC, X::AbstractMatrix)
  return floor(k.bw(X, k))
end

function bandwidth(k::QuadraticSpectralKernel, X::AbstractMatrix)
    return k.bw(X, k)
end

kernel(k::HAC, x::Real) = k.kernel(x)

function Γ(X::AbstractMatrix, j::Int64)
    T, p = size(X)
    Q = zeros(eltype(X), p, p)
    if j>=0
        for h=1:p, s = 1:h
            for t = j+1:T
                @inbounds Q[s, h] = Q[s, h] + X[t, s]*X[t-j, h]
            end
        end
    else
        for h=1:p, s = 1:h
            for t = -j+1:T
                @inbounds Q[s,h] = Q[s ,h] + X[t+j, s]*X[t,h]
            end
        end
    end
    return Q
end

vcov(X::AbstractMatrix, k::VARHAC) = varhac(X, k.imax, k.ilag, k.imodel)

function vcov(X::AbstractMatrix, k::HAC; prewhite=true)
    n, p = size(X)
    !prewhite || ((X, D) = pre_white(X))
    bw = bandwidth(k, X)
    Q  = zeros(eltype(X), p, p)
    for j=-bw:bw
        Base.BLAS.axpy!(kernel(k, j/bw), Γ(X, @compat Int(j)), Q)
    end
    Base.LinAlg.copytri!(Q, 'U')
    if prewhite
        Q[:] = D*Q*D'
    end
    return scale!(Q, 1/n)
end

function vcov(X::AbstractMatrix, k::QuadraticSpectralKernel; prewhite=true)
    n, p = size(X)
    !prewhite || ((X, D) = pre_white(X))
    bw = bandwidth(k, X)
    Q = zeros(eltype(X), p, p)
    for j=-n:n
        Base.BLAS.axpy!(kernel(k, j/bw), Γ(X, @compat Int(j)), Q)
    end
    Base.LinAlg.copytri!(Q, 'U')
    if prewhite
        Q[:] = D*Q*D'
    end
    return scale!(Q, 1/n)
end

vcov(x::DataFrameRegressionModel, k::HAC; args...) = vcov(x.model, k; args...)

function vcov(ll::LinPredModel, k::HAC; args...)
    B = meat(ll, k; args...)
    A = bread(ll)
    scale!(A*B*A, 1/nobs(ll))
end

function meat(l::LinPredModel,  k::HAC; args...)
    u = wrkresidwts(l.rr)
    X = ModelMatrix(l)
    z = X.*u
    vcov(z, k; args...)
end


text:
mutable struct SpglibSpacegroupType
    number::Cint
    international_short::NTuple{11, UInt8}
    international_full::NTuple{20, UInt8}
    international::NTuple{32, UInt8}
    schoenflies::NTuple{7, UInt8}
    hall_symbol::NTuple{17, UInt8}
    choice::NTuple{6, UInt8}
    pointgroup_international::NTuple{6, UInt8}
    pointgroup_schoenflies::NTuple{4, UInt8}
    arithmetic_crystal_class_number::Cint
    arithmetic_crystal_class_symbol::NTuple{7, UInt8}
end

function spg_get_spacegroup_type(hall_number::Int64)
    hall_number = Base.convert(Cint, hall_number)

    return ccall((:spg_get_spacegroup_type, libsymspg), SpglibSpacegroupType,
                (Cint, ),
                hall_number)
end


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_instances_by_tag("ClusterName", ec2_clustername)
println("Stopping EC2 nodes : ", instances)

ec2_stop(instances)
ec2_show_status(instances)


text:
export
    MMCIFDict,
    readmultimmcif,
    writemmcif,
    writemultimmcif

# mmCIF special characters
const quotechars = Set(['\'', '\"'])
const whitespacechars = Set([' ', '\t'])
const missingvals = Set([".", "?"])
const specialchars = Set(['_', '#', '\$', '[', ']', ';'])
# These are technically case-insensitive at any letter but just checking the
#   all-upper and all-lower forms is faster; this also appears in code below
const specialwords = Set(["loop_", "LOOP_", "stop_", "STOP_",
                            "global_", "GLOBAL_"])

# If certain entries should have a certain order of keys, that is specified here
const mmciforder = Dict(
    "_atom_site"=> [
        "group_PDB",
        "id",
        "type_symbol",
        "label_atom_id",
        "label_alt_id",
        "label_comp_id",
        "label_asym_id",
        "label_entity_id",
        "label_seq_id",
        "pdbx_PDB_ins_code",
        "Cartn_x",
        "Cartn_y",
        "Cartn_z",
        "occupancy",
        "B_iso_or_equiv",
        "pdbx_formal_charge",
        "auth_seq_id",
        "auth_comp_id",
        "auth_asym_id",
        "auth_atom_id",
        "pdbx_PDB_model_num",
    ]
)

"""
    MMCIFDict(filepath; gzip=false)
    MMCIFDict(io; gzip=false)
    MMCIFDict()

A macromolecular Crystallographic Information File (mmCIF) dictionary.

Can be accessed using similar functions to a standard `Dict`.
Keys are field names as a `String` and values are always `Vector{String}`, even
for multiple components or numerical data.
To directly access the underlying dictionary of `MMCIFDict` `d`, use
`d.dict`.
Call `MMCIFDict` with a filepath or stream to read the dictionary from that
source.
The keyword argument `gzip` (default `false`) determines if the input is gzipped.
"""
struct MMCIFDict <: AbstractDict{String, Vector{String}}
    dict::Dict{String, Vector{String}}
end

MMCIFDict() = MMCIFDict(Dict())

Base.getindex(mmcif_dict::MMCIFDict, field::AbstractString) = mmcif_dict.dict[field]

function Base.setindex!(mmcif_dict::MMCIFDict,
                    val::Vector{String},
                    field::AbstractString)
    mmcif_dict.dict[field] = val
    return mmcif_dict
end

Base.keys(mmcif_dict::MMCIFDict) = keys(mmcif_dict.dict)
Base.values(mmcif_dict::MMCIFDict) = values(mmcif_dict.dict)
Base.haskey(mmcif_dict::MMCIFDict, key) = haskey(mmcif_dict.dict, key)
Base.get(mmcif_dict::MMCIFDict, key, default) = get(mmcif_dict.dict, key, default)
Base.length(mmcif_dict::MMCIFDict) = length(mmcif_dict.dict)
Base.iterate(mmcif_dict::MMCIFDict) = iterate(mmcif_dict.dict)
Base.iterate(mmcif_dict::MMCIFDict, i) = iterate(mmcif_dict.dict, i)

function Base.show(io::IO, mmcif_dict::MMCIFDict)
    print(io, "mmCIF dictionary with $(length(mmcif_dict)) fields")
end

# Split a mmCIF line into tokens
# See https://www.iucr.org/resources/cif/spec/version1.1/cifsyntax for syntax
function splitline(s::AbstractString)
    tokens = String[]
    in_token = false
    # Quote character of the currently open quote, or ' ' if no quote open
    quote_open_char = ' '
    start_i = 0
    for (i, c) in enumerate(s)
        if c in whitespacechars
            if in_token && quote_open_char == ' '
                in_token = false
                push!(tokens, s[start_i:(i - 1)])
            end
        elseif c in quotechars
            if quote_open_char == ' '
                if in_token
                    throw(ArgumentError("Opening quote in middle of word: $s"))
                end
                quote_open_char = c
                in_token = true
                start_i = i + 1
            elseif c == quote_open_char && (i == length(s) || s[i + 1] in whitespacechars)
                quote_open_char = ' '
                in_token = false
                push!(tokens, s[start_i:(i - 1)])
            end
        elseif c == '#' && !in_token
            return tokens
        elseif !in_token
            in_token = true
            start_i = i
        end
    end
    if in_token
        push!(tokens, s[start_i:end])
    end
    if quote_open_char != ' '
        throw(ArgumentError("Line ended with quote open: $s"))
    end
    return tokens
end

# Get tokens from a mmCIF file
function tokenizecif(f::IO)
    tokens = String[]
    for line in eachline(f)
        if startswith(line, "#")
            continue
        elseif startswith(line, ";")
            token_buffer = [rstrip(line[2:end])]
            for inner_line in eachline(f)
                inner_strip = rstrip(inner_line)
                if inner_strip == ";"
                    break
                end
                push!(token_buffer, inner_strip)
            end
            push!(tokens, join(token_buffer, "\n"))
        else
            append!(tokens, splitline(line))
        end
    end
    return tokens
end

# Get tokens from a mmCIF file corresponding to atom site records only
# This will fail if there is only a single atom record in the file
#   and it is not in the loop format
function tokenizecifstructure(f::IO)
    tokens = String[]
    reading = false
    in_keys = true
    for line in eachline(f)
        if (!reading && !startswith(line, "_atom_site.")) || startswith(line, "#")
            continue
        elseif startswith(line, "_atom_site.")
            reading = true
            push!(tokens, rstrip(line))
        elseif startswith(line, ";")
            in_keys = false
            token_buffer = [rstrip(line[2:end])]
            for inner_line in eachline(f)
                inner_strip = rstrip(inner_line)
                if inner_strip == ";"
                    break
                end
                push!(token_buffer, inner_strip)
            end
            push!(tokens, join(token_buffer, "\n"))
        elseif (!in_keys && startswith(line, "_")) || startswith(line, "loop_") ||
                    startswith(line, "LOOP_")
            break
        else
            in_keys = false
            append!(tokens, splitline(line))
        end
    end
    length(tokens) > 0 ? pushfirst!(tokens, "loop_") : nothing
    return tokens
end

function MMCIFDict(mmcif_filepath::AbstractString; gzip::Bool=false)
    open(mmcif_filepath) do f
        MMCIFDict(f; gzip=gzip)
    end
end

# Read a mmCIF file into a MMCIFDict
function MMCIFDict(f::IO; gzip::Bool=false)
    mmcif_dict = MMCIFDict()
    if gzip
        gz = GzipDecompressorStream(f)
        tokens = tokenizecif(gz)
        close(gz)
    else
        tokens = tokenizecif(f)
    end
    # Data label token is read first
    if length(tokens) == 0
        return mmcif_dict
    end
    data_token = first(tokens)
    mmcif_dict[data_token[1:5]] = [data_token[6:end]]
    return populatedict!(mmcif_dict, tokens[2:end])
end

# Add tokens to a mmCIF dictionary
function populatedict!(mmcif_dict::MMCIFDict, tokens::Vector{String})
    key = ""
    keys = String[]
    loop_flag = false
    i = 0 # Value counter
    n = 0 # Key counter
    for token in tokens
        if token == "loop_" || token == "LOOP_"
            loop_flag = true
            keys = String[]
            i = 0
            n = 0
            continue
        elseif loop_flag
            # The second condition checks we are in the first column
            # Some mmCIF files (e.g. 4q9r) have values in later columns
            #   starting with an underscore and we don't want to read
            #   these as keys
            if startswith(token, "_") && (n == 0 || i % n == 0)
                if i > 0
                    loop_flag = false
                else
                    mmcif_dict[token] = String[]
                    push!(keys, token)
                    n += 1
                    continue
                end
            else
                try
                    push!(mmcif_dict[keys[i % n + 1]], token)
                catch ex
                    # A zero division error means we have not found any keys
                    if isa(ex, DivideError)
                        throw(ArgumentError("Loop keys not found, token: \"$token\""))
                    else
                        rethrow()
                    end
                end
                i += 1
                continue
            end
        end
        if key == ""
            key = token
        else
            # Single values are read in as an array for consistency
            mmcif_dict[key] = [token]
            key = ""
        end
    end
    return mmcif_dict
end

function Base.read(input::IO,
            ::Type{MMCIF};
            structure_name::AbstractString="",
            remove_disorder::Bool=false,
            read_std_atoms::Bool=true,
            read_het_atoms::Bool=true,
            gzip::Bool=false)
    mmcif_dict = MMCIFDict()
    if gzip
        gz = GzipDecompressorStream(input)
        tokens = tokenizecifstructure(gz)
        close(gz)
    else
        tokens = tokenizecifstructure(input)
    end
    populatedict!(mmcif_dict, tokens)
    ProteinStructure(mmcif_dict;
                     structure_name=structure_name,
                     remove_disorder=remove_disorder,
                     read_std_atoms=read_std_atoms,
                     read_het_atoms=read_het_atoms)
end

function ProteinStructure(mmcif_dict::MMCIFDict;
            structure_name::AbstractString="",
            remove_disorder::Bool=false,
            read_std_atoms::Bool=true,
            read_het_atoms::Bool=true)
    # Define ProteinStructure and add to it incrementally
    struc = ProteinStructure(structure_name)
    if haskey(mmcif_dict, "_atom_site.id")
        for i in 1:length(mmcif_dict["_atom_site.id"])
            if (read_std_atoms && mmcif_dict["_atom_site.group_PDB"][i] == "ATOM") ||
                    (read_het_atoms && mmcif_dict["_atom_site.group_PDB"][i] == "HETATM")
                model_n = parse(Int, mmcif_dict["_atom_site.pdbx_PDB_model_num"][i])
                # Create model if required
                if !haskey(models(struc), model_n)
                    struc[model_n] = Model(model_n, struc)
                end
                unsafe_addatomtomodel!(
                    struc[model_n],
                    AtomRecord(mmcif_dict, i),
                    remove_disorder=remove_disorder)
            end
        end
        # Generate lists for iteration
        fixlists!(struc)
    end
    return struc
end

# Constructor from mmCIF ATOM/HETATM line
AtomRecord(d::MMCIFDict, i::Integer) = AtomRecord(
    d["_atom_site.group_PDB"][i] == "HETATM",
    parse(Int, d["_atom_site.id"][i]),
    d["_atom_site.auth_atom_id"][i],
    d["_atom_site.label_alt_id"][i] in missingvals ? ' ' : d["_atom_site.label_alt_id"][i][1],
    d["_atom_site.auth_comp_id"][i],
    d["_atom_site.auth_asym_id"][i],
    parse(Int, d["_atom_site.auth_seq_id"][i]),
    d["_atom_site.pdbx_PDB_ins_code"][i] in missingvals ? ' ' : d["_atom_site.pdbx_PDB_ins_code"][i][1],
    [
        parse(Float64, d["_atom_site.Cartn_x"][i]),
        parse(Float64, d["_atom_site.Cartn_y"][i]),
        parse(Float64, d["_atom_site.Cartn_z"][i])
    ],
    d["_atom_site.occupancy"][i] in missingvals ? 1.0 : parse(Float64, d["_atom_site.occupancy"][i]),
    d["_atom_site.B_iso_or_equiv"][i] in missingvals ? 0.0 : parse(Float64, d["_atom_site.B_iso_or_equiv"][i]),
    d["_atom_site.type_symbol"][i] in missingvals ? "  " : d["_atom_site.type_symbol"][i],
    d["_atom_site.pdbx_formal_charge"][i] in missingvals ? "  " : d["_atom_site.pdbx_formal_charge"][i],
)

# Format a mmCIF data value by enclosing with quotes or semicolon lines where
#   appropriate. See
#   https://www.iucr.org/resources/cif/spec/version1.1/cifsyntax for syntax.
function formatmmcifcol(val::AbstractString, col_width::Integer=length(val))
    # If there is a newline or quotes cannot be contained, use semicolon
    #   and newline construct
    if requiresnewline(val)
        return "\n;$val\n;\n"
    elseif requiresquote(val)
        # Choose quote character
        if occursin("' ", val)
            return rpad("\"$val\"", col_width)
        else
            return rpad("'$val'", col_width)
        end
    # Safe to not quote
    # Numbers must not be quoted
    else
        return rpad(val, col_width)
    end
end

function requiresnewline(val)
    return occursin("\n", val) || (occursin("' ", val) && occursin("\" ", val))
end

function requiresquote(val)
    return occursin(" ", val) || occursin("'", val) || occursin("\"", val) ||
        val[1] in specialchars || startswith(val, "data_") ||
        startswith(val, "DATA_") || startswith(val, "save_") ||
        startswith(val, "SAVE_") || val in specialwords
end

"""
    writemmcif(output, element, atom_selectors...; gzip=false)
    writemmcif(output, mmcif_dict; gzip=false)

Write a `StructuralElementOrList` or a `MMCIFDict` to a mmCIF format file or
output stream.

Atom selector functions can be given as additional arguments - only atoms
that return `true` from all the functions are retained.
The keyword argument `expand_disordered` (default `true`) determines whether to
return all copies of disordered residues and atoms.
The keyword argument `gzip` (default `false`) determines if the output is gzipped.
"""
function writemmcif(filepath::AbstractString, mmcif_dict::MMCIFDict; gzip::Bool=false)
    open(filepath, "w") do output
        writemmcif(output, mmcif_dict; gzip=gzip)
    end
end

function writemmcif(output::IO, mmcif_dict::MMCIFDict; gzip::Bool=false)
    if gzip
        output = GzipCompressorStream(output)
    end
    # Form dictionary where key is first part of mmCIF key and value is list
    #   of corresponding second parts
    key_lists = Dict{String, Vector{String}}()
    data_val = String[]
    for key in keys(mmcif_dict)
        if key == "data_" || key == "DATA_"
            data_val = mmcif_dict[key]
        else
            s = split(key, ".")
            if length(s) == 2
                if s[1] in keys(key_lists)
                    push!(key_lists[s[1]], s[2])
                else
                    key_lists[s[1]] = [s[2]]
                end
            else
                throw(ArgumentError("Invalid key in mmCIF dictionary: $key"))
            end
        end
    end

    # Re-order lists if an order has been specified
    # Not all elements from the specified order are necessarily present
    for (key, key_list) in key_lists
        if key in keys(mmciforder)
            inds = Int[]
            for i in key_list
                f = something(findfirst(isequal(i), mmciforder[key]), 0)
                if f == 0
                    # Unrecognised key - add at end
                    f = length(mmciforder[key]) + 1
                end
                push!(inds, f)
            end
            key_lists[key] = [k for (_, k) in sort(collect(zip(inds, key_list)))]
        end
    end

    # Write out top data_ line
    if data_val != String[]
        println(output, "data_$(first(data_val))\n#")
    else
        println(output, "data_unknown\n#")
    end

    for (key, key_list) in key_lists
        # Pick a sample mmCIF value
        sample_val = mmcif_dict["$key.$(first(key_list))"]
        n_vals = length(sample_val)
        # Check the mmCIF dictionary has consistent list sizes
        for i in key_list
            if length(mmcif_dict["$key.$i"]) != n_vals
                throw(ArgumentError("Inconsistent list sizes in mmCIF dictionary: $key.$i"))
            end
        end
        # If the value has a single component, write as key-value pairs
        if length(sample_val) == 1
            # Find the maximum key length
            m = maximum(length.(key_list))
            for i in key_list
                println(output, "$(rpad("$key.$i", length(key) + m + 4))$(formatmmcifcol(first(mmcif_dict["$key.$i"])))")
            end
        # If the value has more than one component, write as keys then a value table
        else
            println(output, "loop_")
            col_widths = Dict{String, Int}()
            # Write keys and find max widths for each set of values
            # Technically the max of the sum of the column widths is 2048
            for i in key_list
                println(output, "$key.$i")
                col_widths[i] = 0
                for val in mmcif_dict["$key.$i"]
                    len_val = length(val)
                    # If the value requires quoting it will add 2 characters
                    if requiresquote(val) && !requiresnewline(val)
                        len_val += 2
                    end
                    if len_val > col_widths[i]
                        col_widths[i] = len_val
                    end
                end
            end

            # Write the values as rows
            for i in 1:n_vals
                for col in key_list
                    print(output, formatmmcifcol(mmcif_dict["$key.$col"][i], col_widths[col] + 1))
                end
                println(output)
            end
        end
        println(output, "#")
    end
    if gzip
        close(output)
    end
end

function writemmcif(filepath::AbstractString,
                el::StructuralElementOrList,
                atom_selectors::Function...;
                expand_disordered::Bool=true, gzip::Bool=false)
    open(filepath, "w") do output
        writemmcif(output, el, atom_selectors...;
                   expand_disordered=expand_disordered, gzip=gzip)
    end
end

function writemmcif(output::IO,
                el::StructuralElementOrList,
                atom_selectors::Function...;
                expand_disordered::Bool=true, gzip::Bool=false)
    # Ensure multiple models get written out correctly
    loop_el = isa(el, ProteinStructure) ? collectmodels(el) : el
    ats = collectatoms(loop_el, atom_selectors...;
                                        expand_disordered=expand_disordered)
    if length(ats) > 0
        # Create an empty dictionary and add atoms one at a time
        atom_dict = Dict{String, Vector{String}}(
                ["_atom_site.$i"=> String[] for i in mmciforder["_atom_site"]])
        atom_dict["data_"] = [structurename(first(el))]
    else
        # If we are not writing any atoms, don't write the dictionary keys
        atom_dict = Dict{String, Vector{String}}(["data_"=> ["unknown"]])
    end

    for at in ats
        appendatom!(atom_dict, at, string(modelnumber(at)),
            strip(chainid(at)) == "" ? "." : chainid(at),
            string(resnumber(at)), resname(at),
            ishetero(at) ? "HETATM" : "ATOM")
    end

    # Now the MMCIFDict has been generated, write it out to the file
    return writemmcif(output, MMCIFDict(atom_dict); gzip=gzip)
end

# Add an atom record to a growing atom dictionary
function appendatom!(atom_dict, at, model_n, chain_id, res_n, res_name, het)
    push!(atom_dict["_atom_site.group_PDB"], het)
    push!(atom_dict["_atom_site.id"], string(serial(at)))
    push!(atom_dict["_atom_site.type_symbol"], length(element(at)) == 0 ? "?" : element(at))
    push!(atom_dict["_atom_site.label_atom_id"], atomname(at))
    push!(atom_dict["_atom_site.label_alt_id"], altlocid(at) == ' ' ? "." : string(altlocid(at)))
    push!(atom_dict["_atom_site.label_comp_id"], res_name)
    push!(atom_dict["_atom_site.label_asym_id"], "?")
    push!(atom_dict["_atom_site.label_entity_id"], "?")
    push!(atom_dict["_atom_site.label_seq_id"], "?")
    push!(atom_dict["_atom_site.pdbx_PDB_ins_code"], inscode(at) == ' ' ? "?" : string(inscode(at)))
    push!(atom_dict["_atom_site.Cartn_x"], pyfmt(coordspec, round(x(at), digits=3)))
    push!(atom_dict["_atom_site.Cartn_y"], pyfmt(coordspec, round(y(at), digits=3)))
    push!(atom_dict["_atom_site.Cartn_z"], pyfmt(coordspec, round(z(at), digits=3)))
    push!(atom_dict["_atom_site.occupancy"], pyfmt(floatspec, occupancy(at)))
    push!(atom_dict["_atom_site.B_iso_or_equiv"], pyfmt(floatspec, tempfactor(at)))
    push!(atom_dict["_atom_site.pdbx_formal_charge"], length(charge(at)) == 0 ? "?" : charge(at))
    push!(atom_dict["_atom_site.auth_seq_id"], res_n)
    push!(atom_dict["_atom_site.auth_comp_id"], res_name)
    push!(atom_dict["_atom_site.auth_asym_id"], chain_id)
    push!(atom_dict["_atom_site.auth_atom_id"], atomname(at))
    push!(atom_dict["_atom_site.pdbx_PDB_model_num"], model_n)
end

"""
    readmultimmcif(filepath; gzip=false)
    readmultimmcif(io; gzip=false)

Read multiple `MMCIFDict`s from a filepath or stream.
Each `MMCIFDict` in the returned `Dict{String, MMCIFDict}` corresponds to an mmCIF data block from
the input.
An example of such a file is the chemical component dictionary from the Protein Data Bank.
The keyword argument `gzip` (default `false`) determines if the input is gzipped.
"""
function readmultimmcif(filepath::AbstractString; gzip::Bool=false)
    open(filepath) do io
        readmultimmcif(io; gzip=gzip)
    end
end

function readmultimmcif(io::IO; gzip::Bool=false)
    if gzip
        gz = GzipDecompressorStream(io)
        tokens = tokenizecif(gz)
        close(gz)
    else
        tokens = tokenizecif(io)
    end
    idx = findall(t -> startswith(t, "data_"), tokens)
    multicif = Dict{String, MMCIFDict}()
    push!(idx, length(tokens) + 1) # Add end index for last block
    for i in 1:(length(idx) - 1)
        idx_start = idx[i]
        idx_end = idx[i + 1] - 1
        data_token = tokens[idx_start]
        cif = MMCIFDict()
        cif_name = data_token[6:end]
        if haskey(multicif, cif_name)
            throw(ErrorException("More than one mmCIF data block for data_$cif_name"))
        end
        cif["data_"] = [cif_name]
        populatedict!(cif, tokens[(idx_start + 1):idx_end])
        multicif[cif_name] = cif
    end
    return multicif
end

"""
    writemultimmcif(filepath, cifs; gzip=false)
    writemultimmcif(io, cifs; gzip=false)

Write multiple `MMCIFDict`s as a `Dict{String, MMCIFDict}` to a filepath or stream.
The keyword argument `gzip` (default `false`) determines if the output is gzipped.
"""
function writemultimmcif(filepath::AbstractString, cifs::Dict{String, MMCIFDict}; gzip::Bool=false)
    open(filepath, "w") do f
        writemultimmcif(f, cifs; gzip=gzip)
    end
end

function writemultimmcif(io::IO, cifs::Dict{String, MMCIFDict}; gzip::Bool=false)
    if gzip
        io = GzipCompressorStream(io)
    end
    for (k, c) in cifs
        c_key = c["data_"][1]
        if k != c_key
            @warn "writemultimmcif: MMCIFDict for key \"$k\" has different \"data_\" key (\"$c_key\")"
        end
        writemmcif(io, c)
        write(io, "##\n\n")
    end
    if gzip
        close(io)
    end
end


text:
using Cthulhu
using REPL
using InteractiveUtils
using Test
using Random
using StaticArrays

function firstassigned(specializations::Core.SimpleVector)
    # the methodinstance may not be first (annoying)
    for i = 1:length(specializations)
        if isassigned(specializations, i)
            return specializations[i]
        end
    end
    return nothing
end
function firstassigned(specializations)
    mis = []
    Base.visit(specializations) do mi
        isempty(mis) && push!(mis, mi)
    end
    return mis[1]
end

function process(@nospecialize(f), @nospecialize(TT); optimize=true)
    (interp, mi) = Cthulhu.mkinterp(f, TT)
    (ci, rt, infos, slottypes) = Cthulhu.lookup(interp, mi, optimize)
    Cthulhu.preprocess_ci!(ci, mi, optimize, Cthulhu.CthulhuConfig(dead_code_elimination=true))
    ci, infos, mi, rt, slottypes
end

function find_callsites_by_ftt(@nospecialize(f), @nospecialize(TT); optimize=true)
    ci, infos, mi, _, slottypes = process(f, TT; optimize=optimize)
    callsites = Cthulhu.find_callsites(ci, infos, mi, slottypes)
end

function test()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

let callsites = find_callsites_by_ftt(test, Tuple{})
    @test length(callsites) >= 4
end

let callsites = find_callsites_by_ftt(test, Tuple{}; optimize=false)
    @test length(callsites) == 4
end

@testset "Expr heads" begin
    # Check that the Expr head (:invoke or :call) is preserved
    @noinline twice(x::Real) = 2x
    calltwice(c) = twice(c[1])

    callsites = find_callsites_by_ftt(calltwice, Tuple{Vector{Float64}})
    @test length(callsites) == 1 && callsites[1].head === :invoke
    io = IOBuffer()
    print(io, callsites[1])
    @test occursin("invoke twice(::Float64)::Float64", String(take!(io)))

    callsites = find_callsites_by_ftt(calltwice, Tuple{Vector{AbstractFloat}})
    @test length(callsites) == 1 && callsites[1].head === :call
    io = IOBuffer()
    print(io, callsites[1])
    @test occursin("call twice(::AbstractFloat)", String(take!(io)))

    # Note the failure of `callinfo` to properly handle specialization
    @test_broken Cthulhu.callinfo(Tuple{typeof(twice), AbstractFloat}, AbstractFloat) isa Cthulhu.MultiCallInfo
end

# Check that we see callsites that are the rhs of assignments
@noinline bar_callsite_assign() = nothing
function foo_callsite_assign()
    x = bar_callsite_assign()
    x
end
let callsites = find_callsites_by_ftt(foo_callsite_assign, Tuple{}; optimize=false)
    @test length(callsites) == 1
end

@eval function call_rt()
    S = $(Core.Compiler.return_type)(+, Tuple{Int, Int})
end
let callsites = find_callsites_by_ftt(call_rt, Tuple{}; optimize=false)
    @test length(callsites) == 1
end

function call_by_apply(args...)
    identity(args...)
end
let callsites = find_callsites_by_ftt(call_by_apply, Tuple{Tuple{Int}}; optimize=false)
    @test length(callsites) == 1
end

if Base.JLOptions().check_bounds == 0
Base.@propagate_inbounds function f(x)
    @boundscheck error()
end
g(x) = @inbounds f(x)
h(x) = f(x)

let (CI, _, _, _, _) = process(g, Tuple{Vector{Float64}})
    @test length(CI.stmts) == 3
end

let (CI, _, _, _, _) = process(h, Tuple{Vector{Float64}})
    @test length(CI.stmts) == 2
end
end

# Something with many methods, but enough to be under the match limit
g_matches(a::Int, b::Int) = a+b
g_matches(a::Float64, b::Float64) = a+b
f_matches(a, b) = g_matches(a, b)
let callsites = find_callsites_by_ftt(f_matches, Tuple{Any, Any}; optimize=false)
    @test length(callsites) == 1
    callinfo = callsites[1].info
    @test callinfo isa Cthulhu.MultiCallInfo
end

# Failed return_type
only_ints(::Integer) = 1
return_type_failure(::T) where T = Base._return_type(only_ints, Tuple{T})
let callsites = find_callsites_by_ftt(return_type_failure, Tuple{Float64}, optimize=false)
    @test length(callsites) == 1
    callinfo = callsites[1].info
    @test callinfo isa Cthulhu.ReturnTypeCallInfo
    callinfo = callinfo.called_mi
    @test callinfo isa Cthulhu.MultiCallInfo
    @test length(callinfo.callinfos) == 0
end

# tasks
ftask() = @sync @async show(io, "Hello")
let callsites = find_callsites_by_ftt(ftask, Tuple{})
    task_callsites = filter(c->c.info isa Cthulhu.TaskCallInfo, callsites)
    @test !isempty(task_callsites)
    @test filter(c -> c.info.ci isa Cthulhu.FailedCallInfo, task_callsites) == []
end

##
# Union{f, g}
##
callf(f::F, x) where F = f(x)
let callsites = find_callsites_by_ftt(callf, Tuple{Union{typeof(sin), typeof(cos)}, Float64})
    @test !isempty(callsites)
    if length(callsites) == 1
        @test first(callsites).info isa Cthulhu.MultiCallInfo
        callinfos = first(callsites).info.callinfos
        @test !isempty(callinfos)
        mis = map(Cthulhu.get_mi, callinfos)
        @test any(mi->mi.def.name == :cos, mis)
        @test any(mi->mi.def.name == :sin, mis)
    else
        @test all(cs->cs.info isa Cthulhu.MICallInfo, callsites)
    end
end

function toggler(toggle)
    if toggle
        g = sin
    else
        g = cos
    end
    g(0.0)
end
let callsites = find_callsites_by_ftt(toggler, Tuple{Bool})
    @test !isempty(callsites)
    if length(callsites) == 1
        @test first(callsites).info isa Cthulhu.MultiCallInfo
        callinfos = first(callsites).info.callinfos
        @test !isempty(callinfos)
        mis = map(Cthulhu.get_mi, callinfos)
        @test any(mi->mi.def.name == :cos, mis)
        @test any(mi->mi.def.name == :sin, mis)
    else
        @test all(cs->cs.info isa Cthulhu.MICallInfo, callsites)
    end
end

@testset "Varargs" begin
    function fsplat(::Type{Int}, a...)
        z = zero(Int)
        for v in a
            z += v
        end
        return z
    end
    gsplat1(T::Type, a...) = fsplat(T, a...)   # does not force specialization
    hsplat1(A) = gsplat1(eltype(A), A...)
    io = IOBuffer()
    iolim = Cthulhu.TextWidthLimiter(io, 80)
    for (Atype, haslen) in ((Tuple{Tuple{Int,Int,Int}}, true),
                            (Tuple{Vector{Int}}, false),
                            (Tuple{SVector{3,Int}}, true))
        let callsites = find_callsites_by_ftt(hsplat1, Atype; optimize=false)
            @test !isempty(callsites)
            cs = callsites[end]
            @test cs isa Cthulhu.Callsite
            mi = cs.info.mi
            @test mi.specTypes.parameters[end] === (haslen ? Int : Vararg{Int})
            Cthulhu.show_callinfo(iolim, cs.info)
            @test occursin("…", String(take!(io))) != haslen
        end
    end
end

@testset "MaybeUndef" begin
    function undef(b::Bool)
        b || @goto final_step
        str = randstring(8)
        @label final_step
        return str*"end"
    end
    @test isa(undef(true), String)
    @test_throws UndefVarError undef(false)
    cs = find_callsites_by_ftt(undef, Tuple{Bool})[end]
    @test cs.head === :invoke
    @test cs.info.mi.def == which(string, (String,String))
end

like_cat(dims, xs::AbstractArray{T}...) where T = like_cat_t(T, xs...; dims=dims)
like_cat_t(::Type{T}, xs...; dims) where T = T

let callsites = find_callsites_by_ftt(like_cat, Tuple{Val{3}, Vararg{Matrix{Float32}}})
    @test length(callsites) == 1
    cs = callsites[1]
    @test cs isa Cthulhu.Callsite
    mi = cs.info.mi
    @test mi.specTypes.parameters[4] === Type{Float32}
end

@testset "warntype variables" begin
    src, rettype = code_typed(identity, (Any,); optimize=false)[1]
    io = IOBuffer()
    ioctx = IOContext(io, :color=>true)
    Cthulhu.cthulhu_warntype(ioctx, src, rettype, :none, true)
    str = String(take!(io))
    @test occursin("x\e[91m\e[1m::Any\e[22m\e[39m", str)
end

# issue #132
f132(w, dim) = [i == dim ? w[i]/2 : w[i]/1 for i in eachindex(w)]
interp, mi = Cthulhu.mkinterp(f132, (Vector{Int}, Int))
@test isa(mi, Core.MethodInstance)   # just check that the above succeeded

## Functions for "backedges & treelist"
# The printing changes when the functions are defined inside the testset
fbackedge1() = 1
fbackedge2(x) = x > 0 ? fbackedge1() : -fbackedge1()
fst1(x) = backtrace()
@inline fst2(x) = fst1(x)
@noinline fst3(x) = fst2(x)
@inline fst4(x) = fst3(x)
fst5(x) = fst4(x)

@testset "backedges and treelist" begin
    @test fbackedge2(0.2) == 1
    @test fbackedge2(-0.2) == -1
    mspec = @which(fbackedge1()).specializations
    mi = isa(mspec, Core.SimpleVector) ? mspec[1] : mspec.func
    root = Cthulhu.treelist(mi)
    @test Cthulhu.count_open_leaves(root) == 2
    @test root.data.callstr == "fbackedge1()"
    @test root.children[1].data.callstr == " fbackedge2(::Float64)"

    # issue #114
    unspecva(@nospecialize(i::Int...)) = 1
    @test unspecva(1, 2) == 1
    mi = firstassigned(first(methods(unspecva)).specializations)
    root = Cthulhu.treelist(mi)
    @test occursin("Vararg", root.data.callstr)

    # treelist for stacktraces
    tree = Cthulhu.treelist(fst5(1.0))
    @test match(r"fst1 at .*:\d+ => fst2 at .*:\d+ => fst3\(::Float64\) at .*:\d+", tree.data.callstr) !== nothing
    @test length(tree.children) == 1
    child = tree.children[1]
    @test match(r" fst4 at .*:\d+ => fst5\(::Float64\) at .*:\d+", child.data.callstr) !== nothing
end

@testset "ascend" begin
    # This tests only the non-interactive "look up the caller" portion
    callee(x) = 2x
    function caller(x)
        val = 0.0
        val += callee(3); line1 = @__LINE__
        val += callee(3.0); line2 = @__LINE__
        val += callee(x); line3 = @__LINE__
        val = sum([val])     # FIXME: without this line, `lookup` fails because codeinst.inferred === nothing
        return val, line1, line2, line3
    end
    _, line1, line2, line3 = caller(7)
    micaller = Cthulhu.get_specialization(caller, Tuple{Int})
    micallee_Int = Cthulhu.get_specialization(callee, Tuple{Int})
    micallee_Float4 = Cthulhu.get_specialization(callee, Tuple{Float64})
    info, lines = only(Cthulhu.find_caller_of(micallee_Int, micaller))
    @test info == (:caller, Symbol(@__FILE__), 0) && lines == [line1, line3]
    info, lines = only(Cthulhu.find_caller_of(micallee_Float4, micaller))
    @test info == (:caller, Symbol(@__FILE__), 0) && lines == [line2]

    # Detection in optimized (post-inlining) code
    @noinline nicallee(x) = 2x
    midcaller(x) = nicallee(x), @__LINE__
    function outercaller(x)
        val, line2 = midcaller(x); line1 = @__LINE__
        val = sum([val])     # FIXME: without this line, `lookup` fails because codeinst.inferred === nothing
        return val, line1, line2
    end
    _, line1, line2 = outercaller(7)
    micaller = Cthulhu.get_specialization(outercaller, Tuple{Int})
    micallee = Cthulhu.get_specialization(nicallee, Tuple{Int})
    callerinfo = Cthulhu.find_caller_of(micallee, micaller)
    @test length(callerinfo) == 2
    info, lines = callerinfo[1]
    @test info == (:outercaller, Symbol(@__FILE__), 0)
    @test lines == [line1]
    info, lines = callerinfo[2]
    @test info == (:midcaller, Symbol(@__FILE__), 1)
    @test lines == [line2]
end

##
# Cthulhu config test
##
let config = Cthulhu.CthulhuConfig(enable_highlighter=false)
    for lexer in ["llvm", "asm"]
        @test sprint() do io
            Cthulhu.highlight(io, "INPUT", lexer, config)
        end == "INPUT\n"
    end
end

let config = Cthulhu.CthulhuConfig(
    highlighter = `I_am_hoping_this_command_does_not_exist`,
    enable_highlighter = true,
)
    for lexer in ["llvm", "asm"]
        @test begin
            @test_logs (:warn, r"Highlighter command .* does not exist.") begin
                sprint() do io
                    Cthulhu.highlight(io, "INPUT", lexer, config)
                end == "INPUT\n"
            end
        end
    end
end

let config = Cthulhu.CthulhuConfig(
    # Implementing `cat` in Julia:
    highlighter = `$(Base.julia_cmd()) -e "write(stdout, read(stdin))" --`,
    enable_highlighter = true,
)
    for lexer in ["llvm", "asm"]
        @test sprint() do io
            Cthulhu.highlight(io, "INPUT", lexer, config)
        end == "INPUT\n"
        @test sprint() do io
            Cthulhu.highlight(io, "INPUT\n", lexer, config)
        end == "INPUT\n"
    end
end

###
### Printer test:
###
ci, infos, mi, rt, slottypes = process(test, Tuple{});
io = IOBuffer()
Cthulhu.cthulu_typed(io, :none, ci, rt, mi, true, false)
str = String(take!(io))
print(str)
# test by bounding the number of lines printed
@test_broken count(isequal('\n'), str) <= 50


text:
export Waveform, WaveformHeader

struct RawWaveformHeader
    version::UInt16 # Version of the header
    flags::UInt64 # bit field with flags
    measurement_uid::UInt32 # Unique id for the measurement
    scan_counter::UInt32 # Number of the acquisition after this waverform
    time_stamp::UInt32 # Time stamp of the waveform
    number_of_samples::UInt16 # Number of samples in the waveform
    channels::UInt16 # Number of channels per sample
    sample_time_us::Float32 # Time between samples in micro second
    waveform_id::UInt16 # Id matchin the types specified in the xml header
end

struct WaveformHeader
    version::UInt16 # Version of the header
    flags::UInt64 # bit field with flags
    measurement_uid::UInt32 # Unique id for the measurement
    scan_counter::UInt32 # Number of the acquisition after this waverform
    time_stamp::UInt32 # Time stamp of the waveform
    sample_time_us::Float32 # Time between samples in micro second
    waveform_id::UInt16 # Id matchin the types specified in the xml header
end

struct Waveform
    header::WaveformHeader
    data::Array{UInt32,2}
end
function RawWaveformHeader(wav::Waveform)
    names = fieldnames(WaveformHeader)
    fields = Dict((n,getfield(wav.header,n)) for n in names)

    fields[:number_of_samples] = size(wav.data)[1]
    fields[:channels] = size(wav.data)[2]
    RawWaveformHeader((fieldnames(RawWaveformHeader) .|> n -> fields[n])...)

end

@generated function WaveformHeader(raw::RawWaveformHeader)
    names = fieldnames(WaveformHeader)
    args = names .|> s -> Expr(:call,:getfield,:raw,QuoteNode(s))
    :(WaveformHeader($(args...)))
end

write(io::IO, header::RawWaveformHeader) = fieldnames(RawWaveformHeader) .|> getfield $ header .|> MRD.unsafe_write $ io
read(io::IO, ::Type{RawWaveformHeader}) = RawWaveformHeader((fieldnames(RawWaveformHeader) .|> fieldtype $ RawWaveformHeader .|> MRD.unsafe_read $ io)...)

function read(io::IO, ::Type{Waveform})
    header = MRD.read(io::IO, RawWaveformHeader)
    data = Array{UInt32,2}(undef, header.number_of_samples, header.channels)
    Base.read!(io, data)
    return Waveform(WaveformHeader(header), data)
end

function write(io::IO, wav::Waveform)
    write(io, RawWaveformHeader(wav))
    Base.write(io,wav.data)
end


text:
using Statistics
using LinearAlgebra
using Distributions

include(string(pwd(),"/src/gensys.jl"))
include(string(pwd(),"/hmc/diffs_estrut_v3.jl"))
include(string(pwd(),"/misc/matrix_no_kalman.jl"))

# See Scjmitt-Grohé paper on evaluating likelihoods without the Kalman filter to get the notation

function dSx(h,Sx,dh,dQ)
    m = size(h,1)
    Kmm = commutation_matrix(m,m)
    D = duplication_matrix(m)

    b1 = (I(m^2) - kron(h,h))
    b2 = kron(h*Sx,I(m)) + kron(I(m),h*Sx)*Kmm
    res = -inv(b1)*(b2*dh' - D*dQ')

    return res
end

function kron_and_sum(M,k)
    res = zeros(size(M).^2)
    for j in 1:k
        term = kron(M^(k-j),(M')^j)
        res += term
    end
    return res
end

function dSy(g,h,Sx,dSx,dh,j) #j is the power of h
    m = size(h,1)
    Kmm = commutation_matrix(m,m)
    h_pow = h^j
    ks = kron_and_sum(h,j)
    b1 = kron(g,g*h*Sx)*ks
    b2 = kron(g*h_pow*Sx,g)
    b3 = kron(g*h_pow,g*h)
    res = (b1*Kmm+b2)*dh' + b3*dSx
    return res
end

function diff_ll(P,y,dP)
    P_inv = inv(P)
    return -1/2*vec(P_inv)'*dP + 1/2*kron(y'*P_inv,y'*P_inv)*dP
end

function log_like_dsge(par,data;kalman_tol = 1e-10)
    #order to par
    #alfa
    #beta
    #epsilon
    #theta
    #sig
    #sigma: this is the std dev of the innovation
    #phi
    #phi_pi
    #phi_y
    #rho_v

    #data will have t x p dimension: lines are periods p are variables

    alfa = par[1]
    bet = par[2]
    epsilon = par[3]
    theta = par[4]
    sig = par[5]
    #par 6 is coded bellow see line 56
    phi = par[7]
    phi_pi = par[8]
    phi_y = par[9]
    rho_v = par[10]

    THETA = (1-alfa)/(1-alfa+alfa*epsilon)
    lamb = (1-theta)*(1-bet*theta)/theta*THETA
    kappa = lamb*(sig+(phi+alfa)/(1-alfa))

    nobs = size(data,1)
    l = size(data,2)

    GAMMA_0 = [bet  0    0  0;
               1    sig  0  0;
               0    0    0  0;
               0    0    0  1]

    GAMMA_1 = [ 1       -kappa  0   0;
                0        sig    1   0;
               -phi_pi  -phi_y  1  -1;
                0        0      0   rho_v]

    PSI = [0; 0; 0; 1]

    PI = [bet  0;
          1    sig;
          0    0;
          0    0]

    p = size(GAMMA_1,1) #number of endogenous vars

    sol = gensys(GAMMA_0,GAMMA_1,PSI,PI; verbose = false)

    if sum(sol.eu) != 2
        return -Inf, repeat([0],length(par))
    end

    #Sig = zeros(p,p)
    #Sig[4,4] = par[6]

    G = zeros(1,p)
    G[1,2] = 1

    A = sol.Theta1
    R = [0] .+ 1e-8
    Q = par[6]^2*sol.Theta2*sol.Theta2'

    dA,dB = diff_mod(par)
    #dB = dB[1:4,:]
    dA = dA'*I
    dG = zeros(10,4)
    dQ = (kron(sol.Theta2,I(p)) + kron(I(p),sol.Theta2)*commutation_matrix(4,1))*dB
    dQ = pinv(duplication_matrix(4))*dQ
    dQ = dQ'
    dR = zeros(size(par,1),1)

    Sx = solve_lyapunov_vec(A,Q)

    S = build_variance(G,A,Q,nobs)
    S = copy(S)

    dSx_mat = dSx(A,Q,dA,dQ)
    dSy_foo(j) = dSy(G,A,Sx,dSx_mat,dA,j)

    d_vecS = map(dSy_foo,1:(nobs))

    sel_mat = Toeplitz(1:nobs,1:nobs)
    sel_mat = Int.(sel_mat)

    dS_mat = d_vecS[sel_mat]
    dS_mat = vec(dS_mat)
    dS_mat = mapreduce(x->dS_mat[x],vcat,1:nobs^2)

    dll = diff_ll(S,data,dS_mat)

    dist = MvNormal(S)

    llh = logpdf(dist,data)

    return llh,dll
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using TranscodingStreams
using Test
using Pkg

# Tool tests
# ----------

using TranscodingStreams:
    Buffer, Memory,
    bufferptr, buffersize, buffermem,
    marginptr, marginsize, marginmem,
    readbyte!, writebyte!,
    #=ismarked,=# mark!, unmark!, reset!,
    makemargin!, emptybuffer!

@testset "Buffer" begin
    buf = Buffer(1024)
    @test buf isa Buffer
    @test length(buf.data) == 1024

    data = Vector{UInt8}(b"foobar")
    buf = Buffer(data)
    @test buf isa Buffer
    @test bufferptr(buf) === pointer(data)
    @test buffersize(buf) === 6
    @test buffermem(buf) === Memory(pointer(data), 6)
    @test marginptr(buf) === pointer(data) + 6
    @test marginsize(buf) === 0
    @test marginmem(buf) === Memory(pointer(data)+6, 0)

    buf = Buffer(2)
    writebyte!(buf, 0x34)
    writebyte!(buf, 0x9f)
    @test buffersize(buf) === 2
    @test readbyte!(buf) === 0x34
    @test readbyte!(buf) === 0x9f
    @test buffersize(buf) === 0

    buf = Buffer(16)
    @test !TranscodingStreams.ismarked(buf)
    @test mark!(buf) == 1
    @test TranscodingStreams.ismarked(buf)
    @test unmark!(buf)
    @test !TranscodingStreams.ismarked(buf)
    @test !unmark!(buf)

    buf = Buffer(16)
    mark!(buf)
    writebyte!(buf, 0x34)
    writebyte!(buf, 0x99)
    reset!(buf)
    @test !TranscodingStreams.ismarked(buf)
    @test readbyte!(buf) === 0x34
    @test readbyte!(buf) === 0x99

    buf = Buffer(16)
    @test makemargin!(buf, 0) === 16
    writebyte!(buf, 0x34)
    @test makemargin!(buf, 0) === 15
    writebyte!(buf, 0x99)
    @test makemargin!(buf, 20) === 20
    emptybuffer!(buf)
    @test makemargin!(buf, 0) === 22
end

@testset "Memory" begin
    data = Vector{UInt8}(b"foobar")
    mem = TranscodingStreams.Memory(pointer(data), sizeof(data))
    @test mem isa TranscodingStreams.Memory
    @test mem.ptr === pointer(data)
    @test mem.size === length(mem) === UInt(sizeof(data))
    @test lastindex(mem) === 6
    @test mem[1] === UInt8('f')
    @test mem[2] === UInt8('o')
    @test mem[3] === UInt8('o')
    @test mem[4] === UInt8('b')
    @test mem[5] === UInt8('a')
    @test mem[6] === UInt8('r')
    @test_throws BoundsError mem[7]
    @test_throws BoundsError mem[0]
    mem[1] = UInt8('z')
    @test mem[1] === UInt8('z')
    mem[3] = UInt8('!')
    @test mem[3] === UInt8('!')
    @test_throws BoundsError mem[7] = 0x00
    @test_throws BoundsError mem[0] = 0x00

    data = Vector{UInt8}(b"foobar")
    mem = TranscodingStreams.Memory(data)
    @test mem isa TranscodingStreams.Memory
    @test mem.ptr == pointer(data)
    @test mem.size == sizeof(data)
end

@testset "Stats" begin
    stats = TranscodingStreams.Stats(1,2,3,4)
    @test repr(stats) ==
    """
    TranscodingStreams.Stats:
      in: 1
      out: 2
      transcoded_in: 3
      transcoded_out: 4"""
    @test stats.in == 1
    @test stats.out == 2
    @test stats.transcoded_in == 3
    @test stats.transcoded_out == 4
end

@testset "Utils" begin
    @test TranscodingStreams.splitkwargs(
        [(:foo, 1), (:bar, true), (:baz, :ok)], (:foo,)) ==
        ([(:foo, 1)], [(:bar, true), (:baz, :ok)])
end

# Codec tests
# -----------

include("codecnoop.jl")
include("codecinvalid.jl")
include("codecquadruple.jl")

# Test third-party codec packages.
for pkg in ["CodecZlib", "CodecXz", "CodecZstd", "CodecBase"]
    Pkg.test(pkg)
end


text:
using BenchmarkTools
using Statistics
using UltraDark
using UltraDark: take_steps!

nthreads = Threads.nthreads()

resol = try
    parse(Int64, ARGS[1])
catch
    @warn "Resolution not supplied, using default"
    64
end

Δt = 0.1
n_steps = 10

grids = Grids(1.0, resol)

# b = @benchmark take_steps!($grids, $1.0, $Δt, 10, $(OutputConfig("output", 1:2)), $(Config.constant_scale_factor))
# time_mean = mean(b.times / 1e9 / n_steps)

# Run once to ensure function is precompiled
take_steps!(grids, 1.0, Δt, n_steps, OutputConfig("output", 1:2), Config.constant_scale_factor)

# Collect data
res = @timed take_steps!(grids, 1.0, Δt, n_steps, OutputConfig("output", 1:2), Config.constant_scale_factor)

time_mean = mean(res[2] / n_steps)

println("$nthreads, $resol, $time_mean")


text:
push!(LOAD_PATH,"../src")
using Documenter: Documenter, makedocs, deploydocs, doctest, DocMeta
using LabeledArrays
using Test
using Literate

docdir = joinpath(dirname(pathof(LabeledArrays)),"../docs/src/")
mandir = joinpath(docdir,"man")

DocMeta.setdocmeta!(LabeledArrays, :DocTestSetup, quote
    using LabeledArrays
end; recursive=true)

doctest(LabeledArrays; fix=true)

makedocs(;
         source = docdir,
         modules=[LabeledArrays],
         authors="Gregor Kappler",
         # repo="https://github.com/gkappler/CombinedParsers.jl/blob/{commit}{path}#L{line}",
         sitename="LabeledArrays.jl",
         # format=Documenter.HTML(;
         #                        prettyurls=get(ENV, "CI", "false") == "true",
         #                        canonical="https://gkappler.github.io/CombinedParsers.jl",
         #                        assets=String[],
         #                        ),
         pages=[
             "Home" => "index.md",
             # "Library" => Any[
             #     "Public" => "lib/public.md",
             #     ##"Internals" => "lib/internals.md"
             # ],
             # "Developer Guide" => "developer.md"
         ]
         )

deploydocs(;
    repo="github.com/gkappler/LabeledArrays.jl",
)


text:
#
# The main file focuses on defining types and the dimensionless form
# of the equations and parameters.
#
# The functions to convert the model results to physical units are all
# given in the physical.jl file, while utility functions to calculate
# stellar parameters for a spectrum (in particular black body ones) are
# in stars.jl
#

module DustyHII

using DifferentialEquations
using QuadGK
using NumericalIntegration
using Unitful, UnitfulRecipes, UnitfulAstro
using Dierckx

abstract type AbstractDustyHII end
abstract type AbstractModelParameters end
abstract type AbstractPhysicalModel end

"A struct that contains the parameters (γ, β) of the Draine 2011 model."
struct Draine2011Parameters <: AbstractModelParameters
γ::Real
β::Real
end

"""
Draine2011(solution, parameterss, u0, ϕ0, τ0, ymax, ymin)

A type containing the results of solving the differential equations for the Draine 2011 model

The solution is stored in `solution`, which is a `ODESolution`.
The parameters in `parameters`, which is a `Draine2011Parameters`
The initial conditions in `u0`, `ϕ0`, `τ0`
The range of validity in `ymax` and `ymin`.
"""
struct Draine2011 <: AbstractDustyHII
solution::ODESolution
parameters::Draine2011Parameters
u0::Real
ϕ0::Real
τ0::Real
ymin::Real
ymax::Real
end

"""
Draine2011Physical(model, T, Q0, hnui, σ)

A type containing the necessary information to give a physical realisation
of a model. It encapsulates the abstract unitless model solution as well
as the physical characteristics needed to create a model.

TODO: Support self-gravity
TODO (maybe): allow for type-checking for the extra parameters.
"""
mutable struct Draine2011Physical <: AbstractPhysicalModel
model::Draine2011
T
Q0
hnui
σ
end

# I would like a pretty printing method as well.
#
# This is far from being ready, it is here mostly as a placeholder.
function Base.show(io::IO, m::Draine2011)
    println(io, "Dusty HII region model")
    ymin=m.ymin ; ymax=m.ymax
    println(io, "    Defined for $ymin < y < $ymax")
end

# The actual equations are defined here
include("equations.jl")
# Physical conversion routines are kept in this file.
include("physical.jl")
# Stellar quantities are calculated here
include("stars.jl")

end # module


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/nvector/nvector_serial.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function N_VNew_Serial(vec_length::sunindextype)
    ccall((:N_VNew_Serial, libsundials_nvecserial), N_Vector, (sunindextype,), vec_length)
end

function N_VNewEmpty_Serial(vec_length::sunindextype)
    ccall((:N_VNewEmpty_Serial, libsundials_nvecserial), N_Vector, (sunindextype,), vec_length)
end

function N_VMake_Serial(vec_length::sunindextype, v_data)
    ccall((:N_VMake_Serial, libsundials_nvecserial), N_Vector, (sunindextype, Ptr{realtype}), vec_length, v_data)
end

function __N_VCloneVectorArray_Serial(count::Cint, w::N_Vector)
    ccall((:N_VCloneVectorArray_Serial, libsundials_nvecserial), Ptr{N_Vector}, (Cint, N_Vector), count, w)
end

function N_VCloneVectorArray_Serial(count, w)
    __w = convert(NVector, w)
    __N_VCloneVectorArray_Serial(convert(Cint, count), convert(N_Vector, __w))
end

function __N_VCloneVectorArrayEmpty_Serial(count::Cint, w::N_Vector)
    ccall((:N_VCloneVectorArrayEmpty_Serial, libsundials_nvecserial), Ptr{N_Vector}, (Cint, N_Vector), count, w)
end

function N_VCloneVectorArrayEmpty_Serial(count, w)
    __w = convert(NVector, w)
    __N_VCloneVectorArrayEmpty_Serial(convert(Cint, count), convert(N_Vector, __w))
end

function __N_VDestroyVectorArray_Serial(vs, count::Cint)
    ccall((:N_VDestroyVectorArray_Serial, libsundials_nvecserial), Void, (Ptr{N_Vector}, Cint), vs, count)
end

function N_VDestroyVectorArray_Serial(vs, count)
    __N_VDestroyVectorArray_Serial(vs, convert(Cint, count))
end

function __N_VGetLength_Serial(v::N_Vector)
    ccall((:N_VGetLength_Serial, libsundials_nvecserial), sunindextype, (N_Vector,), v)
end

function N_VGetLength_Serial(v)
    __v = convert(NVector, v)
    __N_VGetLength_Serial(convert(N_Vector, __v))
end

function __N_VPrint_Serial(v::N_Vector)
    ccall((:N_VPrint_Serial, libsundials_nvecserial), Void, (N_Vector,), v)
end

function N_VPrint_Serial(v)
    __v = convert(NVector, v)
    __N_VPrint_Serial(convert(N_Vector, __v))
end

function __N_VPrintFile_Serial(v::N_Vector, outfile)
    ccall((:N_VPrintFile_Serial, libsundials_nvecserial), Void, (N_Vector, Ptr{Cint}), v, outfile)
end

function N_VPrintFile_Serial(v, outfile)
    __v = convert(NVector, v)
    __N_VPrintFile_Serial(convert(N_Vector, __v), outfile)
end

function __N_VGetVectorID_Serial(v::N_Vector)
    ccall((:N_VGetVectorID_Serial, libsundials_nvecserial), N_Vector_ID, (N_Vector,), v)
end

function N_VGetVectorID_Serial(v)
    __v = convert(NVector, v)
    __N_VGetVectorID_Serial(convert(N_Vector, __v))
end

function __N_VCloneEmpty_Serial(w::N_Vector)
    ccall((:N_VCloneEmpty_Serial, libsundials_nvecserial), N_Vector, (N_Vector,), w)
end

function N_VCloneEmpty_Serial(w)
    __w = convert(NVector, w)
    __N_VCloneEmpty_Serial(convert(N_Vector, __w))
end

function __N_VClone_Serial(w::N_Vector)
    ccall((:N_VClone_Serial, libsundials_nvecserial), N_Vector, (N_Vector,), w)
end

function N_VClone_Serial(w)
    __w = convert(NVector, w)
    __N_VClone_Serial(convert(N_Vector, __w))
end

function __N_VDestroy_Serial(v::N_Vector)
    ccall((:N_VDestroy_Serial, libsundials_nvecserial), Void, (N_Vector,), v)
end

function N_VDestroy_Serial(v)
    __v = convert(NVector, v)
    __N_VDestroy_Serial(convert(N_Vector, __v))
end

function __N_VSpace_Serial(v::N_Vector, lrw, liw)
    ccall((:N_VSpace_Serial, libsundials_nvecserial), Void, (N_Vector, Ptr{sunindextype}, Ptr{sunindextype}), v, lrw, liw)
end

function N_VSpace_Serial(v, lrw, liw)
    __v = convert(NVector, v)
    __N_VSpace_Serial(convert(N_Vector, __v), lrw, liw)
end

function __N_VGetArrayPointer_Serial(v::N_Vector)
    ccall((:N_VGetArrayPointer_Serial, libsundials_nvecserial), Ptr{realtype}, (N_Vector,), v)
end

function N_VGetArrayPointer_Serial(v)
    __v = convert(NVector, v)
    __N_VGetArrayPointer_Serial(convert(N_Vector, __v))
end

function __N_VSetArrayPointer_Serial(v_data, v::N_Vector)
    ccall((:N_VSetArrayPointer_Serial, libsundials_nvecserial), Void, (Ptr{realtype}, N_Vector), v_data, v)
end

function N_VSetArrayPointer_Serial(v_data, v)
    __v = convert(NVector, v)
    __N_VSetArrayPointer_Serial(v_data, convert(N_Vector, __v))
end

function __N_VLinearSum_Serial(a::realtype, x::N_Vector, b::realtype, y::N_Vector, z::N_Vector)
    ccall((:N_VLinearSum_Serial, libsundials_nvecserial), Void, (realtype, N_Vector, realtype, N_Vector, N_Vector), a, x, b, y, z)
end

function N_VLinearSum_Serial(a, x, b, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    __N_VLinearSum_Serial(a, convert(N_Vector, __x), b, convert(N_Vector, __y), convert(N_Vector, __z))
end

function __N_VConst_Serial(c::realtype, z::N_Vector)
    ccall((:N_VConst_Serial, libsundials_nvecserial), Void, (realtype, N_Vector), c, z)
end

function N_VConst_Serial(c, z)
    __z = convert(NVector, z)
    __N_VConst_Serial(c, convert(N_Vector, __z))
end

function __N_VProd_Serial(x::N_Vector, y::N_Vector, z::N_Vector)
    ccall((:N_VProd_Serial, libsundials_nvecserial), Void, (N_Vector, N_Vector, N_Vector), x, y, z)
end

function N_VProd_Serial(x, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    __N_VProd_Serial(convert(N_Vector, __x), convert(N_Vector, __y), convert(N_Vector, __z))
end

function __N_VDiv_Serial(x::N_Vector, y::N_Vector, z::N_Vector)
    ccall((:N_VDiv_Serial, libsundials_nvecserial), Void, (N_Vector, N_Vector, N_Vector), x, y, z)
end

function N_VDiv_Serial(x, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    __N_VDiv_Serial(convert(N_Vector, __x), convert(N_Vector, __y), convert(N_Vector, __z))
end

function __N_VScale_Serial(c::realtype, x::N_Vector, z::N_Vector)
    ccall((:N_VScale_Serial, libsundials_nvecserial), Void, (realtype, N_Vector, N_Vector), c, x, z)
end

function N_VScale_Serial(c, x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VScale_Serial(c, convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VAbs_Serial(x::N_Vector, z::N_Vector)
    ccall((:N_VAbs_Serial, libsundials_nvecserial), Void, (N_Vector, N_Vector), x, z)
end

function N_VAbs_Serial(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VAbs_Serial(convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VInv_Serial(x::N_Vector, z::N_Vector)
    ccall((:N_VInv_Serial, libsundials_nvecserial), Void, (N_Vector, N_Vector), x, z)
end

function N_VInv_Serial(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VInv_Serial(convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VAddConst_Serial(x::N_Vector, b::realtype, z::N_Vector)
    ccall((:N_VAddConst_Serial, libsundials_nvecserial), Void, (N_Vector, realtype, N_Vector), x, b, z)
end

function N_VAddConst_Serial(x, b, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VAddConst_Serial(convert(N_Vector, __x), b, convert(N_Vector, __z))
end

function __N_VDotProd_Serial(x::N_Vector, y::N_Vector)
    ccall((:N_VDotProd_Serial, libsundials_nvecserial), realtype, (N_Vector, N_Vector), x, y)
end

function N_VDotProd_Serial(x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __N_VDotProd_Serial(convert(N_Vector, __x), convert(N_Vector, __y))
end

function __N_VMaxNorm_Serial(x::N_Vector)
    ccall((:N_VMaxNorm_Serial, libsundials_nvecserial), realtype, (N_Vector,), x)
end

function N_VMaxNorm_Serial(x)
    __x = convert(NVector, x)
    __N_VMaxNorm_Serial(convert(N_Vector, __x))
end

function __N_VWrmsNorm_Serial(x::N_Vector, w::N_Vector)
    ccall((:N_VWrmsNorm_Serial, libsundials_nvecserial), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWrmsNorm_Serial(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __N_VWrmsNorm_Serial(convert(N_Vector, __x), convert(N_Vector, __w))
end

function __N_VWrmsNormMask_Serial(x::N_Vector, w::N_Vector, id::N_Vector)
    ccall((:N_VWrmsNormMask_Serial, libsundials_nvecserial), realtype, (N_Vector, N_Vector, N_Vector), x, w, id)
end

function N_VWrmsNormMask_Serial(x, w, id)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __id = convert(NVector, id)
    __N_VWrmsNormMask_Serial(convert(N_Vector, __x), convert(N_Vector, __w), convert(N_Vector, __id))
end

function __N_VMin_Serial(x::N_Vector)
    ccall((:N_VMin_Serial, libsundials_nvecserial), realtype, (N_Vector,), x)
end

function N_VMin_Serial(x)
    __x = convert(NVector, x)
    __N_VMin_Serial(convert(N_Vector, __x))
end

function __N_VWL2Norm_Serial(x::N_Vector, w::N_Vector)
    ccall((:N_VWL2Norm_Serial, libsundials_nvecserial), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWL2Norm_Serial(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __N_VWL2Norm_Serial(convert(N_Vector, __x), convert(N_Vector, __w))
end

function __N_VL1Norm_Serial(x::N_Vector)
    ccall((:N_VL1Norm_Serial, libsundials_nvecserial), realtype, (N_Vector,), x)
end

function N_VL1Norm_Serial(x)
    __x = convert(NVector, x)
    __N_VL1Norm_Serial(convert(N_Vector, __x))
end

function __N_VCompare_Serial(c::realtype, x::N_Vector, z::N_Vector)
    ccall((:N_VCompare_Serial, libsundials_nvecserial), Void, (realtype, N_Vector, N_Vector), c, x, z)
end

function N_VCompare_Serial(c, x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VCompare_Serial(c, convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VInvTest_Serial(x::N_Vector, z::N_Vector)
    ccall((:N_VInvTest_Serial, libsundials_nvecserial), Cint, (N_Vector, N_Vector), x, z)
end

function N_VInvTest_Serial(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VInvTest_Serial(convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VConstrMask_Serial(c::N_Vector, x::N_Vector, m::N_Vector)
    ccall((:N_VConstrMask_Serial, libsundials_nvecserial), Cint, (N_Vector, N_Vector, N_Vector), c, x, m)
end

function N_VConstrMask_Serial(c, x, m)
    __c = convert(NVector, c)
    __x = convert(NVector, x)
    __m = convert(NVector, m)
    __N_VConstrMask_Serial(convert(N_Vector, __c), convert(N_Vector, __x), convert(N_Vector, __m))
end

function __N_VMinQuotient_Serial(num::N_Vector, denom::N_Vector)
    ccall((:N_VMinQuotient_Serial, libsundials_nvecserial), realtype, (N_Vector, N_Vector), num, denom)
end

function N_VMinQuotient_Serial(num, denom)
    __num = convert(NVector, num)
    __denom = convert(NVector, denom)
    __N_VMinQuotient_Serial(convert(N_Vector, __num), convert(N_Vector, __denom))
end


text:
function ly_generate(pkgname::String; kwargs...)
    preset = deepcopy(LY_PKGTEMPLATE)
    push!(preset.authors, "The Contributors of $pkgname")
    t = Template(; preset..., kwargs...)
    generate(pkgname, t)
end

function ly_register(
    repo_name::AbstractString,
    commit::AbstractString,
    push::Bool = true,
    branch::AbstractString = "master",
    kwargs...,
)
    try
        Base.SHA1(commit)
    catch e
        throw(ArgumentError("\"$commit\" is not a valid SHA-1"))
    end

    url = "https://github.com/Lyceum/$(repo_name).git"

    with_tempdir() do _
        git = create_git_cmd(LY_GITCONFIG)

        run(`$git clone $url $(pwd())`)
        run(`$git checkout master`)
        if !success(`$git merge-base --is-ancestor $commit HEAD`)
            error("$commit not found in master branch")
        end
        run(`$git checkout $commit`)

        pkg = read_project(projectfile_path(pwd()))

        registry_deps = map(reg -> reg.url, Pkg.Types.collect_registries())
        LY_REGISTRY.url in registry_deps || error("Please add LyceumRegistry: $LY_REGISTRY")

        rbrn = register(
            url,
            pkg,
            string(treehash(pwd()));
            registry = LY_REGISTRY.url,
            registry_deps = registry_deps,
            gitconfig = LY_GITCONFIG,
            push = push,
            branch = branch,
            kwargs...,
        )

        haskey(rbrn.metadata, "error") && error(rbrn.metadata["error"])

        return rbrn
    end
end

function incversion!(pkg::Union{Module,AbstractString}, which::Symbol; prerelease = (), build = ())
    pkg = pkg isa Module ? pkgdir(pkg) : pkg
    isdirty(pkg) && error("$pkg is dirty")

    projectfile = projectfile_path(pkg)
    project = read_project(projectfile)
    project.version = incversion(project.version, which, prerelease = prerelease)
    write_project(project, projectfile)

    git = create_git_cmd(LY_GITCONFIG, path = pkg)
    message = "New version: v$(project.version)"
    run(`$git add Project.toml`)
    run(`$git commit -qm $message`)

    return nothing
end

function incversion(version::Union{String,VersionNumber}, which::Symbol; prerelease= (), build = ())
    v = version isa String ? VersionNumber(version) : version
    major = v.major
    minor = v.minor
    patch = v.patch

    if which === :major
        major += 1
        minor = 0
        patch = 0
    elseif which === :minor
        minor += 1
        patch = 0
    elseif which === :patch
        patch += 1
    else
        error("which must be :major, :minor, or :patch. Got $which")
    end

    v.prerelease != () && @warn "prerelease not empty: $pre. Ignoring."
    prerelease = prerelease === :dev ? "DEV" : prerelease
    prerelease = prerelease isa Tuple ? prerelease : (prerelease, )

    v.build != () && @warn "build not empty: $build. Ignoring."
    build = build isa Tuple ? build : (build, )

    vnew = VersionNumber(major, minor, patch, prerelease, build)
    @info "Old version: $v. New version: $vnew"

    return vnew
end