text:
import Mercury as Hg
using ZMQ
using Sockets
using Test
using BenchmarkTools
if !isdefined(@__MODULE__, :TestMsg)
    include("jlout/test_msg_pb.jl")
end
ENV["JULIA_DEBUG"] = "Mercury"

@testset "ZmqSubscriber" begin
    Hg.reset_sub_count()
    ctx = ZMQ.Context()
    addr = ip"127.0.0.1"
    port = 5555

    @testset "Construction" begin
        sub = Hg.ZmqSubscriber(ctx, addr, port)

        @test isopen(sub)
        @test sub.port == port
        @test sub.ipaddr == addr
        @test sub.name == "subscriber_1"
        close(sub)
        @test !isopen(sub.socket)

        sub = Hg.ZmqSubscriber(ctx, addr, string(port))
        @test sub.port == port
        @test sub.ipaddr == addr
        @test sub.name == "subscriber_2"
        close(sub)

        sub = Hg.ZmqSubscriber(ctx, string(addr), port)
        @test sub.port == port
        @test sub.ipaddr == addr
        @test sub.name == "subscriber_3"
        close(sub)

        # Create 2 subscribers
        @test_nowarn begin
            sub1 = Hg.ZmqSubscriber(ctx, string(addr), port)
            sub2 = Hg.ZmqSubscriber(ctx, string(addr), port)
            close(sub1)
            close(sub2)
        end
    end

    @testset "Simple Pub/Sub" begin
        # Test simple pub/sub
        sub = Hg.ZmqSubscriber(ctx, addr, port)
        @test isopen(sub)
        msg = TestMsg(x = 10, y = 11, z = 12)

        pub = Hg.ZmqPublisher(ctx, addr, port)
        msg_out = TestMsg(x = 1, y = 2, z = 3)
        @test msg.x == 10
        @test msg.y == 11
        @test msg.z == 12
        sleep(0.2)  # needed to give time to set up publisher?

        # publish and receive a message
        @test !sub.flags.hasreceived
        i = 1
        for i = 1:100
            Hg.publish(pub, msg_out)
            if Hg.receive(sub, msg)
                break
            end
        end
        @test i < 100
        @test sub.flags.hasreceived

        # make sure the correct message was received
        @test msg.x == 1
        @test msg.y == 2
        @test msg.z == 3
        close(pub)
        close(sub)
    end

    @testset "Exchanging Byte Vectors" begin
        # Test simple pub/sub
        sub = Hg.ZmqSubscriber(ctx, addr, port)
        @test isopen(sub)
        pub = Hg.ZmqPublisher(ctx, addr, port)
        msg_out = Vector{UInt8}("Hello, world!")
        msg_in = zeros(UInt8, length(msg_out))
        sleep(0.2)  # needed to give time to set up publisher?

        # publish and receive a message
        @test !sub.flags.hasreceived
        i = 1
        for i = 1:100
            Hg.publish(pub, msg_out)
            if Hg.receive(sub, msg_in)
                break
            end
        end
        @test i < 100
        @test sub.flags.hasreceived

        # make sure the correct message was received
        @test all(msg_in .== msg_out)
        close(pub)
        close(sub)
    end

    @testset "Closing" begin
        sub = Hg.ZmqSubscriber(ctx, addr, port, name = "TestSub")
        pub = Hg.ZmqPublisher(ctx, addr, port, name = "TestPub")
        msg_in = TestMsg(x = 10, y = 11, z = 12)
        msg_out = TestMsg(x = 1, y = 2, z = 3)

        # Close the task by waiting for a receive
        cnt = 0
        timeout = 5.0 # seconds
        @test Hg.publish_until_receive(pub, sub, msg_out, msg_in, timeout)
        @test msg_in.x == msg_out.x

        # Should be able to close at any time
        close(sub)
        sleep(0.1)  # wait a little bit to allow for the lock to be acquired
        @test !isopen(sub)
        @test isopen(pub)
        close(pub)
    end

    @testset "Subscribe performance" begin
        do_publish = Threads.Atomic{Bool}(true)

        function pub_message(pub)
            msg_out = TestMsg(x = 1, y = 2, z = 3)
            do_publish
            i = 0
            while (do_publish[])
                msg_out.x = i
                Hg.publish(pub, msg_out)
                i += 1
                sleep(0.001)
            end
        end

        sub = Hg.ZmqSubscriber(ctx, addr, port)
        msg = TestMsg(x = 0, y = 0, z = 0)

        pub = Hg.ZmqPublisher(ctx, addr, port, name = "TestPub")
        msg_out = TestMsg(x = 10, y = 11, z = 12)
        pub_task = @async pub_message(pub)
        @test !istaskdone(pub_task)
        @test !istaskfailed(pub_task)
        # do_publish[] = false
        @test isopen(pub)

        @test isopen(sub)

        # Retrieve 2 messages to make sure the publisher is working
        Hg.receive(sub, msg)
        x_prev = msg.x
        sleep(1.0)
        Hg.receive(sub, msg)
        x_new = msg.x
        @test x_new - x_prev > 50

        # Benchmark the receive
        b = @benchmark Hg.receive($sub, $msg)
        @test maximum(b.gctimes) == 0  # no garbage collection
        @test b.memory == 0            # no dynamic memory allocations

        # Close task and sockets
        do_publish[] = false
        wait(pub_task)
        close(sub)
        close(pub)
    end

    @testset "Testing Subscriber Conflate" begin
        ## Test receive performance
        function pub_message(pub)
            rate = 100  # Publishing at 100 Hz
            lrl = Hg.LoopRateLimiter(rate)

            msg_out = TestMsg(x = 1, y = 2, z = 3)
            global do_publish
            i = 0
            Hg.@rate while (do_publish)
                msg_out.x = i
                Hg.publish(pub, msg_out)
                i += 1
                sleep(0.001)
            end lrl
        end

        sub = Hg.ZmqSubscriber(ctx, addr, port, name = "TestSub")
        pub = Hg.ZmqPublisher(ctx, addr, port, name = "TestPub")
        msg = TestMsg(x = 10, y = 11, z = 12)

        # Publish message in a separate task (really fast)
        global do_publish = true
        pub_task = @task pub_message(pub)
        schedule(pub_task)
        @test !istaskdone(pub_task)

        Hg.receive(sub, msg)
        first_rec = msg.x
        sleep(0.5)
        Hg.receive(sub, msg)
        second_rec = msg.x
        @test second_rec > first_rec + 5

        do_publish = false
        sleep(0.1)
        @test istaskdone(pub_task)
        close(pub)
        close(sub)
    end

end

@testset "Published/Subscribed ZMQ Messages" begin
    Hg.reset_sub_count()
    ctx = ZMQ.Context()
    addr = ip"127.0.0.1"
    port = 5555

    sub = Hg.ZmqSubscriber(ctx, addr, port, name = "TestSub")
    pub = Hg.ZmqPublisher(ctx, addr, port, name = "TestPub")
    msg = TestMsg(x = 10, y = 11, z = 12)
    msg_out = TestMsg(x = 1, y = 2, z = 3)

    submsg = Hg.SubscribedMessage(msg, sub)
    pubmsg = Hg.PublishedMessage(msg_out, pub)

    @test !Hg.getflags(submsg.sub).hasreceived

    has_heard = false
    while (!has_heard)
        sleep(0.001)
        has_heard = Hg.receive(submsg)
        Hg.publish(pubmsg)
    end
    @test has_heard

    got_new = false
    Hg.on_new(submsg) do msg
        @test msg isa Hg.MercuryMessage
        got_new = true
    end
    @test got_new

    @test isopen(sub)
    @test isopen(pub)
    Hg.forceclose(sub)
    close(pub)
    sleep(0.2)  # give a little bit of time to close the task
end


text:
abstract type Objective end

mutable struct TrackingObjective{Q,U} <: Objective
    q::Vector{Q}
    u::Vector{U}
end

function TrackingObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.dim.q})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.dim.u})) for t = 1:H])
    return TrackingObjective(q, u)
end

function gradient!(res, obj::TrackingObjective, core, traj, ref_traj)
    for t = 1:traj.H
        # Cost function
        delta!(core.Δq[t], traj.q[t+2], ref_traj.q[t+2])
        delta!(core.Δu[t], traj.u[t], ref_traj.u[t])

        res.q2[t] .+= obj.q[t] * core.Δq[t]
        res.u1[t] .+= obj.u[t] * core.Δu[t]
    end
end

function hessian!(hess, obj::TrackingObjective)
    for t = 1:length(obj.u)
        # Cost function
        hess.obj_q2[t] .+= obj.q[t]
        hess.obj_u1[t] .+= obj.u[t]
    end
end

mutable struct TrackingVelocityObjective{Q,V,U} <: Objective
    q::Vector{Q}
    v::Vector{V}
    u::Vector{U}
end

function TrackingVelocityObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.dim.q})) for t = 1:H],
    v = [Diagonal(zeros(SizedVector{model.dim.q})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.dim.u})) for t = 1:H])
    return TrackingVelocityObjective(q, v, u)
end

function gradient!(res, obj::TrackingVelocityObjective, core, traj, ref_traj)
    for t = 1:traj.H
        # Cost function
        delta!(core.Δq[t], traj.q[t+2], ref_traj.q[t+2])
        delta!(core.Δu[t], traj.u[t], ref_traj.u[t])

        res.q2[t] .+= obj.q[t] * core.Δq[t]
        res.u1[t] .+= obj.u[t] * core.Δu[t]

        # velocity
        res.q2[t] .+= obj.v[t] * (traj.q[t+2] - traj.q[t+1])
        t == 1 && continue
        res.q2[t-1] .-= obj.v[t] * (traj.q[t+2] - traj.q[t+1])
    end
end

function hessian!(hess, obj::TrackingVelocityObjective)
    for t = 1:length(obj.u)
        # Cost function
        hess.obj_q2[t] .+= obj.q[t]
        hess.obj_u1[t] .+= obj.u[t]

        # velocity
        hess.obj_q2[t] .+= obj.v[t]
        t == 1 && continue
        hess.obj_q2[t-1] .+= obj.v[t]
        hess.obj_q1q2[t-1] .-= obj.v[t]
        hess.obj_q2q1[t-1] .-= obj.v[t]
    end
end

ext: jl


text:
using DataFrames
using CSV
using Plots
#%%
data_s = open("sor_fluido.dat")
df_s=CSV.read(data_s, header=false, delim=" ")
x_s=df_s[1]
y_s=df_s[2]
z_s=transpose(reshape(df_s[3],(129, 65)))
x=range(0,stop=2,length=129)
y=range(0,stop=1,length=65)
#%%

plot(x,y,z_s, fill=true, label="sor")
xlabel!("x")
ylabel!("y")
title!("sor")
plot!(dpi=600)
savefig("sor_fluido.png")

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
module EnlilGrids

using MeshCore.Exports
using RandomTweaks
using StaticArrays
using Statistics

include("gmsh.jl")
export load_gmsh_file
export gmsh_do_physicalnames
export gmsh_do_entities
export gmsh_do_nodes
export gmsh_do_elements
export load_gmsh
export import_gmsh

end


text:
function mergesorted(outfname::AbstractString, infnames::AbstractString...)
    dfs = HDFSClient(TeraSort.DFSHOST, TeraSort.DFSPORT, UserGroupInformation(TeraSort.DFSUSER))
    for infname in infnames
        exists(dfs, infname) || (return 0)
    end
    fptrs = [open(dfs, infname, "r") for infname in infnames]
    wptr = open(dfs, outfname, "w")
    try
        mergesorted(wptr, fptrs...)
    finally
        for fptr in fptrs
            close(fptr)
        end
        close(wptr)
    end
    nothing
end

function mergesorted(wfile::IO, fptrs::IO...)
    println("merging sorted files...")
    inio = [fptrs...]
    nio = length(inio)
    lines = map(x->readline(x), inio)
    while length(inio) > 0
        (str,pos) = findmin(lines)
        write(wfile, str)
        if eof(inio[pos])
            splice!(inio, pos)
            splice!(lines, pos)
        else
            lines[pos] = readline(inio[pos])
        end
    end
    println("merged $(length(fptrs)) files")
    nothing
end

function merge(outfname::AbstractString, infnames::AbstractString...)
    dfs = HDFSClient(TeraSort.DFSHOST, TeraSort.DFSPORT, UserGroupInformation(TeraSort.DFSUSER))
    for infname in infnames
        exists(dfs, infname) || (return 0)
    end
    fptrs = [open(dfs, infname, "r") for fname in infnames]
    wptr = open(dfs, outfname, "w")
    try
        merge(wptr, fptrs...)
    finally
        for fptr in fptrs
            close(fptr)
        end
        close(wptr)
    end
    nothing
end

function merge(wptr::IO, fptrs::IO...)
    for fptr in fptrs
        d = readall(fptr)
        write(wptr, d)
        if d[end] != '\n'
            write(wptr, '\n')
        end
    end
    nothing
end


text:
answers = open("customs_answers.txt") do f
    map(row -> split(row, "\n"), split(read(f, String), "\n\n"))
end

println(sum(map(answer -> length(union(Set.(answer)...)), answers)))

println(sum(map(answer -> length(intersect(Set.(answer)...)), answers)))


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

mutable struct ConnectionMonitorResultProperties <: SwaggerModel
    source::Any # spec type: Union{ Nothing, ConnectionMonitorSource } # spec name: source
    destination::Any # spec type: Union{ Nothing, ConnectionMonitorDestination } # spec name: destination
    autoStart::Any # spec type: Union{ Nothing, Bool } # spec name: autoStart
    monitoringIntervalInSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: monitoringIntervalInSeconds
    endpoints::Any # spec type: Union{ Nothing, Vector{ConnectionMonitorEndpoint} } # spec name: endpoints
    testConfigurations::Any # spec type: Union{ Nothing, Vector{ConnectionMonitorTestConfiguration} } # spec name: testConfigurations
    testGroups::Any # spec type: Union{ Nothing, Vector{ConnectionMonitorTestGroup} } # spec name: testGroups
    outputs::Any # spec type: Union{ Nothing, Vector{ConnectionMonitorOutput} } # spec name: outputs
    notes::Any # spec type: Union{ Nothing, String } # spec name: notes

    function ConnectionMonitorResultProperties(;source=nothing, destination=nothing, autoStart=true, monitoringIntervalInSeconds=nothing, endpoints=nothing, testConfigurations=nothing, testGroups=nothing, outputs=nothing, notes=nothing)
        o = new()
        validate_property(ConnectionMonitorResultProperties, Symbol("source"), source)
        setfield!(o, Symbol("source"), source)
        validate_property(ConnectionMonitorResultProperties, Symbol("destination"), destination)
        setfield!(o, Symbol("destination"), destination)
        validate_property(ConnectionMonitorResultProperties, Symbol("autoStart"), autoStart)
        setfield!(o, Symbol("autoStart"), autoStart)
        validate_property(ConnectionMonitorResultProperties, Symbol("monitoringIntervalInSeconds"), monitoringIntervalInSeconds)
        setfield!(o, Symbol("monitoringIntervalInSeconds"), monitoringIntervalInSeconds)
        validate_property(ConnectionMonitorResultProperties, Symbol("endpoints"), endpoints)
        setfield!(o, Symbol("endpoints"), endpoints)
        validate_property(ConnectionMonitorResultProperties, Symbol("testConfigurations"), testConfigurations)
        setfield!(o, Symbol("testConfigurations"), testConfigurations)
        validate_property(ConnectionMonitorResultProperties, Symbol("testGroups"), testGroups)
        setfield!(o, Symbol("testGroups"), testGroups)
        validate_property(ConnectionMonitorResultProperties, Symbol("outputs"), outputs)
        setfield!(o, Symbol("outputs"), outputs)
        validate_property(ConnectionMonitorResultProperties, Symbol("notes"), notes)
        setfield!(o, Symbol("notes"), notes)
        o
    end
end # type ConnectionMonitorResultProperties

const _property_map_ConnectionMonitorResultProperties = Dict{Symbol,Symbol}(Symbol("source")=>Symbol("source"), Symbol("destination")=>Symbol("destination"), Symbol("autoStart")=>Symbol("autoStart"), Symbol("monitoringIntervalInSeconds")=>Symbol("monitoringIntervalInSeconds"), Symbol("endpoints")=>Symbol("endpoints"), Symbol("testConfigurations")=>Symbol("testConfigurations"), Symbol("testGroups")=>Symbol("testGroups"), Symbol("outputs")=>Symbol("outputs"), Symbol("notes")=>Symbol("notes"))
const _property_types_ConnectionMonitorResultProperties = Dict{Symbol,String}(Symbol("source")=>"ConnectionMonitorSource", Symbol("destination")=>"ConnectionMonitorDestination", Symbol("autoStart")=>"Bool", Symbol("monitoringIntervalInSeconds")=>"Int32", Symbol("endpoints")=>"Vector{ConnectionMonitorEndpoint}", Symbol("testConfigurations")=>"Vector{ConnectionMonitorTestConfiguration}", Symbol("testGroups")=>"Vector{ConnectionMonitorTestGroup}", Symbol("outputs")=>"Vector{ConnectionMonitorOutput}", Symbol("notes")=>"String")
Base.propertynames(::Type{ ConnectionMonitorResultProperties }) = collect(keys(_property_map_ConnectionMonitorResultProperties))
Swagger.property_type(::Type{ ConnectionMonitorResultProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorResultProperties[name]))}
Swagger.field_name(::Type{ ConnectionMonitorResultProperties }, property_name::Symbol) =  _property_map_ConnectionMonitorResultProperties[property_name]

function check_required(o::ConnectionMonitorResultProperties)
    true
end

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


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

mutable struct BackendAddressPoolPropertiesFormat <: SwaggerModel
    backendIPConfigurations::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceIPConfiguration2} } # spec name: backendIPConfigurations
    loadBalancerBackendAddresses::Any # spec type: Union{ Nothing, Vector{LoadBalancerBackendAddress} } # spec name: loadBalancerBackendAddresses
    loadBalancingRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: loadBalancingRules
    outboundRule::Any # spec type: Union{ Nothing, SubResource } # spec name: outboundRule
    outboundRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: outboundRules
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function BackendAddressPoolPropertiesFormat(;backendIPConfigurations=nothing, loadBalancerBackendAddresses=nothing, loadBalancingRules=nothing, outboundRule=nothing, outboundRules=nothing, provisioningState=nothing)
        o = new()
        validate_property(BackendAddressPoolPropertiesFormat, Symbol("backendIPConfigurations"), backendIPConfigurations)
        setfield!(o, Symbol("backendIPConfigurations"), backendIPConfigurations)
        validate_property(BackendAddressPoolPropertiesFormat, Symbol("loadBalancerBackendAddresses"), loadBalancerBackendAddresses)
        setfield!(o, Symbol("loadBalancerBackendAddresses"), loadBalancerBackendAddresses)
        validate_property(BackendAddressPoolPropertiesFormat, Symbol("loadBalancingRules"), loadBalancingRules)
        setfield!(o, Symbol("loadBalancingRules"), loadBalancingRules)
        validate_property(BackendAddressPoolPropertiesFormat, Symbol("outboundRule"), outboundRule)
        setfield!(o, Symbol("outboundRule"), outboundRule)
        validate_property(BackendAddressPoolPropertiesFormat, Symbol("outboundRules"), outboundRules)
        setfield!(o, Symbol("outboundRules"), outboundRules)
        validate_property(BackendAddressPoolPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type BackendAddressPoolPropertiesFormat

const _property_map_BackendAddressPoolPropertiesFormat = Dict{Symbol,Symbol}(Symbol("backendIPConfigurations")=>Symbol("backendIPConfigurations"), Symbol("loadBalancerBackendAddresses")=>Symbol("loadBalancerBackendAddresses"), Symbol("loadBalancingRules")=>Symbol("loadBalancingRules"), Symbol("outboundRule")=>Symbol("outboundRule"), Symbol("outboundRules")=>Symbol("outboundRules"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_BackendAddressPoolPropertiesFormat = Dict{Symbol,String}(Symbol("backendIPConfigurations")=>"Vector{NetworkInterfaceIPConfiguration2}", Symbol("loadBalancerBackendAddresses")=>"Vector{LoadBalancerBackendAddress}", Symbol("loadBalancingRules")=>"Vector{SubResource}", Symbol("outboundRule")=>"SubResource", Symbol("outboundRules")=>"Vector{SubResource}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ BackendAddressPoolPropertiesFormat }) = collect(keys(_property_map_BackendAddressPoolPropertiesFormat))
Swagger.property_type(::Type{ BackendAddressPoolPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BackendAddressPoolPropertiesFormat[name]))}
Swagger.field_name(::Type{ BackendAddressPoolPropertiesFormat }, property_name::Symbol) =  _property_map_BackendAddressPoolPropertiesFormat[property_name]

function check_required(o::BackendAddressPoolPropertiesFormat)
    true
end

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


text:
using Plots, LaTeXStrings, Printf, Dates
using JLD

include("tester.jl")

function time2str(t)
    ms = Millisecond(round(1000 * t))
    dt = convert(DateTime, ms)
    tt = convert(Time, dt)
    return string(tt)
end

function runtests(
    n,
    γ;
    gaussian_noise = false,
    tol = 1e-1,
    maxfgcalls = 100_000,
    useXold = true,
    printlevel = 0,
)
    methods = (:IR, :IER, :IAPG)
    results = Dict{Symbol, Vector{Float64}}()

    printlevel == 0 && @printf(
        "%4s %5s %7s %5s %6s %9s %9s %9s %7s\n",
        "n",
        "γ",
        "method",
        "k",
        "fgs",
        "rp",
        "rd",
        "ε",
        "time"
    )

    U, G, H, ncm = genprob(
        n,
        γ,
        gaussian_noise = gaussian_noise,
        maxfgcalls = maxfgcalls,
    )

    if useXold
        #X, y = CorNewton3(G)
        X = load("CorNewton3solns/n$n-γ$γ.jld", "X")
    end

    H2 = ncm.H2
    H2.data .= H .^ 2
    L = fronorm(H2, ncm.proj.work)

    for method in methods
        if method == :IR
            τ, α, σ = 0.95, 0.0, 1.0
        elseif method == :IER
            τ, α, σ = 1.0, 19/L, 1.0
        elseif method == :IAPG
            τ, α, σ = 1.0, 0.0, 1.0
        end

        if useXold
            ncm.Xold .= X
        end
        printlevel == 0 && @printf("%4d %5.2f %7s ", n, γ, method)
        t = @elapsed success, k = ncm(
            G,
            H,
            method = method,
            τ = τ,
            α = α,
            σ = σ,
            tol = tol,
            useXold = useXold,
            maxfgcalls = maxfgcalls,
            printlevel = printlevel,
        )
        fgcount = ncm.res.fgcountRef[]
        results[method] = ncm.res.resvals[1:fgcount]
        rp = ncm.res.rpRef[]
        rd = ncm.res.rdRef[]
        ε  = ncm.res.εRef[]
        if printlevel == 0
            @printf("%5d %6d %9.2e %9.2e %9.2e %7.1f",
                    k, fgcount, rp, rd, ε, t)
            success || @printf(" <----- FAILED")
            @printf("\n")
        end
    end

    return results
end

function makeplot(results, tol)

    plt = plot(
        yaxis = :log,
        ylims = [tol, 1e+2],
        xlabel = "function evaluations",
        ylabel = L"\max\{r_p,r_d\}",
        size = (900, 600),
        ls = :auto,
        lc = :black,
    )

    method = (:IR, :IER, :IAPG)
    label  = ("I-FISTA", "IE-FISTA", "IA-FISTA")

    for i = 1:length(method)
        res = results[method[i]]
        plot!(plt, res, label=label[i], ls=:auto, lc=:black)

        # Put an "x" on the plot to indicate a method died
        fgs = length(res)
        finalres = res[end]
        if finalres > tol
            plot!(plt, [fgs], [finalres], m=:x, mc=:black, label=false)
        end
    end

    return plt
end

function test(
    n,
    γ;
    gaussian_noise = false,
    tol = 1e-1,
    maxfgcalls = 100_000,
    useXold = true,
)
    results = runtests(
        n,
        γ,
        gaussian_noise = gaussian_noise,
        tol = tol,
        maxfgcalls = maxfgcalls,
        useXold = useXold,
    )

    plt = makeplot(results, tol)
    filename = @sprintf("n%d-γ%.2f.pdf", n, γ)
    savefig(plt, "../figs/$filename")

    return nothing
end

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

t = @elapsed begin
    for n = 100:100:800
        for γ = 0.1:0.1:1.0
            test(n, γ)
        end
    end
end
println(time2str(t))

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA


text:
#############################################################
#
# Friends & Followers Functions
#
#############################################################

function post_friendships_create(; options=Dict{AbstractString, AbstractString}())

    r = post_oauth("https://api.twitter.com/1.1/friendships/create.json", options)

    #Return array of type USERS
    return r.status == 200 ? to_USERS(Requests.json(r)) : r

end

function post_friendships_destroy(; options=Dict{AbstractString, AbstractString}())

    r = post_oauth("https://api.twitter.com/1.1/friendships/destroy.json", options)

    #Return array of type USERS
    return r.status == 200 ? to_USERS(Requests.json(r)) : r

end

function post_friendships_update(; options=Dict{AbstractString, AbstractString}())

    r = post_oauth("https://api.twitter.com/1.1/friendships/destroy.json", options)

    return r.status == 200 ? Requests.json(r) : r

end

function get_friends_list(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/friends/list.json", options)

    #Return users array
    return r.status == 200 ? to_USERS(Requests.json(r)["users"]) : r

end

function get_followers_list(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/followers/list.json", options)

    #Return users array
    return r.status == 200 ? to_USERS(Requests.json(r)["users"]) : r

end


text:
export GitTreeEntry,
       entry_byname, entry_bypath, walk

type GitTreeEntry{T<:GitObject}
    name::ByteString
    id::Oid
    filemode::Cint
    owns::Bool
end

Oid(te::GitTreeEntry)  = te.id
name(te::GitTreeEntry) = te.name

Base.isless(te1::GitTreeEntry, te2::GitTreeEntry) = isless(Oid(te1), Oid(te2))

Base.length(t::GitTree) =
    int(ccall((:git_tree_entrycount, libgit2), Csize_t, (Ptr{Void},), t))

Base.filemode(te::GitTreeEntry) = te.filemode

let
    function tree_entry_name(ptr::Ptr{Void})
        nptr = ccall((:git_tree_entry_name, libgit2), Ptr{UInt8}, (Ptr{Void},), ptr)
        return bytestring(nptr)
    end

    function tree_entry_oid(ptr::Ptr{Void})
        idptr = ccall((:git_tree_entry_id, libgit2), Ptr{UInt8}, (Ptr{Void},), ptr)
        return Oid(idptr)
    end

    function tree_entry_type(ptr::Ptr{Void})
        t = ccall((:git_tree_entry_type, libgit2), Cint, (Ptr{Void},), ptr)
        t == GitConst.OBJ_BLOB   && return GitBlob
        t == GitConst.OBJ_COMMIT && return GitCommit
        t == GitConst.OBJ_TAG    && return GitTag
        t == GitConst.OBJ_TREE   && return GitTree
        throw(ArgumentError("Unknown tree_entry_type constant $(t)"))
    end

    function tree_entry_filemode(ptr::Ptr{Void})
        ccall((:git_tree_entry_filemode, libgit2), Cint, (Ptr{Void},), ptr)
    end

    function GitTreeEntry(ptr::Ptr{Void}, owns::Bool=false)
        @assert ptr != C_NULL
        ty   = tree_entry_type(ptr)
        name = tree_entry_name(ptr)
        id   = tree_entry_oid(ptr)
        fm   = tree_entry_filemode(ptr)
        owns && ccall((:git_tree_entry_free, libgit2), Void, (Ptr{Void},), ptr)
        return GitTreeEntry{ty}(name, id, fm, owns)
    end
end

function entry_byname(t::GitTree, filename::AbstractString)
    eptr = ccall((:git_tree_entry_byname, libgit2), Ptr{Void},
                 (Ptr{Void}, Ptr{UInt8}), t, filename)
    return eptr == C_NULL ? nothing :
                            GitTreeEntry(eptr, false)
end

function entry_bypath(t::GitTree, path::AbstractString)
    entry_ptr = Ptr{Void}[0]
    @check ccall((:git_tree_entry_bypath, libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{UInt8}), entry_ptr, path)
    return entry_ptr[1] == C_NULL ? nothing :
                                    GitTreeEntry(entry_ptr[1], true)
end

function entry_byindex(t::GitTree, idx::Integer)
    entry_ptr = ccall((:git_tree_entry_byindex, libgit2), Ptr{Void},
                      (Ptr{Void}, Csize_t), t, idx - 1)
    return entry_ptr == C_NULL ? nothing :
                                 GitTreeEntry(entry_ptr, false)
end

function entry_byid(t::GitTree, id::Oid)
    entry_ptr = ccall((:git_tree_entry_byid, libgit2), Ptr{Void},
                      (Ptr{Void}, Ptr{Oid}), t, &id)
    return entry_ptr == C_NULL ? nothing :
                                 GitTreeEntry(entry_ptr, false)
end

Base.getindex(t::GitTree, entry::Integer) = entry_byindex(t, entry)
Base.getindex(t::GitTree, entry::AbstractString)  = entry_byname(t, entry)
Base.getindex(t::GitTree, entry::Oid)     = entry_byid(t, entry)

Base.start(t::GitTree) = begin
    te = entry_byindex(t, 1)
    return (1, te)
end

Base.done(t::GitTree, state) = state[1] > length(t)

Base.next(t::GitTree, state) = begin
    nidx = state[1] + 1
    (state[2], (nidx, entry_byindex(t, nidx)))
end

function foreach(::Type{GitBlob}, t::GitTree)
    @task for te in t
        isa(te, GitTreeEntry{GitBlob}) && produce(te)
    end
end

function foreach(::Type{GitTree}, t::GitTree)
    @task for te in t
        isa(te, GitTreeEntry{GitTree}) && produce(te)
    end
end

function cb_treewalk(root::Ptr{UInt8}, entry::Ptr{Void}, data::Ptr{Void})
    try
        produce(bytestring(root), GitTreeEntry(entry))
        return GitErrorConst.GIT_OK
    catch err
        return GitErrorConst.ERROR
    end
end

const c_cb_treewalk = cfunction(cb_treewalk, Cint,
                                (Ptr{UInt8}, Ptr{Void}, Ptr{Void}))

function walk(t::GitTree, order=:postorder)
    mode = order == :postorder ? GitConst.TREEWALK_POST :
           order == :preorder  ? GitConst.TREEWALK_PRE  :
           throw(ArgumentError("walk order can be :preorder or :postorder, got :$order"))
    @task ccall((:git_tree_walk, libgit2), Cint,
                (Ptr{Void}, Cint, Ptr{Void}, Ptr{Void}),
                t, mode, c_cb_treewalk, C_NULL)
end

function walk(::Type{GitBlob}, t::GitTree, order::Symbol=:postorder)
    @task for res in walk(t, order)
        isa(res[2], GitTreeEntry{GitBlob}) && produce(res)
    end
end

function walk(f::Function, ::Type{GitBlob}, t::GitTree, order::Symbol=:postorder)
    for res in walk(GitBlob, t, order)
        f(res)
    end
end

function walk(::Type{GitTree}, t::GitTree, order::Symbol=:postorder)
    @task for res in walk(t, order)
        isa(res[2], GitTreeEntry{GitTree}) && produce(res)
    end
end

function walk(f::Function, ::Type{GitTree}, t::GitTree, order::Symbol=:postorder)
    for res in walk(GitTree, t, order)
        f(res)
    end
end

#------- Tree Builder -------
type GitTreeBuilder
    ptr::Ptr{Void}
    repo::GitRepo

    function GitTreeBuilder(ptr::Ptr{Void}, r::GitRepo)
        @assert ptr != C_NULL
        tb = new(ptr, r)
        finalizer(tb, free!)
        return tb
    end
end

GitTreeBuilder(r::GitRepo) = begin
    tbptr = Ptr{Void}[0]
    @check ccall((:git_treebuilder_create, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), tbptr, C_NULL)
    return GitTreeBuilder(tbptr[1], r)
end

free!(tb::GitTreeBuilder) = begin
    if tb.ptr != C_NULL
        ccall((:git_treebuilder_free, libgit2), Void, (Ptr{Void},), tb.ptr)
        tb.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, tb::GitTreeBuilder) = tb.ptr

# the number of entries in a tree builder
Base.length(tb::GitTreeBuilder) = begin
    return int(ccall((:git_treebuilder_entrycount, libgit2), Cuint, (Ptr{Void},), tb))
end

# add / update an entry to the builder
Base.insert!(tb::GitTreeBuilder, filename::AbstractString, id::Oid, filemode::Int) = begin
    if !(id in tb.repo)
        throw(ArgumentError("OID $id does not exist in the Object Database"))
    end
    @check ccall((:git_treebuilder_insert, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{UInt8}, Ptr{Oid}, Cint),
                 C_NULL, tb, filename, &id, filemode)
    return tb
end

# get an entry from the builder from its filename
Base.getindex(tb::GitTreeBuilder, filename::AbstractString) = begin
    entry_ptr = ccall((:git_treebuilder_get, libgit2), Ptr{Void},
                      (Ptr{Void}, Ptr{UInt8}), tb, filename)
    return GitTreeEntry(ptr)
end

function cb_tbfilter(te_ptr::Ptr{Void}, payload::Ptr{Void})
    te = GitTreeEntry(te_ptr)
    func = unsafe_pointer_to_objref(payload)::Vector{Any}
    return bool(func(te)) ? Cint(1) : Cint(0)
end

const c_cb_tbfilter = cfunction(cb_tbfilter, Cint, (Ptr{Void}, Ptr{Void}))

# selectively remote entries in the tree
Base.filter!(f::Function, tb::GitTreeBuilder) = begin
    f_ptr = pointer_from_objref(f)
    ccall((:git_treebuilder_filter, libgit2), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{Void}), tb, c_cb_tbfilter, f_ptr)
    return tb
end

# clear all entries in the builder
clear!(tb::GitTreeBuilder) =
    ccall((:git_treebuilder_clear, libgit2), Void, (Ptr{Void},), tb)

# write the contents of the tree builder as a tree object
function write!(tb::GitTreeBuilder)
    id_ptr = [Oid()]
    @check ccall((:git_treebuilder_write, libgit2), Cint,
                 (Ptr{Oid}, Ptr{Void}, Ptr{Void}), id_ptr, tb.repo, tb)
    return id_ptr[1]
end


text:
export AbstractWrappedArray, basetype
export WTranspose, WAdjoint, WSymmetric, WUpperTriangular

abstract type AbstractWrappedArray{T,N,P,S} <: AbstractArray{T,N}; end

basetype(::Type{P}) where P<:AbstractWrappedArray{T,N,X,S} where {T,N,X,S} = S
basetype(::Type{P}) where P = P
basetype(x::T) where T = basetype(T)

SparseArrays.issparse(::Type{P}) where P<:AbstractWrappedArray = issparse(basetype(P))
SparseArrays.issparse(::Type{<:Union{SparseVector,SparseMatrixCSC}}) = true
SparseArrays.issparse(::Type) = false
SparseArrays.issparse(x::T) where T<:AbstractArray = issparse(T)

# the following concrete wrapper types are demo replacements for the standard ones.
# an additional type parameter is appended to the original
# actually the original definitions in LinearAlgebra and Base have to be modified

# Transpose
struct WTranspose{T,P,S} <:AbstractWrappedArray{T,2,P,S}
    parent::P
    WTranspose(A::P) where P<:AbstractArray{T,N} where {T,N} = new{T,P,basetype(P)}(A)
end
Base.size(x::WTranspose) = reverse(size(x.parent))
Base.getindex(x::WTranspose,I,J,K...) = getindex(x.parent, J, I, K...)

# Adjoint
struct WAdjoint{T,P,S} <:AbstractWrappedArray{T,2,P,S}
    parent::P
    WAdjoint(A::P) where P<:AbstractArray{T,N} where {T,N} = new{T,P,basetype(P)}(A)
end
Base.size(x::WAdjoint) = reverse(size(x.parent))
Base.getindex(x::WAdjoint,I,J,K...) = adjoint(getindex(x.parent, J, I, K...))

# Symmetric
struct WSymmetric{T,P,S} <:AbstractWrappedArray{T,2,P,S}
    parent::P
    uplo::Char
    WSymmetric(A::P,ul::Symbol) where P<:AbstractArray{T,N} where {T,N} = new{T,P,basetype(P)}(A, ul == :U ? 'U' : 'L')
end
Base.size(x::WSymmetric) = size(x.parent)
Base.getindex(x::WSymmetric,I...) = getindex(Symmetric(x.parent, x.uplo == 'U' ? :U : :L), I...)

# UpperTriangular
struct WUpperTriangular{T,P,S} <:AbstractWrappedArray{T,2,P,S}
    parent::P
    UpperTriangular(A::P) where P<:AbstractArray{T,N} where {T,N} = new{T,P,basetype(P)}(A)
end
Base.size(x::WUpperTriangular) = size(x.parent)
Base.getindex(x::WUpperTriangular,I...) = getindex(UpperTriangular(x.parent), I...)

# SubArray


text:
include("./Compressions.jl")
export Algebraic, Exponential, Trigonometric

include("./Couplings.jl")
export generateCoupling, modeContribution, modifiedMode

include("./ModeFunctions.jl")
export Mode, FourierMode, HardBoxMode, SoftBoxMode, SoftBoxExpMode, GaussMode, ExpMode

include("./WavePacketStates.jl")
export WavePacketState, Displaced, NonDisplaced,
  Coherent, ArbitraryState, Fock, SqueezedVacuum, createState,
  DisplacedArbitraryState, DisplacedFock

abstract type _WavePacket{T <: WavePacketState} end

"""
    WavePacket(mode::Mode, state::WavePacketState)

Describes the _input_ wave packet.
A `WavePacket` consists of any [`Mode`](@ref) and any [`WavePacketState`](@ref).
"""
mutable struct WavePacket{S} <: _WavePacket{S}
  mode::Mode
  state::S
end

"""
    ContinuousWave(α::Union{Number, Function})

A special form of a [`Coherent`](@ref) state. For constant `α` this allows to solver to perform additional
optimisations. Additionally, expresses the intend of using a non-normalisable input mode.
"""
mutable struct ContinuousWave <: _WavePacket{Coherent}
  α::Number
end

function ContinuousWave(α::Function)
	warning = "`ContinuousWave(α::Function)` doesn't define gᵢ, gₒ, or norm, " *
		"use `WavePacket(Mode(α), Coherent(1.0))` instead"
  mode = Mode(α, # Since we do not want to assume `α` to be normalised.
    t -> begin @warn(warning); NaN end,
    t -> begin @warn(warning); NaN end,
		t -> begin @warn(warning); NaN end
  )
  WavePacket(mode, Coherent(1.0))
end

NoPacket() = ContinuousWave(0.0)


text:
"""
%%
% CMG, Copyright (c) 2008-2020  Ioannis Koutis, Gary Miller
%
% The CMG solver is distributed under the terms of the GNU General Public
% Lincense Version 3.0 of the Free Software Foundation.
% CMG is also available under other licenses; contact authors for details.

%% List of functions in this file:
%
% cmg_precondition:    main function
%
% validate_input
% make_preconditioner
%
% Combinatorial Functions
% steiner_group
% forest_components_
% split_forest_
%
% Numerical Functions
% ldl_
% preconditioner_
%sd

"""

## Prerequisite
# import Pkg
# Pkg.add("LinearAlgebra", "SparseArrays", "Laplacians","Statistics")
# using LinearAlgebra, SparseArrays, Laplacians, Statistics

"""
   cmg_precondition(A::SparseMatrixCSC)

For an abritrary b-side in the null space of A, the
system Ax = b can be solved by:
   x = pcg(A, b, tol, iter, pfun)

# Arguments
- `A::SparseMatrixCSC`: SDD matrix (positive off-diagonals not supported currently)

# Outputs
- `pfun`: Preconditioning function for A
- `H`: Preconditioning hierarchy used internally in pfun
- `flag`: failure status variable

# Example
include("cmg_precondition.jl")
A = lap(grid2(100,100))
@time pfun, H, flag = cmg_precondition(A);
"""

mutable struct Chol
    ld::SparseMatrixCSC
    ldT::SparseMatrixCSC
    d::Array{Float64,1}
    p::Array{Int64,1}
    invp::Array{Int64,1}
end

Chol() = Chol(sparse([]), sparse([]), sparse([]), [], [])

mutable struct HierarchyLevel
    sd::Bool
    islast::Bool
    iterative::Bool
    A::SparseMatrixCSC
    invD::Array{Float64}
    cI::Array{Int}
    nc::Int
    R::SparseMatrixCSC
    repeat::Int
    chol::Chol
end

HierarchyLevel() = HierarchyLevel(true, 0, 1, sparse([]), [], [], 0, sparse([]), 0, Chol())

function cmg_precondition(A::SparseMatrixCSC)

    msg = undef
    flag = undef

    if size(A, 1) < 500 # handle small input
        msg = "Input matrix is small. Solve Ax=B with A\\b."
    else # validate input
        (A_, flag) = validate_input(A)

        if flag == 1
            msg = "The input matrix must be symmetric."
        elseif flag == 2
            msg = "The current version of CMG does not support positive off-diagonals."
        end
    end

    if typeof(msg) != UndefInitializer
        println("\n", msg, "\n")
        pfun = undef;
        H = undef;
        return pfun, H, flag
    end

   # initialize H
    H = HierarchyLevel[]

   # initialize a hierarchy level

    push!(H, HierarchyLevel())

   # A_ is Laplacian
    if size(A_, 1) > size(A, 1)
        H[1].sd = true # original matrix is strongly dominant
    else
        H[1].sd = false
    end

   # construct hierarchy
    loop = true
    j = 1
    h_nnz = 0
    n = 0 # initialize n
    # build up H
    while loop
        n = size(A_, 1)

      # direct method for small size
        if (n < 500)
            H[j].iterative = 0
            break
        end

        dA_ = diag(A_)
        (cI, ~) = steiner_group(A_, dA_)
        nc = maximum(cI)
        H[j].islast = 0
        H[j].A = A_ # !
        H[j].invD = 1 ./ (2 * dA_) # !
        H[j].cI = cI
        H[j].nc = nc
        H[j].R = sparse(cI, 1:n, ones(n), nc, n) # added for efficiency

        # check for full contraction
        if nc == 1
            H[j].islast = 1
            H[j].iterative = 1
        end

        # check for hierarchy stagnation for potentially bad reasons
        h_nnz = h_nnz + nnz(A_)
        if (nc >= n - 1) || (h_nnz > 5 * nnz(H[1].A))
            H[j].islast = 1
            H[j].iterative = 1
            flag = 3 # indicates stagnation
            @warn "CMG convergence may be slow due to matrix density. Future versions of CMG will eliminate this problem."
            break
        end

        Rt = sparse(cI, 1:n, 1, nc, n) # ! take out double
        A_ = Rt * (H[j].A) * Rt'

        j += 1
        push!(H, HierarchyLevel())
    end

    # code for last hierarchy level
    if flag == 0
        H[j].islast = 1
        H[j].iterative = 0
        H[j].A = A_[1:n - 1,1:n - 1]
      # L, D, p = ldl(H[j].A, 'vector')
        L, D, p = ldl_(H[j].A)
        H[j].chol.ld = L
        H[j].chol.ldT = L'
        H[j].chol.d = (1 / diag(D))' # ! 1.?
        H[j].chol.p = p # x = A*y => y(p) = LT\(H{j}.d.*(L\x(p)));
    end

    # determine number of recursive calls between levels
    for k in 1:(j - 2)
        H[k].repeat = max(floor(nnz(H[k].A) / nnz(H[k + 1].A) - 1), 1)
    end

    H[j - 1].repeat = max(floor(nnz(H[j - 1].A) / nnz(H[j].chol.ld) - 1), 1)

   # H = cell2mat(H);

   # create precondition function
    pfun = make_preconditioner(H)
    return pfun, H, flag
end

"""
   function (A0, flag) = validate_input(A::SparseMatrixCSC)
Input validation.

# Arguments
-`A`: Laplacian, possibly extended by one coordinate when A is SD

# Outputs
-`A0`
-`flag`: failure flag code indicating issue with input validity
"""

function validate_input(A::SparseMatrixCSC)
    flag = 0

   # check symmetry
    if !issymmetric(A)
        flag = 1
    end

   # detect strict dominance
    n = size(A, 1)
    dA = diag(A)
    sA = sum(A, dims = 2)
    sAp = (sA + abs.(sA)) / 2
    sd = sAp ./ dA .> 1e-13 # if sd(i) = 1 then row i is numerically strictly dominant

   # A is a sparse matrix
    (i, j, v) = findnz(A)

   # check for positive off-diagonals that are not currently supported
    if maximum(v[i .!= j]) > 0
        flag = 2
    end

   # augment by extra coordinate if strictly dominant
    if maximum(sd)
        ex_v = -sAp[sd]
        exd = length(ex_v)
        exd_c = nonzeros(sd) # get coordinates
        i = Int[i;(n + 1) * ones(exd); exd_c; n + 1]
        j = Int[j;exd_c; (n + 1) * ones(exd);n + 1]
        v = [v; ex_v; ex_v; -sum(ex_v)]
        A0 = sparse(i, j, v, n + 1, n + 1)
    else
        A0 = A
    end

    return A0, flag
end

"""
   function(cI, nc) = steiner_group(A, dA_)
Steiner groups
# Arguments
-`A`: Laplacian
"""
function steiner_group(A, dA_)

    (M, c1) = findmin(A, 1) # c1 represents a unimodal forest #! dims = 1

   # adjust c1 to array of indices
    C = Array{Int,1}(undef, length(c1))

    for i in 1:length(c1)
        C[i] = c1[i][1]
    end

    split_forest_!(C)

    efd = abs.(M' ./ dA_)

    if minimum(efd) < 1 / 8 # low effective degree nodes found
        C = update_groups_(A, C, dA_)
    end
    cI, nc, ~ = forest_components_(C)
end

"""
   function C = split_forest_(C1::Array{Int})
decompose unimodal forest into low conductance components
"""
function split_forest_!(C::Array{Int64})
    n = length(C)
    ancestors = zeros(Int, n)
    indegree = zeros(Int, n + 2)
    visited = falses(n) # logical sparse array

    walkbuffer = zeros(Int, 20);
    newancestorbuff = zeros(Int, 20)

   # compute indegrees
    for i in 1:n
        indegree[C[i]] = indegree[C[i]] + 1
    end

   # partition into clusters of small diameter

    for j in 1:n
        jwalk = j
        startwalk = true

        while (startwalk && (indegree[jwalk] == 0 ) && !visited[jwalk])
            startwalk = false
            ancestors_in_path = 0 # change over C-CMG
            k = 1
            walkbuffer[k] = jwalk
            newancestorbuff[k] = 0
            while (k <= 6 || visited[jwalk]) # +1 for c-indexing adjust
                jwalk = C[jwalk]
                walkterminated = (jwalk == walkbuffer[k]) || ((k > 1) && (jwalk == walkbuffer[k - 1]))
                if walkterminated
                    break; # while
                end
                k = k + 1
                walkbuffer[k] = jwalk
                if visited[jwalk]
                    newancestorbuff[k] = ancestors_in_path
                else
                    ancestors_in_path = ancestors_in_path + 1
                    newancestorbuff[k] = ancestors_in_path
                end
            end

            if k > 6 # large diameter - cut
                middlek = Int(ceil(k / 2))
                C[walkbuffer[middlek]] = walkbuffer[middlek] # cut middle edge
                indegree[walkbuffer[middlek + 1]] = indegree[walkbuffer[middlek + 1]] - 1 # update indegree

                for ik in (middlek + 1):k
                    ancestors[walkbuffer[ik]] = ancestors[walkbuffer[ik]] - ancestors[walkbuffer[middlek]]
                end

            # update ancestors and visited flag
                for ik in 1:middlek
                    visited[walkbuffer[ik]] = true
                    ancestors[walkbuffer[ik]] = ancestors[walkbuffer[ik]] + newancestorbuff[ik]
                end

            # set first vertex in new walk
                jwalk = walkbuffer[middlek + 1]
                startwalk = true
            end # end cut procedure

         # commit walk changes
            if !startwalk
                for ik in 1:k
                    ancestors[walkbuffer[ik]] = ancestors[walkbuffer[ik]] + newancestorbuff[ik]
                    visited[walkbuffer[ik]] = true
                end
            end
        end # outer while
    end

   # tree partition into clusters of high conductance
    for j in 1:n
        jwalk = j
        startwalk = true

        while startwalk && (indegree[jwalk] == 0)
            startwalk = false
            jwalkb = jwalk
            cut_mode = false
            # initialize new_front
            new_front = 0
            removed_ancestors = []

            while true
                jwalka = C[jwalk]
                walkterminated = (jwalka == jwalk) || (jwalka == jwalkb)
                if walkterminated
                    break; # while
                end

                if (!cut_mode && (ancestors[jwalk] > 2) && (ancestors[jwalka] - ancestors[jwalk] > 2)) # possibly low conductance - make cut
                    C[jwalk] = jwalk # cut edge
                    indegree[jwalka] = indegree[jwalka] - 1
                    removed_ancestors = ancestors[jwalk]
                    new_front = jwalka
                    cut_mode = true
                end # end making cut

                jwalkb = jwalk
                jwalk = jwalka
                if cut_mode
                    ancestors[jwalk] = ancestors[jwalk] - removed_ancestors
                end
            end
            if cut_mode
                startwalk = true
                jwalk = new_front
            end
        end
    end
end # split_forest_

"""
   function C1 = update_groups_(A, C, dA_)
update groups based on nodes with low effective degree
# Arguments
-`A`: Laplacian
"""
function update_groups_(A, C, dA_)
    n = size(C, 1)
    B = zeros(n)

   # B[j] is the total tree weight incident to node j
    for i in 1:n
        if C[i] != i
            B[i] = A[i,C[i]] + B[i]
            B[C[i]] = A[i,C[i]] + B[C[i]]
        end
    end

    ndx = findall(x->x > 0.125, B ./ dA_)
    C[ndx] = Array{Int32}(ndx)

    C1 = C

    return C1
end # update_groups_

"""
   function[cI, nc, csizes] = forest_components_(C)
forest components, connected components in unimodal forest
# Arguments
-`C`: unimodal tree
"""

function forest_components_(C)
    n = size(C, 1)
    cI = zeros(Int, n)
    cSizes = zeros(Int, n)
    buffer = zeros(Int, 100)

    ccI = 1
    for j in 1:n
        bufferI = 1
        jwalk = j

      # tree walk
        while cI[jwalk] == 0
            cI[jwalk] = ccI
            buffer[bufferI] = jwalk
            bufferI = bufferI + 1

            if bufferI == size(buffer, 1) # for memory efficiency
                buffer_ = zeros(Int, min(2 * size(buffer, 1), n))
                buffer_[1:size(buffer, 1)] = buffer
                buffer = buffer_
            end

            jwalk = C[jwalk]
        end # while

        bufferI = bufferI - 1
        en = C[jwalk] # end node
        if cI[en] != ccI
            cI[buffer[1:bufferI]] .= cI[en] # ! .=
        else
            ccI = ccI + 1
        end
        cSizes[en] = cSizes[en] + bufferI
    end # for
    if cSizes[ccI] == 0
        ccI = ccI - 1
    end

    nc = ccI
    cSizes = cSizes[1:ccI]

    return cI, nc, cSizes
end

"""
   function pfun = make_preconditioner(H)
Make preconditioner
"""
function make_preconditioner(H::Array{HierarchyLevel,1})

    if !H[1].sd
        pfun = b->preconditioner_(b, H, 1)
        return pfun
    end

    if H[1].sd
        pfun = b->preconditioner_sd(b, H)
        return pfun
    end

end

"""
   function x = preconditioner_(b, H::Array{HierarchyLevel,1}, level)
preconditioner
"""
function preconditioner_(b::Array{Float64}, H::Array{HierarchyLevel,1}, level::Int64)
    n = size(H[level].A, 1)

   # base case I: last level direct
    if (H[level].islast && !H[level].iterative)
        b = b[1:n]
        x = zeros(Float64, n + 1)

        L = H[level].chol.ld
        LT = H[level].chol.ldT
        di = H[level].chol.d
        p = H[level].chol.p

      # back substitution needed here
        x[p] = LT \ (di .* (L \ b[p])) # ! .*
        return x
    end

   # base case II: last level iterative
    if (H[level].islast && H[level].iterative)
        invD = H[level].invD
        x = invD .* b
        return x
    end

   # main cycle
   # unpack variables for efficiency
    invD = H[level].invD
    A = H[level].A
    repeat = H[level].repeat
    cI = H[level].cI
    R = H[level].R

    r = Array{Float64}(undef, size(A, 1))
    b_small = Array{Float64}(undef, size(R, 1))
    z = Array{Float64}(undef, size(R, 1))

    for j in 1:repeat
      # jacobi pre-smooth
        if j == 1
            x = invD .* b # previous x = 0
        else
            x = x + invD .* (b - A * x)
        end

      # residual
        r = b - A * x

      # projection of residual onto smaller level
      # b_small = spzeros(nc)
      # for k in 1:n
      #     b_small[cI[k]] = b_small[cI[k]] + r[k]
      # end
        b_small = R * r

        z = preconditioner_(b_small, H, level + 1)

      # interpolation
        x = x + z[cI]    #z[cI] does interpolation

      # Jacobi post-smooth
        x = x + invD .* (b - A * x)
    end
    return x
end

"""
   function x = preconditioner_sd(H::Array{HierarchyLevel}, b)
preconditioner sd
"""
function preconditioner_sd(b::Array{Float64}, H::Array{HierarchyLevel})
    n = length(b)
    push!(b, -sum(b))

    x = b->preconditioner_(b, H, 1)

    x = x[1:n] + x[n + 1]
    return x
end

"""
    function[L, D, p] = ldl_(A)
ldl Factorization with min degree heuristic
"""
function ldl_(A::SparseMatrixCSC)
    n = size(A, 1)
    D = zeros(n)

    # trivial min degree heuristic: sorting degrees
    (iA, jA, ~) = findnz(A)
    Ao = sparse(iA, jA, ones(length(iA)), n, n)
    dA = sum(Ao, dims = 1)
    p = sortperm(dA[:])

    A = A[p,p]
    L = sparse(1:n, 1:n, ones(n)) # identity matrix
    for i in 1:n
        vi = A[i + 1:n,i]
        di = A[i,i]
        D[i] = di
        L[i + 1:n,i] = vi / di
        Z = vi * vi' / di
        (iz, jz, vz) = findnz(Z)
        Z = sparse(iz .+ i, jz .+ i, vz, n, n)
        A = A - Z
    end

    D = sparse(1:n, 1:n, D, n, n)
    return L, D, p
end


text:
using QXTns
using Test
using TestSetExtensions

@testset ExtendedTestSet "QXTns.jl" begin
    @includetests ARGS
end


text:
# Closest Correlation matrix test matrix

# Original problem has the following format:
# min_X   1/2 ||X-C||^2
# s.t.    Xii = 1
#         X ⪴ 0

using COSMO, Test, LinearAlgebra, SparseArrays, Random

# this function creates a matrix A that slices out the diagonal entries Xii of a vectorized square matrix x=vec(X)
function create_diagonal_extractor(n)
  A = spzeros(n, n^2)
  A[1,1 ] = 1
  for iii = 2:n-1
    col = (iii - 1) * (n + 1)
    A[iii, col + 1] = 1
  end
  A[n, n^2] = 1
  return A
end

function diagonal_equal_one(A, atol)
  for iii = 1:size(A, 1)
    (abs(A[iii, iii] - 1) > atol) && return false
  end
  return true
end

rng = Random.MersenneTwister(12345)
xMin = -1.
xMax = 1.
n = 50
C = xMin .+ randn(rng, n, n)*(xMax - xMin)
c = vec(C)

n2 = n^2
m = n + n2

P = sparse(1.0I, n2, n2)
q = -vec(C)
r = 0.5*vec(C)' * vec(C)

A1 = create_diagonal_extractor(n)
b1 = -ones(n)
cs1 = COSMO.Constraint(A1, b1, COSMO.ZeroSet)

A2 = sparse(1.0I, n2, n2)
b2 = zeros(n2)
cs2 = COSMO.Constraint(A2, b2, COSMO.PsdCone)
constraints = [cs1; cs2]

model = COSMO.Model()
assemble!(model, P, q, constraints)

res = COSMO.optimize!(model)

Xsol = reshape(res.x, n, n)

@testset "Closest Correlation Matrix - SDP Problems" begin
  @test res.status == :Solved
  @test diagonal_equal_one(Xsol, 1e-5)
  @test minimum(eigen(Xsol).values) > -1e-3

end
nothing


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    int xTESTOUTPUT = 7;
    int wTESTOUTPUT = 7;
    void m() {
      final Condition c = new Condition();
      int y = ({L /c H}) 42;
      int z = 0;
      if (y > 0) {
	  c.set();
	  z = 4; // level of z should be H, since the PC should have been raised since c has been set.
      }
      this.xTESTOUTPUT = y; // At this point, the level of y should be H, since c has been set.
      this.wTESTOUTPUT = z;
    }
    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
"""
n=define(;numControls=1,X0=[0.,1],XF=[0.,-1.],XL=[0.,-Inf],XU=[1/9,Inf],CL=[-Inf],CU=[Inf])
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/1/2017, Last Modified: 4/13/2018 \n
Citations: \n
----------\n
Initially Influenced by: S. Hughes.  steven.p.hughes@nasa.gov
Source: DecisionVector.m [located here](https://sourceforge.net/p/gmat/git/ci/264a12acad195e6a2467cfdc68abdcee801f73fc/tree/prototype/OptimalControl/LowThrust/@DecisionVector/)
-------------------------------------------------------------------------------------\n
"""
function define(;
                numStates::Int64 = 0,
                numControls::Int64 = 0,
                X0 = fill(NaN,numStates,),
                XF = fill(NaN,numStates,),
                XL = fill(NaN,numStates,),
                XU = fill(NaN,numStates,),
                CL = fill(NaN,numControls,),
                CU = fill(NaN,numControls,),
                XS = fill(1.,numStates,),
                CS = fill(1.,numControls,)
                )

  n = NLOpt()

  # validate input
  if numControls <= 0
      error("numControls must be > 0","\n")
  end
  if numStates <= 0
      error("numStates must be > 0","\n")
  end
  if length(X0) != numStates
    error(string("\n Length of X0 must match number of states \n"))
  end
  if length(XF) != numStates
    error(string("\n Length of XF must match number of states \n"))
  end
  if length(XL) != numStates
    error(string("\n Length of XL must match number of states \n"))
  end
  if length(XU) != numStates
    error(string("\n Length of XU must match number of states \n"))
  end
  if length(XS) != numStates
    error(string("\n Length of XS must match number of states \n"))
  end
  if length(CL) != numControls
    error(string("\n Length of CL must match number of controls \n"))
  end
  if length(CU) != numControls
    error(string("\n Length of CU must match number of controls \n"))
  end
  if length(CS) != numControls
    error(string("\n Length of CS must match number of controls \n"))
  end

  n.ocp.state = initState(numStates)
  n.ocp.control = initControl(numControls)
  n.ocp.X0 = X0
  n.ocp.X0_tol = NaN*X0
  n.ocp.XF = XF
  n.ocp.XF_tol = NaN*XF
  n.ocp.XL = XL
  n.ocp.XU = XU
  n.ocp.CL = CL
  n.ocp.CU = CU
  n.ocp.XS = XS
  n.ocp.CS = CS
  n.f.ocp.defined = true

  return n
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/9/2017, Last Modified: 12/02/2019 \n
Last Modified by: Alexander Buck
-------------------------------------------------------------------------------------\n
"""
function defineSolver!(n::NLOpt,kw)

 if typeof(kw)!=Dict{Symbol,Symbol}
   kw = Dict(kw)
 end

  # get the name of the solver
  if haskey(kw,:name); n.s.ocp.solver.name=get(kw,:name,0); end
  if try_import(n.s.ocp.solver.name)
  else error(string("could not import ",n.s.ocp.solver.name) )
  end

  # see if the user would like to use a standard set of solver settings for mpc
  if haskey(kw,:mpc_defaults); mpc_defaults=get(kw,:mpc_defaults,0); else mpc_defaults=false; end
  if mpc_defaults
    if n.s.ocp.solver.name==:Ipopt
      n.s.ocp.solver.settings=_Ipopt_MPC;
    elseif n.s.ocp.solver.name==:KNITRO
      n.s.ocp.solver.settings=_KNITRO_MPC;
    else
      error(string("solver ",n.s.sover.name, " not defined"))
    end
  else # default solver settings
    if n.s.ocp.solver.name==:Ipopt  # NOTE this should already have been done by default, but could get messed up is user is playing with options
      n.s.ocp.solver.settings=_Ipopt_defaults;
    elseif n.s.ocp.solver.name==:KNITRO
      n.s.ocp.solver.settings=_KNITRO_defaults;
    else
      error(string("solver ", n.s.sover.name, " not defined"))
    end
  end

  # modify additional defaults individually
  for (key,value) in kw
    if haskey(n.s.ocp.solver.settings,key)
      n.s.ocp.solver.settings[key]=value
    elseif key!=:name && key!=:mpc_defaults # ignore the name and default settings option TODO could remove them from the Dict
      error(string(" \n Unknown key: ", kw, " for ", n.s.ocp.solver.name, " used in defineSolver!() \n "))
    end
  end

  if n.s.ocp.solver.name==:Ipopt
    setsolver(n.ocp.mdl,Ipopt.IpoptSolver(;max_cpu_time=n.s.ocp.solver.settings[:max_cpu_time],
                               print_level=n.s.ocp.solver.settings[:print_level],
                               warm_start_init_point=n.s.ocp.solver.settings[:warm_start_init_point],
                               max_iter=n.s.ocp.solver.settings[:max_iter],
                               tol=n.s.ocp.solver.settings[:tol],
                               dual_inf_tol=n.s.ocp.solver.settings[:dual_inf_tol],
                               constr_viol_tol=n.s.ocp.solver.settings[:constr_viol_tol],
                               compl_inf_tol=n.s.ocp.solver.settings[:compl_inf_tol],
                               acceptable_tol=n.s.ocp.solver.settings[:acceptable_tol],
                               acceptable_constr_viol_tol=n.s.ocp.solver.settings[:acceptable_constr_viol_tol],
                               acceptable_dual_inf_tol=n.s.ocp.solver.settings[:acceptable_dual_inf_tol],
                               acceptable_compl_inf_tol=n.s.ocp.solver.settings[:acceptable_compl_inf_tol],
                               acceptable_obj_change_tol=n.s.ocp.solver.settings[:acceptable_obj_change_tol],
                               diverging_iterates_tol=n.s.ocp.solver.settings[:diverging_iterates_tol]))
  elseif n.s.ocp.solver.name==:KNITRO
    setsolver(n.ocp.mdl,KnitroSolver(;outlev=n.s.ocp.solver.settings[:outlev],
                                 feastol=n.s.ocp.solver.settings[:feastol],
                                 feastol_abs=n.s.ocp.solver.settings[:feastol_abs],
                                 ftol=n.s.ocp.solver.settings[:ftol],
                                 ftol_iters=n.s.ocp.solver.settings[:ftol_iters],
                                 infeastol=n.s.ocp.solver.settings[:infeastol],
                                 maxfevals=n.s.ocp.solver.settings[:maxfevals],
                                 maxit=n.s.ocp.solver.settings[:maxit],
                                 maxtime_cpu=n.s.ocp.solver.settings[:maxtime_cpu],
                                 maxtime_real=n.s.ocp.solver.settings[:maxtime_real],
                                 opttol=n.s.ocp.solver.settings[:opttol],
                                 opttol_abs=n.s.ocp.solver.settings[:opttol_abs],
                                 xtol=n.s.ocp.solver.settings[:xtol],
                                 xtol_iters=n.s.ocp.solver.settings[:xtol_iters],
                                 algorithm=n.s.ocp.solver.settings[:algorithm],
                                 bar_initpt=n.s.ocp.solver.settings[:bar_initpt],
                                 bar_murule=n.s.ocp.solver.settings[:bar_murule],
                                 bar_penaltycons=n.s.ocp.solver.settings[:bar_penaltycons],
                                 bar_penaltyrule=n.s.ocp.solver.settings[:bar_penaltyrule],
                                 bar_switchrule=n.s.ocp.solver.settings[:bar_switchrule],
                                 linesearch=n.s.ocp.solver.settings[:linesearch],
                                 linsolver=n.s.ocp.solver.settings[:linsolver],
                                 tuner=n.s.ocp.solver.settings[:tuner],
                                 cg_pmem=n.s.ocp.solver.settings[:cg_pmem]))
  else
    error(string("solver ",n.s.sover.name, " not defined"))
  end
  return nothing
end  # function

"""
OCPdef!(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/14/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function OCPdef!(n::NLOpt)
  # state variables
  @variable(n.ocp.mdl,x[1:n.ocp.state.pts,1:n.ocp.state.num])
  n.r.ocp.x = Matrix{Any}(undef, n.ocp.state.pts,n.ocp.state.num)
  for st in 1:n.ocp.state.num
      n.r.ocp.x[:,st] = @NLexpression(n.ocp.mdl, [j=1:n.ocp.state.pts], n.ocp.XS[st]*x[j,st])
  end
#  n.r.ocp.x=x;

  for st in 1:n.ocp.state.num
    # lower state constraint
    if !isnan.(n.ocp.XL[st])
      if n.ocp.mXL[st] == false
        for j in 1:n.ocp.state.pts
          setlowerbound(x[j,st], n.ocp.XL[st]/n.ocp.XS[st])
        end
      else
        for j in 1:n.ocp.state.pts
          setlowerbound(x[j,st], n.ocp.XL_var[st,j]/n.ocp.XS[st])
        end
      end
    end

    # upper state constraint
    if !isnan.(n.ocp.XU[st])
      if n.ocp.XU[st] != false
        for j in 1:n.ocp.state.pts
          setupperbound(x[j,st], n.ocp.XU[st]/n.ocp.XS[st])
        end
      else
        for j in 1:n.ocp.state.pts
          setlowerbound(x[j,st], n.ocp.XU_var[st,j]/n.ocp.XS[st])
        end
      end
    end
  end

  # control variables
  @variable(n.ocp.mdl,u[1:n.ocp.control.pts,1:n.ocp.control.num])
  n.r.ocp.u = Matrix{Any}(undef, n.ocp.control.pts,n.ocp.control.num)
  for ctr in 1:n.ocp.control.num
      n.r.ocp.u[:,ctr] = @NLexpression(n.ocp.mdl, [j=1:n.ocp.control.pts], n.ocp.CS[ctr]*u[j,ctr])
  end
#n.r.ocp.u=u;

  for ctr in 1:n.ocp.control.num
    if !isnan.(n.ocp.CL[ctr])
      for j in 1:n.ocp.control.pts
        setlowerbound(u[j,ctr], n.ocp.CL[ctr]/n.ocp.CS[ctr])
      end
    end
    if !isnan.(n.ocp.CU[ctr])
      for j in 1:n.ocp.control.pts
        setupperbound(u[j,ctr], n.ocp.CU[ctr]/n.ocp.CS[ctr])
      end
    end
  end

  # boundary constraints
  if n.s.ocp.x0slackVariables   # create handles for constraining the entire initial state
    n.r.ocp.x0Con = Array{Any}(undef, n.ocp.state.num,2) # this is so they can be easily reference when doing MPC
  else
    n.r.ocp.x0Con = []
  end

  if n.s.ocp.xFslackVariables # create handles for constraining the entire final state
    n.r.ocp.xfCon = Array{Any}(undef, n.ocp.state.num,2) # this is so they can be easily reference when doing MPC
  else
    n.r.ocp.xfCon = []
  end

  if n.s.ocp.x0slackVariables # currently adding slack variables for all initial states even if there are no constraints on them
    @variable(n.ocp.mdl,x0s[1:n.ocp.state.num])
    n.ocp.x0s = x0s
    n.r.ocp.x0sCon = Array{Any}(undef, n.ocp.state.num)
  end
  if n.s.ocp.xFslackVariables # currently adding slack variables for all final states even if there are no constraints on them
    @variable(n.ocp.mdl,xFs[1:n.ocp.state.num])
    n.ocp.xFs = xFs
    n.r.ocp.xfsCon = Array{Any}(undef, n.ocp.state.num)
  end

  for st in 1:n.ocp.state.num
    if !isnan(n.ocp.X0[st]) # could have a bool for this
      if n.s.ocp.x0slackVariables
        n.r.ocp.x0Con[st,1] = @constraint(n.ocp.mdl, x[1,st]*n.ocp.XS[st] -n.ocp.x0s[st] <=  n.ocp.X0[st])
        n.r.ocp.x0Con[st,2] = @constraint(n.ocp.mdl,-x[1,st]*n.ocp.XS[st] -n.ocp.x0s[st] <= -n.ocp.X0[st])
        if !isnan(n.ocp.X0_tol[st])
          n.r.ocp.x0sCon[st] = @constraint(n.ocp.mdl, n.ocp.x0s[st] <= n.ocp.X0_tol[st])
        end
      else
        n.r.ocp.x0Con = [n.r.ocp.x0Con; @constraint(n.ocp.mdl, x[1,st]*n.ocp.XS[st] == n.ocp.X0[st])]
      end
    end
    if !isnan(n.ocp.XF[st])
      if n.s.ocp.xFslackVariables
        n.r.ocp.xfCon[st,1] = @constraint(n.ocp.mdl, x[end,st]*n.ocp.XS[st] -n.ocp.xFs[st] <=  n.ocp.XF[st] )
        n.r.ocp.xfCon[st,2] = @constraint(n.ocp.mdl,-x[end,st]*n.ocp.XS[st] -n.ocp.xFs[st] <= -n.ocp.XF[st] )
        if !isnan(n.ocp.XF_tol[st])
          n.r.ocp.xfsCon[st] = @constraint(n.ocp.mdl, n.ocp.xFs[st] <= n.ocp.XF_tol[st])
        end
      else
        n.r.ocp.xfCon = [n.r.ocp.xfCon; @constraint(n.ocp.mdl, x[end,st]*n.ocp.XS[st] == n.ocp.XF[st])]
      end
    end
  end

  @NLparameter(n.ocp.mdl,t0_param == 0.0);   # for now we just start at zero
  n.ocp.t0 = t0_param  # NOTE consider making a constraint that t0 < tf

  if n.s.ocp.integrationMethod == :ps
    n.r.ocp.dynCon = [Array{Any}(undef, n.ocp.Nck[int],n.ocp.state.num) for int in 1:n.ocp.Ni]
    dynamics_expr = [Array{Any}(undef, n.ocp.Nck[int],n.ocp.state.num) for int in 1:n.ocp.Ni]

    if n.s.ocp.finalTimeDV
      @variable(n.ocp.mdl, n.s.ocp.tfMin <= tf <=  n.s.ocp.tfMax)
      n.ocp.tf = tf
      create_tV!(n)          # make a time vector
    end

    for int in 1:n.ocp.Ni
      x_int,u_int = intervals(n,int,n.r.ocp.x,n.r.ocp.u)

      # dynamics
      L = size(x_int)[1]-1;
      dx = Matrix{Any}(undef, L,n.ocp.state.num)
      for st in 1:n.ocp.state.num
        dx[:,st] = DiffEq(n,x_int,u_int,L,st)
      end

      for st in 1:n.ocp.state.num # TODO consider multiplying X*D to reduce computations (i.e. remove this for loop for the states)
        if n.s.ocp.integrationScheme == :lgrExplicit
          dynamics_expr[int][:,st] = @NLexpression(n.ocp.mdl, [j in 1:n.ocp.Nck[int]], sum(n.ocp.DMatrix[int][j,i]*x_int[i,st] for i in 1:n.ocp.Nck[int]+1) - ((n.ocp.tf)/2)*dx[j,st]  )
        elseif n.s.ocp.integrationScheme==:lgrImplicit
          dynamics_expr[int][:,st] = @NLexpression(n.ocp.mdl, [j in 1:n.ocp.Nck[int]], x_int[j+1,st] - x_int[1,st] - ((n.ocp.tf)/2)*sum(n.ocp.IMatrix[int][j,i]*dx[i,st] for i in 1:n.ocp.Nck[int]) )
        end
        for j in 1:n.ocp.Nck[int]
          n.r.ocp.dynCon[int][j,st] = @NLconstraint(n.ocp.mdl, 0. == dynamics_expr[int][j,st])
        end
      end

      # additional constraints
      for num in 1:length(n.ocp.NLcon)
         ch = addCon(n,x_int,u_int,L,num)
         newConstraint!(n,ch,Symbol(string("ch",num))) # TODO could let the user name these
      end
    end
  elseif n.s.ocp.integrationMethod == :tm
    n.r.ocp.dynCon = Array{Any}(undef, n.ocp.N,n.ocp.state.num)
    if n.s.ocp.finalTimeDV
     @variable(n.ocp.mdl, n.s.ocp.tfMin <= tf <= n.s.ocp.tfMax)
     n.ocp.tf = tf
     create_tV!(n)          # make a time vector
    end
    n.ocp.dt = n.ocp.tf/n.ocp.N*ones(n.ocp.N,)

    L = size(n.r.ocp.x)[1]
    dx = Array{Any}(undef, L,n.ocp.state.num)
    for st in 1:n.ocp.state.num
      dx[:,st] = DiffEq(n,n.r.ocp.x,n.r.ocp.u,L,st)
    end

    if n.s.ocp.integrationScheme==:bkwEuler
      for st in 1:n.ocp.state.num
        n.r.ocp.dynCon[:,st] = @NLconstraint(n.ocp.mdl, [j in 1:n.ocp.N], n.r.ocp.x[j+1,st] - n.r.ocp.x[j,st] ==  dx[j+1,st]*n.ocp.tf/(n.ocp.N) );
      end
    elseif n.s.ocp.integrationScheme==:trapezoidal
      for st in 1:n.ocp.state.num
        n.r.ocp.dynCon[:,st] = @NLconstraint(n.ocp.mdl, [j in 1:n.ocp.N], n.r.ocp.x[j+1,st] - n.r.ocp.x[j,st] == 0.5*(dx[j,st] + dx[j+1,st])*n.ocp.tf/(n.ocp.N) )
      end
    end

    # additional constraints
    for num in 1:length(n.ocp.NLcon)
       ch = addCon(n,n.r.ocp.x,n.r.ocp.u,L,num)
       newConstraint!(n,ch,Symbol(string("ch",num))) # TODO could let the user name these
    end
  end

  # save constraint data
  newConstraint!(n,n.r.ocp.x0Con,:x0_con)
  newConstraint!(n,n.r.ocp.xfCon,:xf_con)
  newConstraint!(n,n.r.ocp.dynCon,:dyn_con)

  # save the current working directory for navigation purposes
  n.r.mainDir = pwd()
  return nothing
end

"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/1/2017, Last Modified: 12/06/2019 \n
-------------------------------------------------------------------------------------\n
"""
function configure!(n::NLOpt; kwargs... )
  kw = Dict(kwargs)

  # final time
  if !haskey(kw,:finalTimeDV);n.s.ocp.finalTimeDV=false;
  else; n.s.ocp.finalTimeDV=get(kw,:finalTimeDV,0);
  end

  if !haskey(kw,:tf) && !n.s.ocp.finalTimeDV
    error("\n If the final is not a design variable pass it as: (:tf=>Float64(some #)) \n
        If the final time is a design variable, indicate that as: (:finalTimeDV=>true)\n")
  elseif haskey(kw,:tf) && !n.s.ocp.finalTimeDV
    n.ocp.tf = get(kw,:tf,0)
  elseif n.s.ocp.finalTimeDV
    n.ocp.tf = Any
  end

  # x0 slack variables
  if !haskey(kw,:x0slackVariables)
    if !any(isnan.(n.ocp.X0_tol))
      n.s.ocp.x0slackVariables = true
    else
      n.s.ocp.x0slackVariables = false # default
    end
  else; n.s.ocp.x0slackVariables=get(kw,:x0slackVariables,0);  # TODO check to see if user is using tolerances and passed false to configure
  end

  # xF slack variables
  if !haskey(kw,:xFslackVariables)
    if !any(isnan.(n.ocp.XF_tol))
      n.s.ocp.xFslackVariables = true
    else
      n.s.ocp.xFslackVariables = false # default
    end
  else; n.s.ocp.xFslackVariables=get(kw,:xFslackVariables,0);
  end

  # integrationScheme
  if !haskey(kw,:integrationScheme); n.s.ocp.integrationScheme=:lgrExplicit; # default
  else; n.s.ocp.integrationScheme=get(kw,:integrationScheme,0);
  end

  if n.s.ocp.integrationScheme==:lgrExplicit ||  n.s.ocp.integrationScheme==:lgrImplicit
    n.s.ocp.integrationMethod = :ps
  elseif n.s.ocp.integrationScheme==:trapezoidal || n.s.ocp.integrationScheme==:bkwEuler
    n.s.ocp.integrationMethod = :tm
  else
    error("the :integrationScheme that you specified is not currently implemeted \n")
  end

  if n.s.ocp.integrationMethod==:ps
    if haskey(kw,:N)
      error(" \n N is not an appropriate kwargs for :ps methods \n")
    end
    if !haskey(kw,:Nck);n.ocp.Nck = [10,10,10,10]; # default
    else; n.ocp.Nck = get(kw,:Nck,0);
    end
    n.ocp.Ni = length(n.ocp.Nck)

    for int in 1:n.ocp.Ni
        if (n.ocp.Nck[int]<0)
            error("\n Nck must be > 0");
        end
    end
    n.ocp.state.pts = sum(n.ocp.Nck) + 1
    n.ocp.control.pts = sum(n.ocp.Nck)
    n.ocp.Nck_full = [0;cumsum(n.ocp.Nck.+1)]
    n.ocp.Nck_cum = [0;cumsum(n.ocp.Nck)]

    # initialize node data
    if n.s.ocp.integrationScheme==:lgrExplicit ||  n.s.ocp.integrationScheme==:lgrImplicit
      taus_and_weights = [gaussradau(n.ocp.Nck[int]) for int in 1:n.ocp.Ni];
    end
    n.ocp.tau = [taus_and_weights[int][1] for int in 1:n.ocp.Ni]
    n.ocp.w = [taus_and_weights[int][2] for int in 1:n.ocp.Ni]
    createIntervals!(n)
    DMatrix!(n)

  elseif n.s.ocp.integrationMethod==:tm
    if haskey(kw,:Nck)
      error(" \n Nck is not appropriate kwargs for :tm methods \n")
    end
    if !haskey(kw,:N);n.ocp.N = 100; # default
    else; n.ocp.N = get(kw,:N,0);
    end
    n.ocp.state.pts = n.ocp.N + 1
    n.ocp.control.pts = n.ocp.N + 1
  end
  n.ocp.mXL = falses(n.ocp.state.num)
  n.ocp.mXU = falses(n.ocp.state.num)
  n.ocp.XL_var = Matrix{Float64}(undef, n.ocp.state.num, n.ocp.state.pts)
  n.ocp.XU_var = Matrix{Float64}(undef, n.ocp.state.num, n.ocp.state.pts)

  # solver settings
  if !haskey(kw,:solverSettings);SS = Dict((:name=>:Ipopt)); # default
  else; SS = get(kw,:solverSettings,0); SS = Dict(SS);
  end
  defineSolver!(n,SS)

  # optimal control problem
  OCPdef!(n)

  if n.s.ocp.evalCostates
     if n.s.ocp.integrationMethod != :ps
       error("costates are only implmented for :ps methods")
     end
     n.s.ocp.evalConstraints = true
  end
  return nothing
end


text:
struct DoData
    data::Vector{Float64}
end

function (f::DoData)()
    return f.data[1]
end

data_holder = DoData([1,2,3,4,5])
data_holder()

function go()
    something = [1,2,3,454]
    function zz()
        return something[1]
    end

    return zz()
end


text:
# Wrap win32 API function mouse_event() from the User32 dll.
function mouse_event_wrapper(dwFlags,dx,dy,dwData,dwExtraInfo)
    ccall((:mouse_event, "User32"),stdcall,Void,(UInt32,UInt32,UInt32,UInt32,UInt),dwFlags,dx,dy,dwData,dwExtraInfo)
end

function click()
    mouse_event_wrapper(0x2,0,0,0,0)
    mouse_event_wrapper(0x4,0,0,0,0)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JLD2, Test, LinearAlgebra, Printf, Random

macro read(fid, sym)
    if !isa(sym, Symbol)
        error("Second input to @read must be a symbol (i.e., a variable)")
    end
    esc(:($sym = read($fid, $(string(sym)))))
end
macro write(fid, sym)
    if !isa(sym, Symbol)
        error("Second input to @write must be a symbol (i.e., a variable)")
    end
    esc(:(write($fid, $(string(sym)), $sym)))
end

# Seed random so that we get the same values every time
Random.seed!(1337)

# Define variables of different types
x = 3.7
A = reshape(1:15, 3, 5)
A3 = reshape(1:30, 3, 5, 2)
A4 = reshape(1:120, 2, 3, 4, 5)
Aarray = Vector{Float64}[[1.2,1.3],[2.2,2.3,2.4]]
basic_types = Any[UInt8(42), UInt16(42), UInt32(42), UInt64(42), UInt128(42),
                  Int8(42), Int16(42), Int32(42), Int64(42), Int128(42),
                  Float16(42), Float32(42), Float64(42)]
str = "Hello"
str_unicode = "pandapanda🐼panda"
str_embedded_null = "there is a null\0 in the middle of this string"
strings = String["It", "was", "a", "dark", "and", "stormy", "night"]
empty_string = ""
empty_string_array = String[]
empty_array_of_strings = String[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
c = Complex{Float32}(3,7)
cint = 1+im  # issue 108
C = reshape(reinterpret(Complex{Float64}, B), (4,))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
try
    global MyStruct
    mutable struct MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
d = Dict([(syms[1],"aardvark"), (syms[2], "banana")])
oidd = IdDict([(syms[1],"aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
T = UInt8
Tarr = DataType[UInt8, UInt16, UInt32, UInt64, UInt128]
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]]  # issue #123
typevar = Array{Int}[[1]]
typevar_lb = (Vector{U} where U<:Integer)[[1]]
typevar_ub = (Vector{U} where U>:Int)[[1]]
typevar_lb_ub = (Vector{U} where Int<:U<:Real)[[1]]
arr_undef = Vector{Any}(undef, 1)
arr_undefs = Matrix{Any}(undef, 2, 2)
ms_undef = MyStruct(0)
# Unexported type:
version_info = Base.GIT_VERSION_INFO
# Immutable type:
rng = 1:5
# Custom BitsType (#99)
primitive type MyBT 64 end
bt = reinterpret(MyBT, Int64(55))
btarray = fill(bt, 5, 7)
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = view([1:5;], 1:5)
# Array of empty tuples (to test tuple type params)
arr_empty_tuple = Tuple{}[]
struct EmptyImmutable end
emptyimmutable = EmptyImmutable()
arr_emptyimmutable = [emptyimmutable]
empty_arr_emptyimmutable = EmptyImmutable[]
mutable struct EmptyType end
emptytype = EmptyType()
arr_emptytype = [emptytype]
empty_arr_emptytype = EmptyImmutable[]
uninitialized_arr_emptytype = Vector{EmptyType}(undef, 1)
struct EmptyII
    x::EmptyImmutable
end
emptyii = EmptyII(EmptyImmutable())
struct EmptyIT
    x::EmptyType
end
emptyit = EmptyIT(EmptyType())
mutable struct EmptyTI
    x::EmptyImmutable
end
emptyti = EmptyTI(EmptyImmutable())
mutable struct EmptyTT
    x::EmptyType
end
emptytt = EmptyTT(EmptyType())
struct EmptyIIOtherField
    x::EmptyImmutable
    y::Float64
end
emptyiiotherfield = EmptyIIOtherField(EmptyImmutable(), 5.0)
struct EmptyIIType
    x::Type{Float64}
    y::EmptyImmutable
end
emptyiitype = EmptyIIType(Float64, EmptyImmutable())

# Unicode type field names (#118)
mutable struct MyUnicodeStruct☺{τ}
    α::τ
    ∂ₓα::τ
    MyUnicodeStruct☺{τ}(α::τ, ∂ₓα::τ) where τ = new(α, ∂ₓα)
end
unicodestruct☺ = MyUnicodeStruct☺{Float64}(1.0, -1.0)
# Arrays of matrices (#131)
array_of_matrices = Matrix{Int}[[1 2; 3 4], [5 6; 7 8]]
# Tuple of arrays and bitstype
tup = (1, 2, [1, 2], [1 2; 3 4], bt)
# Empty tuple
empty_tup = ()
# Non-pointer-free struct
struct MyImmutable{T}
    x::Int
    y::Vector{T}
    z::Bool
end
nonpointerfree_immutable_1 = MyImmutable(1, [1., 2., 3.], false)
nonpointerfree_immutable_2 = MyImmutable(2, Any[3., 4., 5.], true)
struct MyImmutable2
    x::Vector{Int}
    MyImmutable2() = new()
end
nonpointerfree_immutable_3 = MyImmutable2()
struct Vague
    name::String
end
vague = Vague("foo")
# Immutable with a union of BitsTypes
struct BitsUnion
    x::Union{Int64,Float64}
end
bitsunion = BitsUnion(5.0)
# Immutable with a union of Types
struct TypeUnionField
    x::Type{T} where T<:Union{Int64,Float64}
end
typeunionfield = TypeUnionField(Int64)
# Generic union type field
struct GenericUnionField
    x::Union{Vector{Int},Int}
end
genericunionfield = GenericUnionField(1)
# Array references
arr_contained = [1, 2, 3]
arr_ref = typeof(arr_contained)[]
push!(arr_ref, arr_contained, arr_contained)
# Object references
mutable struct ObjRefType
    x::ObjRefType
    y::ObjRefType
    ObjRefType() = new()
    ObjRefType(x, y) = new(x, y)
end
ref1 = ObjRefType()
obj_ref = ObjRefType(ObjRefType(ref1, ref1), ObjRefType(ref1, ref1))
# Immutable that requires padding between elements in array
struct PaddingTest
    x::Int64
    y::Int8
end
padding_test = PaddingTest[PaddingTest(i, i) for i = 1:8]
# Empty arrays of various types and sizes
empty_arr_1 = Int[]
empty_arr_2 = Matrix{Int}(undef, 56, 0)
empty_arr_3 = Any[]
empty_arr_4 = Matrix{Any}(undef, 0, 97)
# Moderately big dataset (which will be mmapped)
bigdata = [1:1000000;]
# BigFloats and BigInts
bigint = big(3)
bigfloat = big(3.2)
bigints = big(3).^(1:100)
bigfloats = big(3.2).^(1:100)
struct BigFloatIntObject
    bigfloat::BigFloat
    bigint::BigInt
end
bigfloatintobj = BigFloatIntObject(big(pi), big(typemax(UInt128))+1)
# None
none = Union{}
nonearr = Vector{Union{}}(undef, 5)
# nothing/Nothing
scalar_nothing = nothing
vector_nothing = Union{Int,Nothing}[1,nothing]

# some data big enough to ensure that compression is used:
Abig = kron(Matrix{Float64}(I, 10, 10), rand(20, 20))
Bbig = Any[i for i=1:3000]
Sbig = "A test string "^1000

# Bitstype type parameters
mutable struct BitsParams{x}; end
bitsparamfloat  = BitsParams{1.0}()
bitsparambool   = BitsParams{true}()
bitsparamsymbol = BitsParams{:x}()
bitsparamint    = BitsParams{1}()
bitsparamuint   = BitsParams{0x01}()
bitsparamint16  = BitsParams{Int16(1)}()

# Tuple of tuples
tuple_of_tuples = (1, 2, (3, 4, [5, 6]), [7, 8])

# Zero-dimensional arrays
zerod = Array{Int}(undef)
zerod[] = 1
zerod_any = Array{Any}(undef)
zerod_any[] = 1.0+1.0im

# Type with None typed field
struct NoneTypedField{T}
    a::Int
    b::T

    NoneTypedField{T}(a) where T = new(a)
end
nonetypedfield = NoneTypedField{Union{}}(47)

# Cyclic object
mutable struct CyclicObject
    x::CyclicObject

    CyclicObject() = new()
    CyclicObject(x) = new(x)
end
cyclicobject = CyclicObject()
cyclicobject.x = cyclicobject

# SimpleVector
simplevec = Core.svec(1, 2, Int64, "foo")
iseq(x::Core.SimpleVector, y::Core.SimpleVector) = collect(x) == collect(y)

# JLD issue #243
# Type that overloads != so that it is not boolean
mutable struct NALikeType; end
Base.:!=(::NALikeType, ::NALikeType) = NALikeType()
Base.:!=(::NALikeType, ::Nothing) = NALikeType()
Base.:!=(::Nothing, ::NALikeType) = NALikeType()
natyperef = Any[NALikeType(), NALikeType()]

# JLD2 issue #31 (lots of strings)
lotsastrings = fill("a", 100000)

iseq(x,y) = isequal(x,y)
function iseq(x::Array{EmptyType}, y::Array{EmptyType})
    size(x) != size(y) && return false
    for i = 1:length(x)
        def = isassigned(x, i)
        def != isassigned(y, i) && return false
        if def
            iseq(x[i], y[i]) || return false
        end
    end
    return true
end
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(x::MyImmutable, y::MyImmutable) = (isequal(x.x, y.x) && isequal(x.y, y.y) && isequal(x.z, y.z))
iseq(x::Union{EmptyTI,EmptyTT,EmptyIT}, y::Union{EmptyTI,EmptyTT,EmptyIT}) = iseq(x.x, y.x)
iseq(x::NoneTypedField{Union{}}, y::NoneTypedField{Union{}}) = x.a === y.a
iseq(c1::Array, c2::Array) = length(c1) == length(c2) && all(p->iseq(p...), zip(c1, c2))
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in fieldnames(Base.Sys.CPUinfo)
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
iseq(x::MyUnicodeStruct☺, y::MyUnicodeStruct☺) = (x.α == y.α && x.∂ₓα == y.∂ₓα)
iseq(x::Array{Union{}}, y::Array{Union{}}) = size(x) == size(y)
iseq(x::BigFloatIntObject, y::BigFloatIntObject) = (x.bigfloat == y.bigfloat && x.bigint == y.bigint)
iseq(x::T, y::T) where {T<:Union{EmptyType,EmptyImmutable,NALikeType}} = true
iseq(x::BitsParams{T}, y::BitsParams{S}) where {T,S} = (T == S)
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                @show e
                Base.show_backtrace(stdout, catch_backtrace())
                error("error reading ", $(string(sym)))
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
            if !iseq(tmp, $sym)
                written = $sym
                error("For ", $(string(sym)), ", read value $tmp does not agree with written value $written")
            end
        end
    end
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @assert a == b
function checkexpr(a::Expr, b::Expr)
    @assert a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if isa(a.args[i], Expr) && a.args[i].head == :line
            i += 1
            continue
        end
        if isa(b.args[j], Expr) && b.args[j].head == :line
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(a.args) && j >= length(b.args)
end

fn = joinpath(mktempdir(), "test.jld")
for ioty in [JLD2.MmapIO, IOStream], compress in [false, true]
    @info("[$fn]: Using $(ioty), $(compress ? "compressed" : "uncompressed")")
    @info("  Write time:")
    fid = jldopen(fn, true, true, true, ioty, compress=compress)
    @time begin
    @write fid x
    @write fid A
    @write fid A3
    @write fid A4
    @write fid Aarray
    @write fid basic_types
    @write fid str
    @write fid str_unicode
    @write fid str_embedded_null
    @write fid strings
    @write fid empty_string
    @write fid empty_string_array
    @write fid empty_array_of_strings
    @write fid tf
    @write fid TF
    @write fid AB
    @write fid t
    @write fid c
    @write fid cint
    @write fid C
    @write fid emptyA
    @write fid emptyB
    @write fid ms
    @write fid msempty
    @write fid sym
    @write fid syms
    @write fid d
    @write fid oidd
    @write fid ex
    @write fid T
    @write fid Tarr
    @write fid char
    @write fid unicode_char
    @write fid α
    @write fid β
    @write fid vv
    @write fid version_info
    @write fid rng
    @write fid typevar
    @write fid typevar_lb
    @write fid typevar_ub
    @write fid typevar_lb_ub
    @write fid arr_undef
    @write fid arr_undefs
    @write fid ms_undef
    @write fid bt
    @write fid btarray
    @write fid sa_asc
    @write fid sa_utf8
    @write fid subarray
    @write fid arr_empty_tuple
    @write fid emptyimmutable
    @write fid arr_emptyimmutable
    @write fid empty_arr_emptyimmutable
    @write fid emptytype
    @write fid arr_emptytype
    @write fid empty_arr_emptytype
    @write fid uninitialized_arr_emptytype
    @write fid emptyii
    @write fid emptyit
    @write fid emptyti
    @write fid emptytt
    @write fid emptyiiotherfield
    @write fid emptyiitype
    @write fid unicodestruct☺
    @write fid array_of_matrices
    @write fid tup
    @write fid empty_tup
    @write fid nonpointerfree_immutable_1
    @write fid nonpointerfree_immutable_2
    @write fid nonpointerfree_immutable_3
    @write fid vague
    @write fid bitsunion
    @write fid typeunionfield
    @write fid genericunionfield
    @write fid arr_ref
    @write fid obj_ref
    @write fid padding_test
    @write fid empty_arr_1
    @write fid empty_arr_2
    @write fid empty_arr_3
    @write fid empty_arr_4
    @write fid bigdata
    @write fid bigfloat
    @write fid bigint
    @write fid bigfloats
    @write fid bigints
    @write fid bigfloatintobj
    @write fid none
    @write fid nonearr
    @write fid scalar_nothing
    @write fid vector_nothing
    @write fid Abig
    @write fid Bbig
    @write fid Sbig
    @write fid bitsparamint16
    @write fid bitsparamfloat
    @write fid bitsparambool
    @write fid bitsparamsymbol
    @write fid bitsparamint
    @write fid bitsparamuint
    @write fid tuple_of_tuples
    @write fid zerod
    @write fid zerod_any
    @write fid nonetypedfield
    @write fid cyclicobject
    @write fid simplevec
    @write fid natyperef
    @write fid lotsastrings
    end
    close(fid)

    @info("  Read time:")
    fidr = jldopen(fn, false, false, false, ioty)
    @time begin
    @check fidr x
    @check fidr A
    @check fidr A3
    @check fidr A4
    @check fidr Aarray
    @check fidr basic_types
    @check fidr str
    @check fidr str_unicode
    @check fidr str_embedded_null
    @check fidr strings
    @check fidr empty_string
    @check fidr empty_string_array
    @check fidr empty_array_of_strings
    @check fidr tf
    @check fidr TF
    @check fidr AB
    @check fidr t
    @check fidr c
    @check fidr cint
    @check fidr C
    @check fidr emptyA
    @check fidr emptyB
    @check fidr ms
    @check fidr msempty
    @check fidr sym
    @check fidr syms
    @check fidr d
    @check fidr oidd
    exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
    checkexpr(ex, exr)
    @check fidr T
    @check fidr Tarr
    @check fidr char
    @check fidr unicode_char
    @check fidr α
    @check fidr β
    @check fidr vv
    @check fidr version_info
    @check fidr rng
    @check fidr typevar
    @check fidr typevar_lb
    @check fidr typevar_ub
    @check fidr typevar_lb_ub

    # Special cases for reading undefs
    global arr_undef = read(fidr, "arr_undef")
    if !isa(arr_undef, Array{Any, 1}) || length(arr_undef) != 1 || isassigned(arr_undef, 1)
        error("For arr_undef, read value does not agree with written value")
    end
    global arr_undefs = read(fidr, "arr_undefs")
    if !isa(arr_undefs, Array{Any, 2}) || length(arr_undefs) != 4 || any(map(i->isassigned(arr_undefs, i), 1:4))
        error("For arr_undefs, read value does not agree with written value")
    end
    global ms_undef = read(fidr, "ms_undef")
    if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
        error("For ms_undef, read value does not agree with written value")
    end

    @check fidr bt
    @check fidr btarray
    @check fidr sa_asc
    @check fidr sa_utf8
    @check fidr subarray
    @check fidr arr_empty_tuple
    @check fidr emptyimmutable
    @check fidr arr_emptyimmutable
    @check fidr empty_arr_emptyimmutable
    @check fidr emptytype
    @check fidr arr_emptytype
    @check fidr empty_arr_emptytype
    @check fidr uninitialized_arr_emptytype
    @check fidr emptyii
    @check fidr emptyit
    @check fidr emptyti
    @check fidr emptytt
    @check fidr emptyiiotherfield
    @check fidr emptyiitype
    @check fidr unicodestruct☺
    @check fidr array_of_matrices
    @check fidr tup
    @check fidr empty_tup
    @check fidr nonpointerfree_immutable_1
    @check fidr nonpointerfree_immutable_2
    @check fidr nonpointerfree_immutable_3
    vaguer = read(fidr, "vague")
    @test typeof(vaguer) == typeof(vague) && vaguer.name == vague.name
    @check fidr bitsunion
    @check fidr typeunionfield
    @check fidr genericunionfield

    arr = read(fidr, "arr_ref")
    @test arr == arr_ref
    @test arr[1] === arr[2]

    obj = read(fidr, "obj_ref")
    @test obj.x.x === obj.x.y == obj.y.x === obj.y.y
    @test obj.x !== obj.y

    @check fidr padding_test
    @check fidr empty_arr_1
    @check fidr empty_arr_2
    @check fidr empty_arr_3
    @check fidr empty_arr_4
    @check fidr bigdata
    @check fidr bigfloat
    @check fidr bigint
    @check fidr bigfloats
    @check fidr bigints
    @check fidr bigfloatintobj
    @check fidr none
    @check fidr nonearr
    @check fidr scalar_nothing
    @check fidr vector_nothing
    @check fidr Abig
    @check fidr Bbig
    @check fidr Sbig
    @check fidr bitsparamfloat
    @check fidr bitsparambool
    @check fidr bitsparamsymbol
    @check fidr bitsparamint
    @check fidr bitsparamuint
    @check fidr tuple_of_tuples
    @check fidr zerod
    @check fidr zerod_any
    @check fidr nonetypedfield

    obj = read(fidr, "cyclicobject")
    @test obj.x === obj

    @check fidr simplevec
    @check fidr natyperef
    @check fidr lotsastrings
    end
    close(fidr)
end


text:
using BinDeps

@BinDeps.setup

deps = [
    libpng = library_dependency("libpng", aliases = ["libpng","libpng-1.5.14","libpng15","libpng12.so.0","libpng12"])
]

@BinDeps.install Dict(:libpng => :libpng)