text:
using LinearAlgebra, ForwardDiff, Distributions
using Convex, SCS, ECOS
using IterativeSolvers

function κ_no(z)
    return max.(0.0, z)
end

function Jκ_no(z)
    x = zero(z)
    for i = 1:length(x)
        if z[i] >= 0.0
            x[i] = 1.0
        end
    end
    return Diagonal(x)
end

function κ_soc(z)
    z1 = z[1:end-1]
    z2 = z[end]

    z_proj = zero(z)

    if norm(z1) <= z2
        z_proj = copy(z)
    elseif norm(z1) <= -z2
        z_proj = zero(z)
    else
        a = 0.5 * (1.0 + z2 / norm(z1))
        z_proj[1:end-1] = a * z1
        z_proj[end] = a * norm(z1)
    end
    return z_proj
end

function Jκ_soc(z)
    z1 = z[1:end-1]
    z2 = z[end]
    m = length(z)

    if norm(z1) <= z2
        return Diagonal(ones(m))
    elseif norm(z1) <= -z2
        return Diagonal(zeros(m))
    else
        D = zeros(m, m)
        for i = 1:m
            if i < m
                D[i, i] = 0.5 + 0.5 * z2 / norm(z1) - 0.5 * z2 * ((z1[i])^2.0) / norm(z1)^3.0
            else
                D[i, i] = 0.5
            end
            for j = 1:m
                if j > i
                    if j < m
                        D[i, j] = -0.5 * z2 * z1[i] * z1[j] / norm(z1)^3.0
                        D[j, i] = -0.5 * z2 * z1[i] * z1[j] / norm(z1)^3.0
                    elseif j == m
                        D[i, j] = 0.5 * z1[i] / norm(z1)
                        D[j, i] = 0.5 * z1[i] / norm(z1)
                    end
                end
            end
        end
        return D
    end
end

# problem setup
p = 10
n = p + 1
Σ = Diagonal(rand(p))
Σ_sqrt = sqrt(Σ)
norm(Σ_sqrt * Σ_sqrt - Σ)

c = [zeros(p); 1.0]
G1 = [2.0 * Σ_sqrt zeros(p); zeros(1, p) -1.0]
h = [zeros(p); 1.0]
q = [zeros(p); 1.0]
z = 1.0
G2 = [ones(p)' 0.0]
G3 = [-ones(p)' 0.0]

A = [-G1; -q'; G2; G3]
b = [h; z; 1.0; -1.0]
m = size(A, 1)

"Convex.jl"
x = Variable(n)
prob = minimize(c' * x)
prob.constraints += norm(G1 * x + h) <= q' * x + z
prob.constraints += G2 * x <= 1.0
prob.constraints += G3 * x <= -1.0

@time solve!(prob, ECOS.Optimizer)

@show prob.status
@show prob.constraints
@show x.value
# @show prob.constraints[1].dual
# @show prob.constraints[2].dual
# A = [-1.0 * Diagonal(ones(n)); G; -G]
# b = [zeros(n); h; -h]
# m = 2 *_m  + p
k = n + m + 1
n + 2
Q = Array([zeros(n, n) A' c;
           -A zeros(m, m) b;
           -c' -b'      0.0])

function Pc(z)
    z_proj = zero(z)
    z_proj[1:n] = z[1:n]
    z_proj[n .+ (1:p+2)] = κ_soc(z[n .+ (1:p+2)])
    z_proj[n + p + 2 + 1] = κ_no(z[n + p + 2 + 1])
    z_proj[n + p + 2 + 2] = κ_no(z[n + p + 2 + 2])
    z_proj[n + m + 1] = max(0.0, z[n + m + 1])

    return z_proj
end

function F(z)
    ũ = z[1:k]
    u = z[k .+ (1:k)]
    v = z[2 * k .+ (1:k)]

    [(I + Q) * ũ - (u + v);
     u - Pc(ũ - v);
     ũ - u]
end

function Ju(z)
    ũ = z[1:k]
    u = z[k .+ (1:k)]
    v = z[2k .+ (1:k)]

    JP = zeros(m, m)
    JP[1:p+2, 1:p+2] = Jκ_soc(ũ[n .+ (1:m)][1:p+2] - v[n .+ (1:m)][1:p+2])
    JP[p+2 .+ (1:1), p+2 .+ (1:1)] = Jκ_no(ũ[n .+ (1:m)][p+2 .+ (1:1)] - v[n .+ (1:m)][p+2 .+ (1:1)])
    JP[p+2 + 1 .+ (1:1), p+2 + 1 .+ (1:1)] = Jκ_no(ũ[p+2 .+ (1:m)][p+2 + 1 .+ (1:1)] - v[p+2 .+ (1:m)][p+2 + 1 .+ (1:1)])

    if z[k] - z[3k] >= 0.0
        ℓ = 1.0
    else
        ℓ = 0.0
    end

    [-I zeros(n, m) zeros(n) I zeros(n, m) zeros(n) I zeros(n, m) zeros(n);
     zeros(m, n) -JP zeros(m) zeros(m, n) I zeros(m) zeros(m, n) JP zeros(m);
     [zeros(1, n) zeros(1, m) -ℓ zeros(1, n) zeros(1, m) 1.0 zeros(1, n) zeros(1, m) ℓ]]
end

function J(z)
    [(I + Q) -I -I
     Ju(z);
     I -I zeros(k, k)]
end

ũ = zeros(k)
u = zeros(k)
v = zeros(k)
ũ[end] = 1.0
u[end] = 1.0
v[end] = 1.0

z = [ũ; u; v]
z = rand(3k)
F(z)
rank(Ju(z))
J(z)
rank(ForwardDiff.jacobian(F, z))
norm(J(z) - ForwardDiff.jacobian(F, z), Inf)
gmres(J(z), -F(z))

function solve()
    ũ = zeros(k)
    u = zeros(k)
    v = zeros(k)
    ũ[end] = 1.0
    u[end] = 1.0
    v[end] = 1.0

    z = [ũ; u; v]

    extra_iters = 0

    for i = 1:100
        _F = F(z)
        _J = J(z)
        Δ = gmres(_J, 1.0 * _F, abstol = 1.0e-12, maxiter = i + extra_iters)

        iter = 0
        α = 1.0
        while norm(F(z - α * Δ))^2.0 >= (1.0 - 0.001 * α) * norm(_F)^2.0 && α > 1.0e-4
            α = 0.5 * α
            # println("   α = $α")
            iter += 1
            if iter > 100
                @error "line search fail"
                # return z
                x = z[k .+ (1:n)]
                y = z[k + n .+ (1:m)]
                τ = z[k + k]
                κ = z[3k]
                println("τ = $τ")
                println("κ = $κ")
                return x ./ τ, z
            end
        end

        if α <= 1.0e-4
            extra_iters += 1
        end

        println("iter ($i) -> norm: $(norm(F(z)))")

        z .-= α * Δ
    end

    # return z
    x = z[k .+ (1:n)]
    y = z[k + n .+ (1:m)]
    τ = z[k + k]
    κ = z[3k]
    println("τ = $τ")
    println("κ = $κ")
    return x ./ τ, z
end

x_sol, z_sol = solve()
@show norm(x_sol - x.value)

y_sol = z_sol[1:k] - z_sol[2k .+ (1:k)]


text:
JULIA_ENV = joinpath(@__DIR__, "..", "..")
data_root = joinpath(@__DIR__, "..", "..", "data")
data_input_root = joinpath(data_root, "input", "processed")
data_output_root = joinpath(data_root, "output")

worker_list = [("localhost", 1)]
exeflags = `--project="$JULIA_ENV"`
sshflags= `-i path/to/remote/ssh/key`

fmt_string = "[{name} | {date} | {level}]: {msg}"
loglevel = "debug"

ext: jl

max_stars_repo_name: englhardt/des-evaluation

max_issues_repo_name: englhardt/des-evaluation

max_forks_repo_name: englhardt/des-evaluation


text:
module FiniteStateMachines

@inline signalstate(state, prob, signal) = map(tbx -> state(tbx, signal, prob), toolboxes(prob))

"""
    runstate(state, prob)

Run code associated with the state provided. The default fallback is to call
state as a function with prob as an argument. Add specialised methods as desired
to change the defaults.
"""
@noinline function runstate(state, prob)
	signalstate(state, prob, :before)
	result = state(prob)
	signalstate(state, prob, :after)
	return result
end

"""
    runfsm(initialstate, prob)

Run the finite state machine starting with the initialstate.
"""
function runfsm(initialstate, prob)
    state = initialstate
    run = true
    while run
        (run, state) = runstate(state, prob)
    end
end

end # end module


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------
@testset ExtendedTestSet "circuit construction" begin
    N = 3
    cgs = [
        circuit_gate(1, X),
        circuit_gate(2, HadamardGate()),
        circuit_gate(3, Z),
        circuit_gate(1, Y),
        ]
    ref_moments = Moment[
        Moment([circuit_gate(1, X),
                circuit_gate(2, HadamardGate()),
                circuit_gate(3, Z)]),
        Moment([circuit_gate(1, Y)])
    ]
    iwires = [(1,), (2,), (3,)]

    @testset "basic circuit with circuit gate" begin
        c = Circuit{N}(cgs[1])

        @test c[1] ≈ Moment(cgs[1])

        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        @test distribution(ψ, c) ≈ apply(ψ, c.moments)

        @test_throws ErrorException("Circuit does not contain any measurement operators") apply(ψ, c)

        meas = mop.([X, X, X], iwires)

        c = Circuit{N}(cgs[1], meas)
        ψs = apply(ψ, c.moments)

        @test c[1] ≈ Moment(cgs[1])
        @test [dot(ψs, apply(ψs, m)) for m in meas] ≈ apply(ψ, c)
        @test [dot(ψs, apply(ψs, m)) for m in meas] ≈ apply!(s, c)
    end

    @testset "basic circuit construction" begin
        c = Circuit{N}(cgs)

        for i in 1:length(c)
            @test c[i] ≈ ref_moments[i]
        end

        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        @test distribution(ψ, c) ≈ apply(ψ, c.moments)

        @test_throws ErrorException("Circuit does not contain any measurement operators") apply(ψ, c)
        @test_throws ErrorException("Circuit does not contain any measurement operators") apply!(s, c)
    end

    @testset "circuit construction without gates" begin
        meas = mop.([X, X, X], iwires)

        c = Circuit{N}()

        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)

        @test_throws ErrorException("Circuit does not contain any gates") apply(ψ, c)
        @test_throws ErrorException("Circuit does not contain any gates") apply!(s, c)

        c = Circuit{N}(meas)

        @test_throws ErrorException("Circuit does not contain any gates") apply(ψ, c)
        @test_throws ErrorException("Circuit does not contain any gates") apply!(s, c)
    end

    @testset "circuit construction with measurements" begin

        meas = mop.([X, X, X], iwires)
        c = Circuit{N}(cgs, meas)
        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)

        @test distribution(ψ, c) ≈ apply(ψ, c.moments)

        ψs = apply(ψ, cgs)
        @test [dot(ψs, apply(ψs, m)) for m in meas] ≈ apply(ψ, c)
        @test [dot(ψs, apply(ψs, m)) for m in meas] ≈ apply!(s, c)
    end
end

@testset ExtendedTestSet "circuit helper functions" begin
    N = 3
    cgs = [
        circuit_gate(1, X),
        circuit_gate(2, HadamardGate()),
        circuit_gate(3, Z),
        circuit_gate(1, Y),
        ]
    ref_moments = Moment[
        Moment([circuit_gate(1, X),
                circuit_gate(2, HadamardGate()),
                circuit_gate(3, Z)]),
        Moment([circuit_gate(1, Y)])
    ]
    @testset "circuit add measurement operator" begin
        c = Circuit{N}(cgs)

        meas = mop(X, 1)
        add_measurement!(c, meas)

        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        ψs = apply(ψ, cgs)

        @test [dot(ψs, apply(ψs, meas))] ≈ apply(ψ, c)
        @test [dot(ψs, apply(ψs, meas))] ≈ apply!(s, c)
    end

    @testset "circuit add measurement operators" begin
        c = Circuit{N}(cgs)

        meas = mop.([X, X, X], (1,2,3))
        add_measurement!(c, meas)

        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        ψs = apply(ψ, cgs)

        @test [dot(ψs, apply(ψs,m)) for m in meas] ≈ apply(ψ, c)
        @test [dot(ψs, apply(ψs,m)) for m in meas] ≈ apply!(s, c)
    end

    @testset "circuit add measurement operator exceptions" begin
        meas = mop.([X, X, X], (1,2,3))
        c = Circuit{N}(cgs, meas)

        meas2 = mop(Y, 1)
        @test_throws ErrorException("Measurement operators do not commute") add_measurement!(c, meas2)
        @test_throws ErrorException("Measurement operators do not commute") add_measurement!(c, [meas2])

        meas3 = mop(X, 5)
        @test_throws ErrorException("Measurement operators affecting 5 wires provided for circuit of size 3") add_measurement!(c, meas3)
        @test_throws ErrorException("Measurement operators affecting 5 wires provided for circuit of size 3") add_measurement!(c, [meas3])
    end

    @testset "circuit array operations" begin
        m1 = Moment([circuit_gate(3, ZGate(), 1), circuit_gate(2, SGate())])
        m2 = Moment([circuit_gate(2, ZGate()), circuit_gate(1, TdagGate())])
        cgc = [circuit_gate(3, ZGate(), 1), circuit_gate(2, SGate()), circuit_gate(2, ZGate()), circuit_gate(1, TdagGate())]
        c = Circuit{3}(cgc)

        @test all(c[:] .≈ c.moments)
        @test all(c[1:2] .≈ [m1, m2])
        @test c[1:2] == c.moments
        @test isempty(c) == false
        @test m2 ≈ pop!(c)
        @test isempty(c) == false
        @test m1 ≈ pop!(c)
        @test isempty(c) == true
        @test_throws ArgumentError("array must be non-empty") pop!(c)
    end

    @testset "circuit reverse" begin
        meas = mop.([X, X, X], (1,2,3))
        c = Circuit{N}(cgs, meas)
        reverse_c = reverse(c)
        reverse_ref_moments = reverse(reverse.(ref_moments))
        for i in 1:length(c)
            @test reverse_c[i] ≈ reverse_ref_moments[i]
            @test c[i] ≈ ref_moments[i]
        end

        reverse_c = reverse!(c)
        for i in 1:length(c)
            @test reverse_c[i] ≈ reverse_ref_moments[i]
            @test c[i] ≈ reverse_ref_moments[i]
        end
    end
end


text:
# fake a package
module QMC

__precompile__(false)

include("dqmc_framework.jl")

function mc_from_inxml(inxml::AbstractString)
    Random.seed!(1234)
    if SUPPRESS_DQMC_INIT_OUTPUT
        old_stdout = stdout
        (rp, wp) = redirect_stdout()
    end

    @assert isfile(inxml)
    p = Params()
    p.output_file = "$(inxml[1:end-7]).out.h5.running"
    xml2parameters!(p, inxml)
    mc = DQMC(p)
    init!(mc)

    if SUPPRESS_DQMC_INIT_OUTPUT
        redirect_stdout(old_stdout)
        close(wp); close(rp);
    end
    return mc
end

function mc_from_inxml(inxml::AbstractString, initconf)
    mc = mc_from_inxml(inxml)
    mc.p.hsfield .= initconf
    return mc
end

if !haskey(ENV, "JULIA_DQMC")
    ENV["JULIA_DQMC"] = dirname(dirname(@__FILE__)) # parent directory
end

(@isdefined SUPPRESS_DQMC_INIT_OUTPUT) || (const SUPPRESS_DQMC_INIT_OUTPUT = false)

# exportall:
for n in names(@__MODULE__, all=true)
    if Base.isidentifier(n) && n ∉ (Symbol(@__MODULE__), :eval)
        @eval export $n
    end
end

end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
"""
    _get_authors_from_CSV(csv::String) -> Vector{String}

Extract the values of the `Author` column from the CSV file.

# Arguments
- `csv::String`: $(CSV_ARGUMENT)

# Returns
- `Vector{String}`: the authors

# Throws
- Exceptions from: [`get_all`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

_file, io = mktemp()
println(io, HEADER, '\\n', ",,Susanna Kaysen,,1")
flush(io)
file = _file * ".csv"
mv(_file, file)

PDFHighlights.Internal.CSV._get_authors_from_CSV(file) == ["Susanna Kaysen"]

# output

true
```
"""
function _get_authors_from_CSV(csv::String)::Vector{String}
    _, _, authors, _, _ = get_all(csv)
    return authors
end


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

@testset "Distances & Amplitude" begin
  @test distance() == # Euclidean distance
  @test bhattDist() == # Bhattacharyya distance
  @test amplitude() == # amplitude
end

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

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
"""
# Description

A logical node that represents a limit operation. It has two fields:

1. `source`: A node describing a relation that will be limited.
2. `n`: The number of rows to limit to.
"""
struct Limit <: LogicalNode
    source::Any
    n::Int
end

function Base.print(io::IO, node::Limit)
    println(io, "Limit LogicalNode")
    println(io, type(node.source))
    println(io, node.n)
end


text:
const DICTMSG = OrderedDict{Int64, String}
import Base: push!

export RateLimit, isexceeded, numleft

struct RateLimit
    buffer::CircularBuffer{DateTime}
    period::Dates.Second
    function RateLimit(maxMessages::Int64, period::Dates.Second)
        return new(CircularBuffer{DateTime}(maxMessages), period)
    end
end

function push!(this::RateLimit, n::DateTime)
    push!(this.buffer, n)
end

function isexceeded(this::RateLimit, n::DateTime)
    return n - first(this) < this.period
end

"messages per second"
function (this::RateLimit)()::Float64
    dur = (last(this.buffer) - first(this.buffer)).value/1000.0
    nm = length(this.buffer)
    return dur ≈ 0.0 ? nm / this.period.value : nm / dur
end

function numleft(this::RateLimit, n::DateTime)::Int64
    out = capacity(this.buffer) - length(this.buffer)
    for t in this.buffer
        if n - t > this.period
            out += 1
        end
    end
    return out
end

struct FIXIncomingMessages
    login::DICTMSG
    logout::DICTMSG
    executionReports::Dict{String, Vector{DICTMSG}} #type -> reports
    orderCancelReject::Vector{DICTMSG}
    orderReject::Vector{DICTMSG}
    resend::Vector{DICTMSG}
    heartbeat::Vector{DICTMSG}
    msgseqnum::Container{String}
    ratelimit::RateLimit
    function FIXIncomingMessages(ratelimit::RateLimit)
        return new(DICTMSG(),
                DICTMSG(),
                Dict{String, Vector{DICTMSG}}(),
                Vector{DICTMSG}(0),
                Vector{DICTMSG}(0),
                Vector{DICTMSG}(0),
                Vector{DICTMSG}(0),
                Container("0"),
                ratelimit)
    end
end

struct FIXOutgoingMessages
    login::DICTMSG
    logout::DICTMSG
    newOrderSingle::Dict{String, DICTMSG} #client order id ->
    orderCancelRequest::Dict{String, DICTMSG} #client order id ->
    orderStatusRequest::Vector{DICTMSG}
    orderCancelReplaceRequest::Dict{String ,DICTMSG}
    msgseqnum::Container{Int64}
    ratelimit::RateLimit
    function FIXOutgoingMessages(ratelimit::RateLimit)
        return new(DICTMSG(),
                DICTMSG(),
                Dict{String, DICTMSG}(),
                Dict{String, DICTMSG}(),
                Vector{DICTMSG}(0),
                Dict{String ,DICTMSG}(),
                Container(0),
                ratelimit)
    end
end

const ClientOrderID = String
const ExchangeOrderID = String

struct FIXOpenOrders
    client_id::Dict{String, DICTMSG}
    exchng_id::Dict{String, DICTMSG}
    verbose::Bool
    function FIXOpenOrders(verbose::Bool = false)
        return new(Dict{String, DICTMSG}(),
                    Dict{String, DICTMSG}(),
                    verbose)
    end
end

struct FIXMessageManagement
    incoming::FIXIncomingMessages
    outgoing::FIXOutgoingMessages
    open::FIXOpenOrders
    function FIXMessageManagement(ratelimit::RateLimit, verbose::Bool = false)
        return new(FIXIncomingMessages(ratelimit),
                    FIXOutgoingMessages(ratelimit),
                    FIXOpenOrders(verbose))

    end
end

function getPosition(this::FIXIncomingMessages, instrument::String)
    out = 0.0
    if !haskey(this.executionReports, "1")
        return out
    end

    for report in this.executionReports["1"]
        if report[55] ≠ instrument
            continue
        end

        if report[54] == "1"
            out += parse(Float64, report[32])
        else
            out -= parse(Float64, report[32])
        end
    end

    return out
end

function getPositions(this::FIXIncomingMessages)
    out = Dict{String, Float64}()

    if !haskey(this.executionReports, "1")
        return out
    end

    for report in this.executionReports["1"]
        if !haskey(out, report[55])
            out[report[55]] = 0.0
        end

        sz = parse(Float64, report[32])
        if report[54] == "1"
            out[report[55]] = out[report[55]] + sz
        else
            out[report[55]] = out[report[55]] - sz
        end
    end

    return out
end

function onNewMessage(this::FIXIncomingMessages, msg::DICTMSG)
    this.msgseqnum.data = msg[34]
    msg_type = msg[35]

    if msg_type == "8"
        addExecutionReport(this, msg)
    elseif msg_type == "9"
        addOrderCancelReject(this, msg)
    elseif msg_type == "3"
        @printf("[%ls][FIX:ER] received order reject. Reason: %ls\n", now(), msg[58])
        addOrderReject(this, msg)
    elseif msg_type == "0"
        addHeartbeat(this, msg)
    elseif msg_type == "A"
        addLogin(this, msg)
    elseif msg_type == "5"
        addLogout(this, msg)
    elseif msg_type == "2"
        addResend(this, msg)
    else
        @printf("[%ls] Unknown INCOMING message type: %ls\n", now(), msg_type)
    end
end

function addLogin(this::FIXIncomingMessages, msg::DICTMSG)
    [this.login[k] = v for (k, v) in msg]
end

function addLogout(this::FIXIncomingMessages, msg::DICTMSG)
    [this.logout[k] = v for (k, v) in msg]
end

function addExecutionReport(this::FIXIncomingMessages, msg::DICTMSG)
    exec_type = msg[150]
    if !haskey(this.executionReports, exec_type)
        this.executionReports[exec_type] = DICTMSG[]
    end

    push!(this.executionReports[msg[150]], msg)
end

function addOrderCancelReject(this::FIXIncomingMessages, msg::DICTMSG)
    push!(this.orderCancelReject, msg)
end

function addOrderReject(this::FIXIncomingMessages, msg::DICTMSG)
    push!(this.orderReject, msg)
end

function addHeartbeat(this::FIXIncomingMessages, msg::DICTMSG)
    push!(this.heartbeat, msg)
end

function addResend(this::FIXIncomingMessages, msg::DICTMSG)
    push!(this.resend, msg)
end

function onNewMessage(this::FIXOutgoingMessages, msg::DICTMSG)
    msg_type = msg[35]
    this.msgseqnum.data = parse(Int64, msg[34])

    if msg_type == "D"
        addNewOrderSingle(this, msg)
    elseif msg_type == "F"
        addOrderCancelRequest(this, msg)
    elseif msg_type == "H"
        addOrderStatusRequest(this, msg)
    elseif msg_type == "A"
        addLogin(this, msg)
    elseif msg_type == "5"
        addLogout(this, msg)
    elseif msg_type == "G"
        addCancelReplace(this, msg)
    else
        @printf("[%ls] Unknown OUTGOING msg type %ls\n", now(), msg_type)
    end
    return nothing
end

function addLogin(this::FIXOutgoingMessages, msg::DICTMSG)
    [this.login[k] = v for (k, v) in msg]
end

function addLogout(this::FIXOutgoingMessages, msg::DICTMSG)
    [this.logout[k] = v for (k, v) in msg]
end

function addNewOrderSingle(this::FIXOutgoingMessages, msg::DICTMSG)
    this.newOrderSingle[msg[11]] = msg
end

function addOrderCancelRequest(this::FIXOutgoingMessages, msg::DICTMSG)
    this.orderCancelRequest[msg[11]] = msg
end

function addOrderStatusRequest(this::FIXOutgoingMessages, msg::DICTMSG)
    push!(this.orderStatusRequest, msg)
end

function addCancelReplace(this::FIXOutgoingMessages, msg::DICTMSG)
    this.orderCancelReplaceRequest[msg[1]]
end

function onSent(this::FIXMessageManagement, msg::DICTMSG)
    onNewMessage(this.outgoing, msg)
    push!(this.outgoing.ratelimit, now())
    return nothing
end

function onGet(this::FIXMessageManagement, msg::DICTMSG)
    onNewMessage(this.incoming, msg)

    if msg[35] == "8"
        onExecutionReport(this.open, msg)
    end

    return nothing
end

function onCancelReject(this::FIXMessageManagement, msg::DICTMSG)
    addOrderCancelReject(this.incoming, msg)
end

function onOrderReject(this::FIXMessageManagement, msg::DICTMSG)
    addOrderReject(this.incoming, msg)
end

function onNewOrder(this::FIXOpenOrders, msg::DICTMSG)
    this.client_id[msg[11]] = msg
    this.exchng_id[msg[37]] = msg
    return nothing
end

function onFill(this::FIXOpenOrders, msg::DICTMSG)
    exchng_id = msg[37]
    filled = msg[32]
    orig_exec = getOrderByExchangeID(this, exchng_id)
    prev_amount = parse(Float64, orig_exec[38])
    filled_amount = parse(Float64, filled)

    cur_amount = prev_amount - filled_amount
    orig_exec[38] = string(cur_amount)
end

function onDone(this::FIXOpenOrders, msg::DICTMSG)
    deleteOrderByExchangeID(this, msg[37])
end

function onCancel(this::FIXOpenOrders, msg::DICTMSG)
    deleteOrderByExchangeID(this, msg[37])
end

function onReject(this::FIXOpenOrders, msg::DICTMSG)
    deleteOrderByExchangeID(this, msg[37])
end

function getOrderByExchangeID(this::FIXOpenOrders, exchange_id::String)
    return this.exchng_id[exchange_id]
end

function hasOrderByClientID(this::FIXOpenOrders, client_id::String)
    return haskey(this.client_id, client_id)
end

function hasOrderByExchangeID(this::FIXOpenOrders, exchange_id::String)
    return haskey(this.exchng_id, exchange_id)
end

function deleteOrderByExchangeID(this::FIXOpenOrders, exchange_id::String)::Void
    if !hasOrderByExchangeID(this, exchange_id)
        return nothing
    end

    orig_exec = getOrderByExchangeID(this, exchange_id)

    client_id = orig_exec[11]
    exchng_id = orig_exec[37]

    delete!(this.client_id, client_id)
    delete!(this.exchng_id, exchng_id)

    return nothing
end

function onExecutionReport(this::FIXOpenOrders, msg::DICTMSG)
    exec_type = msg[150]
    if exec_type == "0"
        this.verbose && @printf("[%ls][FIX:ER] new order\n", now())
        onNewOrder(this, msg)
    elseif exec_type == "1" #fill
        this.verbose && @printf("[%ls][FIX:ER] fill\n", now())
        onFill(this, msg)
    elseif exec_type == "3"
        this.verbose && @printf("[%ls][FIX:ER] done\n", now())
        onDone(this, msg)
    elseif exec_type == "4"
        this.verbose && @printf("[%ls][FIX:ER] cancel\n", now())
        onCancel(this, msg)
    elseif exec_type == "8"
        onReject(this, msg)
        this.verbose && @printf("[%ls][FIX:ER] reject\n", now())
    else
        @printf("Received unhandled execution report of type %ls\n", exec_type)
    end
end

function getOrders(this::FIXOpenOrders)::Vector{DICTMSG}
    return collect(values(this.client_id))
end

function getOpenOrders(this::FIXMessageManagement)::Vector{DICTMSG}
    return getOrders(this.open)
end

function getNextOutgoingMsgSeqNum(this::FIXMessageManagement)::String
    return getNextOutgoingMsgSeqNum(this.outgoing)
end

function getNextOutgoingMsgSeqNum(this::FIXOutgoingMessages)::String
    return string(this.msgseqnum.data + 1)
end

function getPosition(this::FIXMessageManagement, instrument::String)
    return getPosition(this.incoming, instrument)
end

function getPositions(this::FIXMessageManagement)
    return getPositions(this.incoming)
end


text:
# =========================================================================== #
# Compliant julia 1.x

# Using the following packages
using JuMP, GLPK
using LinearAlgebra

include("loadSPP.jl")
include("setSPP.jl")
include("getfname.jl")
include("experiment.jl")
#include("DM1_1.jl")
#include("DM1_2.jl")
#include("DM2_1.jl")
include("DM2_2.jl")
include("recuit2.jl")
#include("brouillon.jl")
include("noVerbose.jl")
# =========================================================================== #

# Setting the data
#fname = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data/didactic.dat"  # path for a standard config on macOS
#cost, matrix, n, m = loadSPP(fname)

# Proceeding to the optimization
#solverSelected = GLPK.Optimizer
#ip, ip_x = setSPP(solverSelected, cost, matrix)
#println("Solving..."); optimize!(ip)

# Displaying the results
#println("z  = ", objective_value(ip))
#sprint("x  = "); println(value.(ip_x))

# =========================================================================== #

# Collecting the names of instances to solve C:Users/Documents/ORO/Metaheuristiques/DM1_Metaheuristique/
#target = "/comptes/E197494S/DM/Metaheuristique/DM1_Metaheuristique-master/solveSPP-master/Data"            # path for a standard config on Linux
target = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data"
#target = "C:/Users/Emmanuel/Documents/ORO/Métaheurisques/DM1_Metaheuristique/solveSPP-master/Data"                         # path for personal config on windows10
cd
fnames = getfname(target)

#===================================================#

for f in fnames
    println("====================================================================")
    println(f)
    cost, matrix, n, m = loadSPP(f)
    ite = 15
    alphaset = 0
    temps = 2
    nbAlpha = 5

    #for i in 1:5
    #    @time (SOL,z) =Glouton(cost, matrix, n, m)
    #    @time (SOL,z) =GRASP(cost, matrix, n, m, 0.9)
    #    @time (SOL,z) =Glouton(cost, matrix, n, m)
    #    @time proba = ReactiveGRASP(matrix,cost, n, m, 15, 5, 0, 2)

    #end

    zinit = zeros(Int64, ite) # zero
    zls   = zeros(Int64, ite) # zero
    zbest = zeros(Int64, ite) # zero

    # calcule la valeur du pas pour les divisions

    println("Experimentation ReactiveGRASP-SPP avec :")
    println("  nbIterationGrasp  = ", ite)

    println(" ")
    cpt = 0

    # run non comptabilise (afin de produire le code compile)
    #zinit, zls, zbest = graspSPP(allfinstance[1], 0.5, 1)

    print("  ",f," : ")
    (liste_zavg,liste_zmax,liste_zmin,z_rouge,z_vert,ligne_verte) = ReactiveGRASP(matrix,cost, n, m, ite, nbAlpha, alphaset, temps)
    gr = GRASP(cost,matrix,n,m,0.8)
    println(" ")

    println("Zavg : ")
    #println(liste_zavg)
    println("Zmin : ")
    #println(liste_zmin)
    println("Zmax : ")
    #println(liste_zmax)

        #Pkg.add("PyPlot") # Mandatory before the first use of this package
    println(" ");println("  Graphiques de synthese")

    #plotRunGrasp(f, liste_zmin, zls, liste_zmax)
    plotAnalyseGrasp(f, 1:length(liste_zavg), liste_zavg, liste_zmin, liste_zmax )
    plotRunGrasp(f, z_rouge, z_vert, ligne_verte)

     #On fait une solution de base
    (Sol,z, crts) = GRASP(cost, matrix, n, m, 0.75)
    #test Recuit
    println(z)
    test = SA(Sol,n,m,cost,crts,matrix,exchange1_2)
    println("solution du SA :", test[2])#=
#    plotAnalyseGrasp(f, x, zmoy[instancenb,:], zmin[instancenb,:], zmax[instancenb,:] )
#    plotCPUt(allfinstance, tmoy)

using Plots
#@userplot PortfolioComposition

     function f(pc::PortfolioComposition)
        weights, returns = pc.args
        weights = cumsum(weights,dims=2)
        seriestype := :shape
        for c=1:size(weights,2)
            sx = vcat(weights[:,c], c==1 ? zeros(length(returns)) : reverse(weights[:,c-1]))
            sy = vcat(returns, reverse(returns))
            Shape(sx, sy)
        end
    end
    using Random
    tickers = ["0.2", "0.4", "0.6", "0.8","1"]
    N = 10
    D = length(tickers)
    weights = rand(N, D)
    weights ./= sum(weights, dims=2)
    returns = sort!((1:N) + D * randn(N))
    portfoliocomposition(weights, returns, labels=permutedims(tickers))=#
	    (vecZBest, vecZAvg, vecZMin,zBest, zAvg, zMin, t) = expPop("pb_200rnd800.dat")
		x = 1:50
		y = (vecZbest,vecZAvg,vecZMin)
		plot(x,y,label=["zBest,zAvg,zMax"])
end

cd("../")

# ============================================================================ #
#n = 9
#m = 7
#a  = zeros(Float64, n)
#cost = [10, 5, 8, 6, 9, 13, 11, 4, 6]
#matrix = [1 1 1 0 1 0 1 1 0; 0 1 1 0 0 0 0 1 0;0 1 0 0 1 1 0 1 1;0 0 0 1 0 0 0 0 0;1 0 1 0 1 1 0 0 1;0 1 1 0 0 0 1 0 1; 1 0 0 1 1 0 0 1 1]
#a = Glouton(cost, matrix, n, m)
#println("Solution: ",a)


text:
export Sum

struct Sum{K,
           C <: AbstractArray,
           D <: AbstractArray,
           L <:NTuple{K,AbstractOperator}} <: AbstractOperator
  A::L
  bufC::C
  bufD::D
  function Sum(A::L, bufC::C, bufD::D) where {C, D, K, L <: NTuple{K,AbstractOperator}}
    if any([size(a) != size(A[1]) for a in A])
      throw(DimensionMismatch("cannot sum operator of different sizes"))
    end
    if any([codomainType(A[1]) != codomainType(a) for a in A]) ||
      any([  domainType(A[1]) !=   domainType(a) for a in A])
      throw(DomainError(A,"cannot sum operator with different codomains"))
    end
    new{K, C, D, L}(A, bufC, bufD)
  end
end

Sum(A::AbstractOperator) = A

function Sum(A::Vararg{AbstractOperator})
  s = size(A[1],1)
  t = codomainType(A[1])
  bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)

  s = size(A[1],2)
  t = domainType(A[1])
  bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)

  return Sum(A, bufC, bufD)
end

# special cases
Sum(L1::AbstractOperator, L2::Sum{K,C,D}) where {K,C,D} =
Sum((L1,L2.A...), L2.bufC, L2.bufD)

# Mappings

@generated function mul!(y::C, S::Sum{K,C,D}, b::D) where {K,C,D}
  ex = :(mul!(y,S.A[1],b))
  for i = 2:K
    ex = quote
      $ex
      mul!(S.bufC,S.A[$i],b)
    end
    ex = :($ex; y .+= S.bufC)
  end
  ex = quote
    $ex
    return y
  end
end

@generated function mul!(y::D, A::AdjointOperator{Sum{K,C,D,L}}, b::C) where {K,C,D,L}
  ex = :(S = A.A; mul!(y,S.A[1]',b))
  for i = 2:K
    ex = quote
      $ex
      mul!(S.bufD, S.A[$i]', b)
    end
    ex = :($ex; y .+= S.bufD)
  end
  ex = quote
    $ex
    return y
  end
end

# Properties

size(L::Sum) = size(L.A[1])

domainType(S::Sum{K, C, D, L}) where {K,C,D<:AbstractArray,L} =    domainType(S.A[1])
domainType(S::Sum{K, C, D, L}) where {K,C,D<:Tuple        ,L} =   domainType.(Ref(S.A[1]))
codomainType(S::Sum{K, C, D, L}) where {K,C<:AbstractArray,D,L} =  codomainType(S.A[1])
codomainType(S::Sum{K, C, D, L}) where {K,C<:Tuple        ,D,L} = codomainType.(Ref(S.A[1]))

fun_domain(S::Sum)   = fun_domain(S.A[1])
fun_codomain(S::Sum) = fun_codomain(S.A[1])

fun_name(S::Sum) =
length(S.A) == 2 ? fun_name(S.A[1])*"+"*fun_name(S.A[2]) : "Σ"

is_linear(L::Sum)        = all(is_linear.(L.A))
is_null(L::Sum)          = all(is_null.(L.A))
is_diagonal(L::Sum)      = all(is_diagonal.(L.A))
is_full_row_rank(L::Sum) = any(is_full_row_rank.(L.A))
is_full_column_rank(L::Sum) = any(is_full_column_rank.(L.A))

diag(L::Sum) = (+).(diag.(L.A)...,)

# utils
function permute(S::Sum, p::AbstractVector{Int})
  AA = ([permute(A,p) for A in S.A]...,)
  return Sum(AA, S.bufC, ArrayPartition(S.bufD.x[p]...))
end

remove_displacement(S::Sum) = Sum(remove_displacement.(S.A), S.bufC, S.bufD)


text:
include("AngularMomentum.jl")
include("Pseudopotentials.jl")
include("HaldaneSphere.jl")

max_stars_repo_name: mishmash/QHED

max_issues_repo_name: mishmash/QHED

max_forks_repo_name: mishmash/QHED


text:
using MIMMP
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here

C = trunc.(Int,readdlm("C.txt"));
A = trunc.(Int,readdlm("A.txt"));
b = trunc.(Int,readdlm("b.txt"));
d = trunc.(Int,readdlm("d.txt"));
VariableDef = [0 0 400];

@test MIMMPSolver(A,b,C,d,VariableDef,  XLB = zeros(400), XUB = ones(400), Solver = "CB")[1] == 6468
@test MIMMPSolver(A,b,C,d,VariableDef,  XLB = zeros(400), XUB = ones(400), Solver = "DB")[1] == 6468


text:
#TODO: add "expand" method to go from shorter to greater frequencies
#TODO: include functionality for filling the missing values (approx, locf, etc.)

# Date vector

mondays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Monday
tuesdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Tuesday
wednesdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Wednesday
thursdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Thursday
fridays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Friday
saturdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Saturday
sundays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Sunday
# DateTime vector
mondays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Monday
tuesdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Tuesday
wednesdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Wednesday
thursdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Thursday
fridays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Friday
saturdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Saturday
sundays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Sunday

# Time series input

"""
    mondays(x::TS)

Return a time series containing all observations occuring on Mondays of the given `TS` object.
"""
mondays(x::TS) = x[mondays(x.index)]

"""
    tuesdays(x::TS)

Return a time series containing all observations occuring on Tuesdays of the given `TS` object.
"""
tuesdays(x::TS) = x[tuesdays(x.index)]

"""
    wednesdays(x::TS)

Return a time series containing all observations occuring on Wednesdays of the given `TS` object.
"""
wednesdays(x::TS) = x[wednesdays(x.index)]

"""
    thursdays(x::TS)

Return a time series containing all observations occuring on Thursdays of the given `TS` object.
"""
thursdays(x::TS) = x[thursdays(x.index)]

"""
    fridays(x::TS)

Return a time series containing all observations occuring on Fridays of the given `TS` object.
"""
fridays(x::TS) = x[fridays(x.index)]

"""
    saturdays(x::TS)

Return a time series containing all observations occuring on Saturdays of the given `TS` object.
"""
saturdays(x::TS) = x[saturdays(x.index)]

"""
    sundays(x::TS)

Return a time series containing all observations occuring on Sundays of the given `TS` object.
"""
sundays(x::TS) = x[sundays(x.index)]

# If `cal` is true, only returns observations on the last calendar day of the month.
# If `cal` is false, returns all observations for which the next index is a new month.
bow(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== firstdayofweek.(t)) : BitVector([0;diff(week.(t)).!=0])
eow(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== lastdayofweek.(t)) : BitVector([diff(week.(t)).!=0;0])
bom(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== firstdayofmonth.(t)) : BitVector([0;diff(month.(t)).!=0])
eom(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== lastdayofmonth.(t)) : BitVector([diff(month.(t)).!=0;0])
boq(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== firstdayofquarter.(t)) : BitVector([0;diff(quarterofyear.(t)).!=0])
eoq(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== lastdayofquarter(t)) : BitVector([diff(quarterofyear.(t)).!=0;0])
boy(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== firstdayofyear.(t)) : BitVector([0;diff(year.(t)).!=0])
eoy(t::AbstractArray{T,1}; cal::Bool=false) where {T<:TimeType} = cal ? BitVector(t .== lastdayofyear.(t)) : BitVector([diff(year.(t)).!=0;0])

"""
    bow(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the beginnings of the weeks of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the week are returned
- If `cal` is `true`, all observation for which the previous index is a prior week are returned
"""
bow(x::TS; cal::Bool=false) = x[bow(x.index, cal=cal)]

"""
    eow(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the ends of the weeks of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the week are returned
- If `cal` is `true`, all observation for which the next index is a new week are returned
"""
eow(x::TS; cal::Bool=false) = x[eow(x.index, cal=cal)]

"""
    bom(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the beginnings of the months of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the month are returned
- If `cal` is `true`, all observation for which the previous index is a prior month are returned
"""
bom(x::TS; cal::Bool=false) = x[bom(x.index, cal=cal)]

"""
    eom(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the ends of the months of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the month are returned
- If `cal` is `true`, all observation for which the next index is a new month are returned
"""
eom(x::TS; cal::Bool=false) = x[eom(x.index, cal=cal)]

"""
    boq(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the beginnings of the quarters of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the quarter are returned
- If `cal` is `true`, all observation for which the previous index is a prior quarter are returned
"""
boq(x::TS; cal::Bool=false) = x[boq(x.index, cal=cal)]

"""
    eoq(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the ends of the quarters of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the quarter are returned
- If `cal` is `true`, all observation for which the next index is a new quarter are returned
"""
eoq(x::TS; cal::Bool=false) = x[eoq(x.index, cal=cal)]

"""
    boy(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the beginnings of the years of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the year are returned
- If `cal` is `true`, all observation for which the previous index is a prior year are returned
"""
boy(x::TS; cal::Bool=false) = x[boy(x.index, cal=cal)]

"""
    eoy(x::TS; cal::Bool=false)

Return a time series containing all observations occuring at the ends of the years of the input.

- If `cal` is `false`, only observations occurring the last calendar day of the year are returned
- If `cal` is `true`, all observation for which the next index is a new year are returned
"""
eoy(x::TS; cal::Bool=false) = x[eoy(x.index, cal=cal)]

function collapse(x::TS, at::AbstractArray{Bool,1}; fun::Function=last, args...)::TS
    @assert size(at,1) == size(x,1) "Arguments `x` and `at` must have same number of rows."
    idx = findall(at)
    if fun == last
        return x[idx]
    end
    out = zeros((length(idx), size(x,2)))
    @inbounds for j in 1:size(x,2)
        out[1,j] = fun(x.values[1:idx[1], j]; args...)
    end
    @inbounds for j=1:size(x,2), i=2:length(idx)
        a = idx[i-1] + 1
        b = idx[i]
        out[i,j] = fun(x.values[a:b, j]; args...)
    end
    return ts(out, x.index[at], x.fields)
end

"""
    collapse(x::TS, at::Function; fun::Function=last, args...)::TS

Compute the function `fun` over period boundaries defined by `at` (e.g. `eom`, `boq`, etc.) and returng a time series of the results.

Keyword arguments (`args...`) are passed to the function `fun`.
"""
function collapse(x::TS, at::Function; fun::Function=last, args...)::TS
    # Here `at` is a function of the index of `x` which generates a BitArray to be used as above
    return collapse(x, at(x.index), fun=fun; args...)
end

"""
    apply(x::TS, dim::Int=1; fun::Function=sum)

Apply function `fun` across dimension `dim` of a time series `x`.

- If `dim` is 1, then apply `fun` to each row of `x`, returning a time series (`TS`) object of the results
- If `dim` is 2, then apply `fun` to each column of `x`, returning an `Array` of the results
"""
function apply(x::TS, dim::Int=1; fun::Function=sum)
    @assert dim == 1 || dim == 2 "dimension must be either 1 (rows) or 2 (columns)"
    if dim == 1
        n = size(x,1)
        out = zeros(n)
        @inbounds for i in 1:n
            out[i] = fun(x.values[i,:])
        end
        return ts(out, x.index, Symbol(uppercasefirst(string(fun))))
    elseif dim == 2
        k = size(x,2)
        out = zeros(k)
        @inbounds for j in 1:k
            out[j] = fun(x.values[:,j])
        end
        return out'
    else
        error("Argument `dim` must be either 1 (rows) or 2 (columns).")
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ensmallen_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ensmallen")
JLLWrappers.@generate_main_file("ensmallen", UUID("d1a6e830-de90-5fe4-8599-c1a9c9fe1834"))
end  # module ensmallen_jll


text:
#bld_parametric_sin_initial.jl: Parametric sin(): Initial observations
#-------------------------------------------------------------------------------
using CMDimData
using CMDimData.EasyPlot

#==Constants
===============================================================================#
LBL_AXIS_TIME = "Time [s]"
LBL_AXIS_POSITION = "Position [m]"
LBL_AXIS_NORMPOSITION = "Normalized Pos [m/m]"
LBL_AXIS_SPEED = "Speed [m/s]"
lstylesweep = cons(:a, line = set(style=:solid, width=2)) #Default is a bit thin
dfltglyph = cons(:a, glyph = set(shape=:o, size=1.5))

function fnbuild(data)
	plot = cons(:plot, nstrips = 3,
		ystrip1 = set(axislabel=LBL_AXIS_POSITION, striplabel="Sinusoidal response"),
		ystrip2 = set(axislabel=LBL_AXIS_NORMPOSITION, striplabel="Normalized response (All peaks should be ±1)"),
		ystrip3 = set(axislabel=LBL_AXIS_SPEED, striplabel="Rate of change (should be larger for higher frequencies)"),
		xaxis = set(label=LBL_AXIS_TIME)
	)
	push!(plot,
		cons(:wfrm, data.signal, lstylesweep, label="", strip=1),
		cons(:wfrm, data.signal_norm, lstylesweep, label="", strip=2),
		cons(:wfrm, data.rate, lstylesweep, label="", strip=3),
	)
	pcoll = push!(cons(:plotcoll, title="Parametric sin() - Initial Observations"), plot)
	return pcoll
end

#Return EasyPlotBuilder object:
EasyPlot.EasyPlotBuilder(fnbuild)


text:
module TestShow
using Test
using Dates
using CategoricalArrays

@testset "show() for CategoricalPool{String} and its values" begin
    pool = CategoricalPool(["c", "b", "a"])
    opool = CategoricalPool(["c", "b", "a"], true)

    nv1 = CategoricalValue(1, pool)
    nv2 = CategoricalValue(2, pool)
    nv3 = CategoricalValue(3, pool)

    ov1 = CategoricalValue(1, opool)
    ov2 = CategoricalValue(2, opool)
    ov3 = CategoricalValue(3, opool)

    @test sprint(show, pool) == "$CategoricalPool{String,UInt32}([\"c\", \"b\", \"a\"])"
    @test sprint(show, opool) == "$CategoricalPool{String,UInt32}([\"c\", \"b\", \"a\"]) with ordered levels"

    @test sprint(show, nv1) == repr(nv1) == "$CategoricalValue{String,UInt32} \"c\""
    @test sprint(show, nv2) == repr(nv2) == "$CategoricalValue{String,UInt32} \"b\""
    @test sprint(show, nv3) == repr(nv3) == "$CategoricalValue{String,UInt32} \"a\""

    @test sprint(show, ov1) == repr(ov1) =="$CategoricalValue{String,UInt32} \"c\" (1/3)"
    @test sprint(show, ov2) == repr(ov2) == "$CategoricalValue{String,UInt32} \"b\" (2/3)"
    @test sprint(show, ov3) == repr(ov3) =="$CategoricalValue{String,UInt32} \"a\" (3/3)"

    @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "\"c\""
    @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "\"b\""
    @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "\"a\""

    @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "\"c\""
    @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "\"b\""
    @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "\"a\""

    @test sprint(print, nv1) == sprint(print, ov1) == "c"
    @test sprint(print, nv2) == sprint(print, ov2) == "b"
    @test sprint(print, nv3) == sprint(print, ov3) == "a"

    @test string(nv1) == string(ov1) == "c"
    @test string(nv2) == string(ov2) == "b"
    @test string(nv3) == string(ov3) == "a"

    @test String(nv1) == String(ov1) == "c"
    @test String(nv2) == String(ov2) == "b"
    @test String(nv3) == String(ov3) == "a"

    b = IOBuffer()
    @test write(b, nv1) == 1
    @test String(take!(b)) == "c"
    @test write(b, nv2) == 1
    @test String(take!(b)) == "b"
    @test write(b, nv3) == 1
    @test String(take!(b)) == "a"
    @test write(b, ov1) == 1
    @test String(take!(b)) == "c"
    @test write(b, ov2) == 1
    @test String(take!(b)) == "b"
    @test write(b, ov3) == 1
    @test String(take!(b)) == "a"
end

using JSON
@testset "JSON.lower" for pool in (CategoricalPool(["a"]),
                                   CategoricalPool([1]),
                                   CategoricalPool([1.0]))
    v = CategoricalValue(1, pool)
    @test JSON.lower(v) == JSON.lower(get(v))
    @test typeof(JSON.lower(v)) == typeof(JSON.lower(get(v)))
end

using JSON3
using StructTypes
@testset "JSON3.write" begin
    v = CategoricalValue(1, CategoricalPool(["a"]))
    @test JSON3.write(v) === "\"a\""

    v = CategoricalValue(1, CategoricalPool([1]))
    @test JSON3.write(v) === "1"
    @test StructTypes.numbertype(typeof(v)) === Int

    v = CategoricalValue(1, CategoricalPool([2.0]))
    @test JSON3.write(v) === "2.0"
    @test StructTypes.numbertype(typeof(v)) === Float64

    v = CategoricalValue(1, CategoricalPool([BigFloat(3.0,10)]))
    @test JSON3.write(v) === "3.0"
    @test StructTypes.numbertype(typeof(v)) === BigFloat

    v = CategoricalValue(2, CategoricalPool([true,false]))
    @test JSON3.write(v) == "false"
    @test StructTypes.numbertype(typeof(v)) === Bool
end

end


text:
#list of added packages
DataFrames Distributions StatsBase GLM CategoricalArrays Dates NLopt
ForwardDiff Formatting TexTables DataStructures Atom Juno
LinearAlgebra Gadfly HCubature HTTP Measures Formatting ODBC
CodecZlib CSV GZip NLopt CategoricalArrays Serialization Compose
Distributed Random Cairo Fontconfig Compose uCSV Revise StaticArrays

#routines for adding the packages
add Atom Juno
add DataFrames Revise Distributions StatsBase GLM CategoricalArrays Dates NLopt ForwardDiff Formatting
add DataStructures LinearAlgebra Gadfly HCubature HTTP Measures Formatting ODBC CodecZlib
add CSV GZip NLopt CategoricalArrays Serialization Compose Distributed Random Cairo
add Fontconfig Compose uCSV StaticArrays CodecLz4
add https://github.com/clintonTE/Finometrics

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
mutable struct SelfReferential
    o::SelfReferential
    SelfReferential() = (incomplete = new(); incomplete.o = incomplete)
end

alice = SelfReferential()

println(alice === alice)
println(alice === alice.o)
println(alice === alice.o.o)

SelfReferential(o::SelfReferential) = (oo = SelfReferential(); oo.o = o; return oo)

bob = SelfReferential(alice)

println(bob === alice)
println(bob == alice)
println(bob.o == alice)

mutable struct Lazy
    sb
    Lazy(name) = completeit(new(), name)
end

function completeit(o::Lazy, name::AbstractString)
    o.sb = name
    return o
end

cto = Lazy("cto")
println(cto)
println(cto.sb)

# =====OurRational=====

struct OurRational{T<:Integer} <: Real
    numerator::T
    denominator::T
    function OurRational{T}(n::T, d::T) where T<:Integer
        if n == 0 && d == 0
            error("invalid rational: 0//0")
        end
        g = gcd(n, d)
        n = div(n, g)
        d = div(d, g)
        new(n, d)
    end
end

OurRational(n::T, d::T) where {T<:Integer} = OurRational{T}(n, d)

OurRational(n::Integer, d::Integer) = OurRational(promote(n, d)...)

OurRational(n::Integer) = OurRational(n, one(n))

⊘(n::Integer, d::Integer) = OurRational(n, d)

⊘(or::OurRational, k::Integer) = or.numerator ⊘ (or.denominator*k)

⊘(k::Integer, or::OurRational) = (k*or.denominator) ⊘ or.numerator

⊘(k::Complex, m::Real) = complex(real(k) ⊘ m, imag(k) ⊘ m)

⊘(k::Real, m::Complex) = (k*m') ⊘ real(m*m')

function ⊘(p::Complex, q::Complex)
    pq = p*q'
    qq = real(q*q')
    complex(real(pq) ⊘ qq, imag(pq) ⊘ qq)
end

sb = (1+2im) ⊘ (1-2im)
println(typeof(sb))

println(typeof(sb) <: Complex{OurRational})

struct SummedArray{T<:Number,S<:Number}
    data::Vector{T}
    sum::S
    function SummedArray(v::Vector{T}) where T
        S = widen(T)
        new{T,S}(v, sum(S, v))
    end
end

SummedArray(Int32[1,2,3])
# SummedArray(Int32[1,2,3], Int64(6))
# SummedArray(Int32[1,2,3], Int32(6))

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
module StructuresKit

export AISC360_16
include("AISC360_16.jl")
using .AISC360_16

export AISIS10016
include("AISIS10016.jl")
using .AISIS10016

export AISIS10024
include("AISIS10024.jl")
using .AISIS10024

export CUFSM
include("CUFSM.jl")
using .CUFSM

export Eurocode1993
include("Eurocode1993.jl")
using .Eurocode1993

export Geometry
include("Geometry.jl")
using .Geometry

export Mesh
include("Mesh.jl")
using .Mesh

export Connections
include("Connections.jl")
using .Connections

export CrossSection
include("CrossSection.jl")
using .CrossSection

export MaterialModels
include("MaterialModels.jl")
using .MaterialModels

export InternalForces
include("InternalForces.jl")
using .InternalForces

export ThinWalledBeam
include("ThinWalledBeam.jl")
using .ThinWalledBeam

# export Column
# include("Column.jl")
# using .Column

export BeamColumn
include("BeamColumn.jl")
using .BeamColumn

export Visualize
include("Visualize.jl")
using .Visualize

# export PurlinDesigner
# include("PurlinDesigner.jl")
# using .PurlinDesigner

export LSDYNA
include("LSDYNA.jl")
using .LSDYNA

export Truss
include("Truss.jl")
using .Truss

end


text:
module objfunc #NECESSARY module name matching filename without.jl
export objfuncF #NECESSARY the function that the worker will call
function objfuncF(x,y)
    sin(4 * x) + sin(4 * y) + -2 * x + x^2 - 2 * y + y^2
end

end

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

## original example, following here

x=25.6;  y=128.0;
x1=102.4; y1=230.4;
x2=153.6; y2=25.6;
x3=230.4; y3=128.0;

move_to(cr, x, y);
curve_to(cr, x1, y1, x2, y2, x3, y3);

set_line_width(cr, 10.0);
stroke(cr);

set_source_rgba(cr, 1, 0.2, 0.2, 0.6);
set_line_width(cr, 6.0);
move_to(cr,x,y);   line_to(cr,x1,y1);
move_to(cr,x2,y2); line_to(cr,x3,y3);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_curve_to.png");


text:
export BrickletPTCTemperatureCallbackThreshold
struct BrickletPTCTemperatureCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletPTCResistanceCallbackThreshold
struct BrickletPTCResistanceCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletPTCIdentity
struct BrickletPTCIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletPTC
"""
Reads temperatures from Pt100 und Pt1000 sensors
"""
mutable struct BrickletPTC <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletPTC(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_ptc")
        deviceInternal = package.BrickletPTC(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of connected sensor.

If you want to get the temperature periodically, it is recommended
to use the :cb:`Temperature` callback and set the period with
:func:`Set Temperature Callback Period`.
"""
function get_temperature(device::BrickletPTC)
    return device.deviceInternal.get_temperature()
end

export get_resistance
"""
    $(SIGNATURES)

Returns the value as measured by the MAX31865 precision delta-sigma ADC.

The value can be converted with the following formulas:

* Pt100:  resistance = (value * 390) / 32768
* Pt1000: resistance = (value * 3900) / 32768

If you want to get the resistance periodically, it is recommended
to use the :cb:`Resistance` callback and set the period with
:func:`Set Resistance Callback Period`.
"""
function get_resistance(device::BrickletPTC)
    return device.deviceInternal.get_resistance()
end

export set_temperature_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Temperature` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Temperature` callback is only triggered if the temperature has
changed since the last triggering.
"""
function set_temperature_callback_period(device::BrickletPTC, period)
    device.deviceInternal.set_temperature_callback_period(period)
end

export get_temperature_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Temperature Callback Period`.
"""
function get_temperature_callback_period(device::BrickletPTC)
    return device.deviceInternal.get_temperature_callback_period()
end

export set_resistance_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Resistance` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Resistance` callback is only triggered if the resistance has changed
since the last triggering.
"""
function set_resistance_callback_period(device::BrickletPTC, period)
    device.deviceInternal.set_resistance_callback_period(period)
end

export get_resistance_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Resistance Callback Period`.
"""
function get_resistance_callback_period(device::BrickletPTC)
    return device.deviceInternal.get_resistance_callback_period()
end

export set_temperature_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Temperature Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the temperature is *outside* the min and max values"
 "'i'",    "Callback is triggered when the temperature is *inside* the min and max values"
 "'<'",    "Callback is triggered when the temperature is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the temperature is greater than the min value (max is ignored)"
"""
function set_temperature_callback_threshold(device::BrickletPTC, option, min, max)
    device.deviceInternal.set_temperature_callback_threshold(option, min, max)
end

export get_temperature_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Temperature Callback Threshold`.
"""
function get_temperature_callback_threshold(device::BrickletPTC)
    return device.deviceInternal.get_temperature_callback_threshold()
end

export set_resistance_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Resistance Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the temperature is *outside* the min and max values"
 "'i'",    "Callback is triggered when the temperature is *inside* the min and max values"
 "'<'",    "Callback is triggered when the temperature is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the temperature is greater than the min value (max is ignored)"
"""
function set_resistance_callback_threshold(device::BrickletPTC, option, min, max)
    device.deviceInternal.set_resistance_callback_threshold(option, min, max)
end

export get_resistance_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Resistance Callback Threshold`.
"""
function get_resistance_callback_threshold(device::BrickletPTC)
    return device.deviceInternal.get_resistance_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Temperature Reached`,
* :cb:`Resistance Reached`

is triggered, if the threshold

* :func:`Set Temperature Callback Threshold`,
* :func:`Set Resistance Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletPTC, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletPTC)
    return device.deviceInternal.get_debounce_period()
end

export set_noise_rejection_filter
"""
    $(SIGNATURES)

Sets the noise rejection filter to either 50Hz (0) or 60Hz (1).
Noise from 50Hz or 60Hz power sources (including
harmonics of the AC power's fundamental frequency) is
attenuated by 82dB.
"""
function set_noise_rejection_filter(device::BrickletPTC, filter)
    device.deviceInternal.set_noise_rejection_filter(filter)
end

export get_noise_rejection_filter
"""
    $(SIGNATURES)

Returns the noise rejection filter option as set by
:func:`Set Noise Rejection Filter`
"""
function get_noise_rejection_filter(device::BrickletPTC)
    return device.deviceInternal.get_noise_rejection_filter()
end

export is_sensor_connected
"""
    $(SIGNATURES)

Returns *true* if the sensor is connected correctly.

If this function
returns *false*, there is either no Pt100 or Pt1000 sensor connected,
the sensor is connected incorrectly or the sensor itself is faulty.
"""
function is_sensor_connected(device::BrickletPTC)
    return device.deviceInternal.is_sensor_connected()
end

export set_wire_mode
"""
    $(SIGNATURES)

Sets the wire mode of the sensor. Possible values are 2, 3 and 4 which
correspond to 2-, 3- and 4-wire sensors. The value has to match the jumper
configuration on the Bricklet.
"""
function set_wire_mode(device::BrickletPTC, mode)
    device.deviceInternal.set_wire_mode(mode)
end

export get_wire_mode
"""
    $(SIGNATURES)

Returns the wire mode as set by :func:`Set Wire Mode`
"""
function get_wire_mode(device::BrickletPTC)
    return device.deviceInternal.get_wire_mode()
end

export set_sensor_connected_callback_configuration
"""
    $(SIGNATURES)

If you enable this callback, the :cb:`Sensor Connected` callback is triggered
every time a Pt sensor is connected/disconnected.

.. versionadded:: 2.0.2\$nbsp;(Plugin)
"""
function set_sensor_connected_callback_configuration(device::BrickletPTC, enabled)
    device.deviceInternal.set_sensor_connected_callback_configuration(enabled)
end

export get_sensor_connected_callback_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Sensor Connected Callback Configuration`.

.. versionadded:: 2.0.2\$nbsp;(Plugin)
"""
function get_sensor_connected_callback_configuration(device::BrickletPTC)
    return device.deviceInternal.get_sensor_connected_callback_configuration()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletPTC)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletPTC, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Test, Plasmo, MadNLP, MadNLPTests, MadNLPGraph

@testset "MadNLPGraphs test" begin

    n_nodes=2
    M = 200
    d = sin.((0:M*n_nodes).*pi/100)

    graph = OptiGraph()
    @optinode(graph,nodes[1:n_nodes])

    #Node models
    nodecnt = 1
    for (i,node) in enumerate(nodes)
        @variable(node, x[1:M])
        @variable(node, -1<=u[1:M]<=1)
        @constraint(node, dynamics[i in 1:M-1], x[i+1] == x[i] + u[i])
        @objective(node, Min, sum(x[i]^2 - 2*x[i]*d[i+(nodecnt-1)*M] for i in 1:M) +
            sum(u[i]^2 for i in 1:M))
        nodecnt += 1
    end
    n1 = getnode(graph,1)
    @constraint(n1,n1[:x][1] == 0)
    for i=1:n_nodes-1
        @linkconstraint(graph, nodes[i+1][:x][1] == nodes[i][:x][M] + nodes[i][:u][M],attach=nodes[i+1])
    end

    #First node gets initial condition
    kwargs_collection = [
        ("default",Dict(:print_level=>MadNLP.ERROR)),
        ("schur",Dict(:linear_solver=>MadNLPSchur,  :schur_custom_partition=>true,:print_level=>MadNLP.ERROR)),
        ("schwarz",Dict(:linear_solver=>MadNLPSchwarz,:schwarz_custom_partition=>true,:print_level=>MadNLP.ERROR))
    ]
    for (name,kwargs) in kwargs_collection
        @testset "plasmo-$name" begin
            MadNLP.optimize!(graph;kwargs...);
            @test MadNLP.termination_status(graph.optimizer) == MOI.LOCALLY_SOLVED
            @test solcmp([0.0,0.03137979101284875,0.0627286139604959,0.09401553133948139,0.12520966673966746,0.15628023531552773,0.18719657416707308,0.21792817260043182,0.24844470223822107,0.2787160469499936], value.(getnode(graph,1)[:x][1:10]))
            @test solcmp([-0.0627595821831224,-0.06269764605256355,-0.06257383491492904,-0.06238827095686386,-0.06214113730759117,-0.06183267785818605,-0.06146319702088434,-0.06103305942866446,-0.06054268957539879,-0.05999257139692946], dual.(getnode(graph,1)[:dynamics][1:10]))
            @test solcmp([-1.574246826343827e-10,-1.5726901766520537e-10,-1.5695784406668639e-10,-1.5649147434189076e-10,-1.5587037679839438e-10,-1.5509517501550025e-10,-1.541666471342132e-10,-1.5308572497187857e-10,-1.5185349296508524e-10,-1.504711869423148e-10], dual.(UpperBoundRef.(getnode(graph,1)[:u][1:10])))
            @test solcmp([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], dual.(LowerBoundRef.(getnode(graph,1)[:u][1:10])))
        end
    end
end


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

@doc raw"""MetricTarget defines the target value, average value, or average utilization of a specific metric

    IoK8sApiAutoscalingV2beta2MetricTarget(;
        averageUtilization=nothing,
        averageValue=nothing,
        type=nothing,
        value=nothing,
    )

    - averageUtilization::Int32 : averageUtilization is the target value of the average of the resource metric across all relevant pods, represented as a percentage of the requested value of the resource for the pods. Currently only valid for Resource metric source type
    - averageValue::IoK8sApimachineryPkgApiResourceQuantity : averageValue is the target value of the average of the metric across all relevant pods (as a quantity)
    - type::String : type represents whether the metric type is Utilization, Value, or AverageValue
    - value::IoK8sApimachineryPkgApiResourceQuantity : value is the target value of the metric (as a quantity).
"""
mutable struct IoK8sApiAutoscalingV2beta2MetricTarget <: SwaggerModel
    averageUtilization::Any # spec type: Union{ Nothing, Int32 } # spec name: averageUtilization
    averageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: averageValue
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    value::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: value

    function IoK8sApiAutoscalingV2beta2MetricTarget(;averageUtilization=nothing, averageValue=nothing, type=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("averageUtilization"), averageUtilization)
        setfield!(o, Symbol("averageUtilization"), averageUtilization)
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("averageValue"), averageValue)
        setfield!(o, Symbol("averageValue"), averageValue)
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiAutoscalingV2beta2MetricTarget

const _property_map_IoK8sApiAutoscalingV2beta2MetricTarget = Dict{Symbol,Symbol}(Symbol("averageUtilization")=>Symbol("averageUtilization"), Symbol("averageValue")=>Symbol("averageValue"), Symbol("type")=>Symbol("type"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiAutoscalingV2beta2MetricTarget = Dict{Symbol,String}(Symbol("averageUtilization")=>"Int32", Symbol("averageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("type")=>"String", Symbol("value")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2MetricTarget }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2MetricTarget))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2MetricTarget }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2MetricTarget[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2MetricTarget }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2MetricTarget[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2MetricTarget)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

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


text:
# generalise this into coupler-specific function

abstract type AbstractCheck end

struct ConservationCheck{A} <: AbstractCheck
    ρe_tot_atmos::A
    ρe_tot_slab::A
end

function check_conservation_callback(cs, atmos_sim, slab_sim)

    atmos_sim.domain.center_space.vertical_topology.mesh.domain
    Δz_1 =
        parent(ClimaCore.Fields.coordinate_field(atmos_sim.domain.face_space).z)[2] -
        parent(ClimaCore.Fields.coordinate_field(atmos_sim.domain.face_space).z)[1]

    atmos_field = atmos_sim.integrator.u.c.ρe  # J
    slab_field = get_slab_energy(slab_sim, slab_sim.integrator.u.T_sfc) ./ Δz_1  # J  [NB: sum of the boundary field inherits the depth from the first atmospheric layer, which ≂̸ slab depth]

    ρe_tot_atmos = sum(atmos_field)
    ρe_tot_slab = sum(slab_field)

    push!(cs.ρe_tot_atmos, ρe_tot_atmos)
    push!(cs.ρe_tot_slab, ρe_tot_slab)
end

get_slab_energy(slab_sim, T_sfc) = slab_sim.params.ρ .* slab_sim.params.c .* T_sfc .* slab_sim.params.h #[NB: upon initialisation FT(275) was subtracted from T_0]

using ClimaCorePlots
function plot(CS::ConservationCheck)
    diff_ρe_tot_atmos = CS.ρe_tot_atmos .- CS.ρe_tot_atmos[3]
    diff_ρe_tot_slab = (CS.ρe_tot_slab .- CS.ρe_tot_slab[3])
    Plots.plot(diff_ρe_tot_atmos, label = "atmos")
    Plots.plot!(diff_ρe_tot_slab, label = "slab")
    tot = diff_ρe_tot_atmos .+ diff_ρe_tot_slab
    Plots.plot!(tot .- tot[1], label = "tot")
end

function conservation_plot(atmos_sim, slab_sim, solu_atm, solu_slab, figname = "tst_c.png")

    atmos_e = [sum(u.c.ρe) for u in solu_atm] # J

    z = parent(ClimaCore.Fields.coordinate_field(atmos_sim.domain.face_space).z)
    Δz_1 = z[2] - z[1]
    slab_e = [sum(get_slab_energy(slab_sim, u)) for u in solu_slab]  # J  [NB: sum of the boundary field inherits the depth from the first atmospheric layer, which ≂̸ slab depth]

    diff_ρe_tot_atmos = atmos_e .- atmos_e[3]
    diff_ρe_tot_slab = (slab_e .- slab_e[3])
    Plots.plot(diff_ρe_tot_atmos, label = "atmos")
    Plots.plot!(diff_ρe_tot_slab, label = "slab")
    tot = diff_ρe_tot_atmos .+ diff_ρe_tot_slab
    Plots.plot!(tot .- tot[1], label = "tot", xlabel = "time [s]", ylabel = "energy(t) - energy(t=0) [s]")
    Plots.savefig(figname)

end

max_stars_repo_name: CliMA/ClimaCoupler


text:
## Exercise 5-5
## What is the output of the following program? Draw a stack diagram that shows the state of the program when it prints the result.

function recurse(n, s)
    if n == 0
        println(s)
    else
        recurse(n-1, n+s)
    end
end

## recurse(3, 0)
println("Ans 0: ")

println(" n --> 3. s --> 0")
println(" n --> 2. s --> 3")
println(" n --> 1. s --> 5")
println(" n --> 0. s --> 6")

## 1. What would happen if you called this function like this: recurse(-1, 0)?
println("Ans 1: ")

println("if the function was called like this: recurse(-1, 0), it would get stuck in infinite recursion")

## 2. Write a docstring that explains everything someone would need to know in order to use this function (and nothing else).
println("Ans 2: ")

"""
recurse(n::Int, s::Number) -> Nothing

recurses for n number of times and sums it up
adding initial value of sum (s).

@params n: positive integer
@params s: number
"""
function recurse(n, s)
    if n == 0
        println(s)
    else
        recurse(n-1, n+s)
    end
end

recurse(9, 2)

println("End.")


text:
module MDTest

using Test, PIMD

@testset "MD" begin
    phase = HamiltonianPhase([1., 2.], [3., 4.])
    @test phase.q ≈ [1., 2.]
    phase.p .= 0
    @test phase.p ≈ [0., 0.]
end

end


text:
"""
    GraphicConfig

    GraphicConfig(; fps=25.0, store=false)

Config and variables for graphic outputs.
"""
mutable struct GraphicConfig{FPS,TS}
    fps::FPS
    store::Bool
    timestamp::TS
    stampframe::Int
    stoppedframe::Int
end
GraphicConfig(; fps=25.0, store=false, kw...) = GraphicConfig(fps, store, 0.0, 1, 1)

fps(gc::GraphicConfig) = gc.fps
store(gc::GraphicConfig) = gc.store
timestamp(gc::GraphicConfig) = gc.timestamp
stampframe(gc::GraphicConfig) = gc.stampframe
stoppedframe(gc::GraphicConfig) = gc.stoppedframe
setfps!(gc::GraphicConfig, x) = gc.fps = x
function settimestamp!(o::GraphicConfig, frame::Int)
    o.timestamp = time()
    o.stampframe = frame
end

setstoppedframe!(gc::GraphicConfig, frame::Int) = gc.stoppedframe = frame

"""
    GraphicOutput <: Output

Abstract supertype for [`Output`](@ref)s that display the simulation frames.

All `GraphicOutputs` must have a [`GraphicConfig`](@ref) object
and define a [`showframe`](@ref) method.

See [`REPLOutput`](@ref) for an example.

## Keywords:

The default constructor will generate these objects from other keyword arguments
and pass them to the object constructor, which must accept the following:

- `frames`: a `Vector` of simulation frames (`NamedTuple` or `Array`).
- `running`: A `Bool`.
- `extent` an [`Extent`](@ref) object.
- `graphicconfig` a [`GraphicConfig`](@ref)object.

"""
abstract type GraphicOutput{T,F} <: Output{T,F} end

# Generic ImageOutput constructor. Converts an init array to vector of arrays.
function (::Type{T})(
    init::Union{NamedTuple,AbstractMatrix};
    extent=nothing, graphicconfig=nothing, kw...
) where T <: GraphicOutput
    extent = extent isa Nothing ? Extent(; init=init, kw...) : extent
    graphicconfig = graphicconfig isa Nothing ? GraphicConfig(; kw...) : graphicconfig
    T(; frames=[deepcopy(init)], running=false,
      extent=extent, graphicconfig=graphicconfig, kw...)
end

graphicconfig(o::Output) = GraphicConfig()
graphicconfig(o::GraphicOutput) = o.graphicconfig

# Forwarded getters and setters
fps(o::GraphicOutput) = fps(graphicconfig(o))
timestamp(o::GraphicOutput) = timestamp(graphicconfig(o))
stampframe(o::GraphicOutput) = stampframe(graphicconfig(o))
stoppedframe(o::GraphicOutput) = stoppedframe(graphicconfig(o))
isstored(o::GraphicOutput) = store(o)
store(o::GraphicOutput) = store(graphicconfig(o))

setfps!(o::GraphicOutput, x) = setfps!(graphicconfig(o), x)
settimestamp!(o::GraphicOutput, f) = settimestamp!(graphicconfig(o), f)
setstoppedframe!(o::GraphicOutput, f) = setstoppedframe!(graphicconfig(o), f)

# Delay output to maintain the frame rate
delay(o::GraphicOutput, f) =
    sleep(max(0.0, timestamp(o) + (f - stampframe(o))/fps(o) - time()))
isshowable(o::GraphicOutput, f) = true

function storeframe!(o::GraphicOutput, data)
    f = frameindex(o, data)
    if f > length(o)
        _pushgrid!(eltype(o), o)
    end
    if isstored(o)
        _storeframe!(eltype(o), o, data)
    end
    if isshowable(o, currentframe(data))
        showframe(o, data)
    end
    return nothing
end

_pushgrid!(::Type{<:NamedTuple}, o) = push!(o, map(grid -> similar(grid), o[1]))
_pushgrid!(::Type{<:AbstractArray}, o) = push!(o, similar(o[1]))

showframe(o::GraphicOutput, data) = showframe(o, proc(data), data)
showframe(o::GraphicOutput, ::Processor, data) = showframe(map(gridview, grids(data)), o, data)
# Handle NamedTuple for outputs that only accept AbstractArray
showframe(frame::NamedTuple, o::GraphicOutput, data) = showframe(first(frame), o, data)

function initialise!(o::GraphicOutput, data)
    initalisegraphics(o, data)
end
function finalise!(o::GraphicOutput, data)
    _storeframe!(eltype(o), o, data)
    finalisegraphics(o, data)
end

initalisegraphics(o::GraphicOutput, data) = nothing
finalisegraphics(o::GraphicOutput, data) = showframe(o, data)


text:
using PyPlot
using Gen

include("poisson_process.jl")
include("load_data.jl")

# Example from Section 4 of Reversible jump Markov chain Monte Carlo
# computation and Bayesian model determination

# minimum of k draws from uniform_continuous(lower, upper)

# we can sequentially sample the order statistics of a collection of K uniform
# continuous samples on the interval [a, b], by:
# x1 ~ min_uniform_continuous(a, b, K)
# x2 | x1 ~ min_uniform_continuous(x1, b, K-1)
# ..
# xK | x1 .. x_{K-1} ~ min_uniform_continuous(x_{K-1}, b, 1)

struct MinUniformContinuous <: Distribution{Float64} end
const min_uniform_continuous = MinUniformContinuous()

function Gen.logpdf(::MinUniformContinuous, x::Float64, lower::T, upper::U, k::Int) where {T<:Real,U<:Real}
    if x > lower && x < upper
        (k-1) * log(upper - x) + log(k) - k * log(upper - lower)
    else
        -Inf
    end
end

function Gen.random(::MinUniformContinuous, lower::T, upper::U, k::Int) where {T<:Real,U<:Real}
    # inverse CDF method
    p = rand()
    upper - (upper - lower) * (1. - p)^(1. / k)
end

#########
# model #
#########

const K = :k
const EVENTS = :events
const CHANGEPT = :changept
const RATE = :rate

@gen function model(T::Float64)

    # prior on number of change points
    k = @trace(poisson(3.), K)

    # prior on the location of (sorted) change points
    change_pts = Vector{Float64}(undef, k)
    lower = 0.
    for i=1:k
        cp = @trace(min_uniform_continuous(lower, T, k-i+1), (CHANGEPT, i))
        change_pts[i] = cp
        lower = cp
    end

    # k + 1 rate values
    # h$i is the rate for cp$(i-1) to cp$i where cp0 := 0 and where cp$(k+1) := T
    alpha = 1.
    beta = 200.
    rates = Float64[@trace(Gen.gamma(alpha, 1. / beta), (RATE, i)) for i=1:k+1]

    # poisson process
    bounds = vcat([0.], change_pts, [T])
    @trace(piecewise_poisson_process(bounds, rates), EVENTS)
end

function render(trace; ymax=0.02)
    T = get_args(trace)[1]
    k = trace[:k]
    bounds = vcat([0.], sort([trace[(CHANGEPT, i)] for i=1:k]), [T])
    rates = [trace[(RATE, i)] for i=1:k+1]
    for i=1:length(rates)
        lower = bounds[i]
        upper = bounds[i+1]
        rate = rates[i]
        plot([lower, upper], [rate, rate], color="black", linewidth=2)
    end
    points = trace[EVENTS]
    scatter(points, -rand(length(points)) * (ymax/5.), color="black", s=5)
    ax = gca()
    xlim = [0., T]
    plot(xlim, [0., 0.], "--")
    ax[:set_xlim](xlim)
    ax[:set_ylim](-ymax/5., ymax)
end

function show_prior_samples()
    figure(figsize=(16,16))
    T = 40000.
    for i=1:16
        println("simulating $i")
        subplot(4, 4, i)
        (trace, ) = generate(model, (T,), EmptyChoiceMap())
        render(trace; ymax=0.015)
    end
    tight_layout(pad=0)
    savefig("prior_samples.pdf")
end

#############
# rate move #
#############

@gen function rate_proposal(trace)

    # pick a random segment whose rate to change
    i = @trace(uniform_discrete(1, trace[K]+1), :i)

    # propose new value for the rate
    cur_rate = trace[(RATE, i)]
    @trace(uniform_continuous(cur_rate/2., cur_rate*2.), :new_rate)

    nothing
end

# it is an involution because it:
# - maintains i = fwd_choices[:i] constant
# - swaps choices[(RATE, i)] with fwd_choices[:new_rate]

@transform rate_involution() (model_in, aux_in) to (model_out, aux_out) begin
    i = @read(aux_in[:i], :discrete)
    @write(aux_out[:i], i, :discrete)
    new_rate = @read(aux_in[:new_rate], :continuous)
    @write(model_out[(RATE, i)], new_rate, :continuous)
    prev_rate = @read(model_in[(RATE, i)], :continuous)
    @write(aux_out[:new_rate], prev_rate, :continuous)
end

is_involution!(rate_involution)

rate_move(trace) = metropolis_hastings(trace, rate_proposal, (), rate_involution, check=true)

#################
# position move #
#################

@gen function position_proposal(trace)
    k = trace[K]
    @assert k > 0

    # pick a random changepoint to change
    i = @trace(uniform_discrete(1, k), :i)

    lower = (i == 1) ? 0. : trace[(CHANGEPT, i-1)]
    upper = (i == k) ? T : trace[(CHANGEPT, i+1)]
    @trace(uniform_continuous(lower, upper), :new_changept)

    i
end

# it is an involution because it:
# - maintains i = fwd_choices[:i] constant
# - swaps choices[(CHANGEPT, i)] with fwd_choices[:new_changept]

@transform position_involution (model_in, aux_in) to (model_out, aux_out) begin
    i = @read(aux_in[:i], :discrete)
    @write(aux_out[:i], i, :discrete)
    @copy(model_in[(CHANGEPT, i)], aux_out[:new_changept])
    @copy(aux_in[:new_changept], model_out[(CHANGEPT, i)])
end

is_involution!(position_involution)

position_move(trace) = metropolis_hastings(trace, position_proposal, (), position_involution, check=true)

######################
# birth / death move #
######################

const CHOSEN = :chosen
const IS_BIRTH = :is_birth
const NEW_CHANGEPT = :new_changept
const U = :u

@gen function birth_death_proposal(trace)
    T = get_args(trace)[1]
    k = trace[K]

    # if k = 0, then always do a birth move
    # if k > 0, then randomly choose a birth or death move
    isbirth = (k == 0) ? true : @trace(bernoulli(0.5), IS_BIRTH)

    if isbirth
        # pick the segment in which to insert the new changepoint
        # changepoints before move:  | 1     2    3 |
        # new changepoint (i = 2):   |    *         |
        # changepoints after move:   | 1  2  3    4 |
        i = @trace(uniform_discrete(1, k+1), CHOSEN)
        lower = (i == 1) ? 0. : trace[(CHANGEPT, i-1)]
        upper = (i == k+1) ? T : trace[(CHANGEPT, i)]
        @trace(uniform_continuous(lower, upper), NEW_CHANGEPT)
        @trace(uniform_continuous(0., 1.), U)
    else
        # pick the changepoint to be deleted
        # changepoints before move:     | 1  2  3    4 |
        # deleted changepoint (i = 2):  |    *         |
        # changepoints after move:      | 1     2    3 |
        @trace(uniform_discrete(1, k), CHOSEN)
    end
    nothing
end

function new_rates(cur_rate, u, cur_cp, prev_cp, next_cp)
    d_prev = cur_cp - prev_cp
    d_next = next_cp - cur_cp
    @assert d_prev > 0
    @assert d_next > 0
    d_total = d_prev + d_next
    log_cur_rate = log(cur_rate)
    log_ratio = log(1 - u) - log(u)
    prev_rate = exp(log_cur_rate - (d_next / d_total) * log_ratio)
    next_rate = exp(log_cur_rate + (d_prev / d_total) * log_ratio)
    @assert prev_rate > 0.
    @assert next_rate > 0.
    (prev_rate, next_rate)
end

function new_rates_inverse(prev_rate, next_rate, cur_cp, prev_cp, next_cp)
    d_prev = cur_cp - prev_cp
    d_next = next_cp - cur_cp
    @assert d_prev > 0
    @assert d_next > 0
    d_total = d_prev + d_next
    log_prev_rate = log(prev_rate)
    log_next_rate = log(next_rate)
    cur_rate = exp((d_prev / d_total) * log_prev_rate + (d_next / d_total) * log_next_rate)
    u = prev_rate / (prev_rate + next_rate)
    @assert cur_rate > 0.
    (cur_rate, u)
end

# it is an involution because:
# - it switches back and forth between birth move and death move
# - it maintains fwd_choices[CHOSEN] constant (applying it twice will first insert a new
#   changepoint and then remove that same changepoint)
# - new_rates, curried on cp_new, cp_prev, and cp_next, is the inverse of new_rates_inverse.

@transform birth_death_involution (model_in, aux_in) to (model_out, aux_out) begin
    T = get_args(model_in)[1]

    # current number of changepoints
    k = @read(model_in[K], :discrete)

    # if k == 0, then we can only do a birth move
    isbirth = (k == 0) || @read(aux_in[IS_BIRTH], :discrete)

    # if we are a birth move, the inverse is a death move
    if k > 1 || isbirth
        @write(aux_out[IS_BIRTH], !isbirth, :discrete)
    end

    # the changepoint to be added or deleted
    i = @read(aux_in[CHOSEN], :discrete)
    @copy(aux_in[CHOSEN], aux_out[CHOSEN])

    if isbirth
        @tcall(birth(k, i))
    else
        @tcall(death(k, i))
    end
end

@transform birth(k::Int, i::Int) (model_in, aux_in) to (model_out, aux_out) begin
    @write(model_out[K], k+1, :discrete)

    cp_new = @read(aux_in[NEW_CHANGEPT], :continuous)
    cp_prev = (i == 1) ? 0. : @read(model_in[(CHANGEPT, i-1)], :continuous)
    cp_next = (i == k+1) ? T : @read(model_in[(CHANGEPT, i)], :continuous)

    # set new changepoint
    @copy(aux_in[NEW_CHANGEPT], model_out[(CHANGEPT, i)])

    # shift up changepoints
    for j=i+1:k+1
        @copy(model_in[(CHANGEPT, j-1)], model_out[(CHANGEPT, j)])
    end

    # compute new rates
    h_cur = @read(model_in[(RATE, i)], :continuous)
    u = @read(aux_in[U], :continuous)
    (h_prev, h_next) = new_rates(h_cur, u, cp_new, cp_prev, cp_next)

    # set new rates
    @write(model_out[(RATE, i)], h_prev, :continuous)
    @write(model_out[(RATE, i+1)], h_next, :continuous)

    # shift up rates
    for j=i+2:k+2
        @copy(model_in[(RATE, j-1)], model_out[(RATE, j)])
    end
end

@transform death(k::Int, i::Int) (model_in, aux_in) to (model_out, aux_out) begin
    @write(model_out[K], k-1, :discrete)

    cp_deleted = @read(model_in[(CHANGEPT, i)], :continuous)
    cp_prev = (i == 1) ? 0. : @read(model_in[(CHANGEPT, i-1)], :continuous)
    cp_next = (i == k) ? T : @read(model_in[(CHANGEPT, i+1)], :continuous)
    @copy(model_in[(CHANGEPT, i)], aux_out[NEW_CHANGEPT])

    # shift down changepoints
    for j=i:k-1
        @copy(model_in[(CHANGEPT, j+1)], model_out[(CHANGEPT, j)])
    end

    # compute cur rate and u
    h_prev = @read(model_in[(RATE, i)], :continuous)
    h_next = @read(model_in[(RATE, i+1)], :continuous)
    (h_cur, u) = new_rates_inverse(h_prev, h_next, cp_deleted, cp_prev, cp_next)
    @write(aux_out[U], u, :continuous)

    # set cur rate
    @write(model_out[(RATE, i)], h_cur, :continuous)

    # shift down rates
    for j=i+1:k
        @copy(model_in[(RATE, j+1)], model_out[(RATE, j)])
    end
end

is_involution!(birth_death_involution)

birth_death_move(trace) = metropolis_hastings(trace, birth_death_proposal, (), birth_death_involution, check=true)

function mcmc_step(trace)
    (trace, _) = rate_move(trace)
    if trace[K] > 0
        (trace, _) = position_move(trace)
    end
    (trace, _) = birth_death_move(trace)
    trace
end

function simple_mcmc_step(trace)
    (trace, _) = rate_move(trace)
    if trace[K] > 0
        (trace, _) = position_move(trace)
    end
    (trace, _) = metropolis_hastings(trace, k_selection)
    trace
end

function do_mcmc(T, num_steps::Int)
    (trace, _) = generate(model, (T,), observations)
    for iter=1:num_steps
        k = trace[K]
        if iter % 1000 == 0
            println("iter $iter of $num_steps, k: $k")
        end
        trace = mcmc_step(trace)
    end
    trace
end

const k_selection = select(K)

function do_simple_mcmc(T, num_steps::Int)
    (trace, _) = generate(model, (T,), observations)
    for iter=1:num_steps
        k = trace[K]
        if iter % 1000 == 0
            println("iter $iter of $num_steps, k: $k")
        end
        trace = simple_mcmc_step(trace)
    end
    trace
end

########################
# inference experiment #
########################

import Random
Random.seed!(1)

const points = load_data_set()
const T = maximum(points)
const observations = choicemap()
observations[EVENTS] = points

function show_posterior_samples()
    figure(figsize=(16,16))
    for i=1:16
        println("replicate $i")
        subplot(4, 4, i)
        trace = do_mcmc(T, 200)
        render(trace; ymax=0.015)
    end
    tight_layout(pad=0)
    savefig("posterior_samples.pdf")

    figure(figsize=(16,16))
    for i=1:16
        println("replicate $i")
        subplot(4, 4, i)
        trace = do_simple_mcmc(T, 200)
        render(trace; ymax=0.015)
    end
    tight_layout(pad=0)
    savefig("posterior_samples_simple.pdf")
end

function get_rate_vector(trace, test_points)
    k = trace[K]
    cps = [trace[(CHANGEPT, i)] for i=1:k]
    hs = [trace[(RATE, i)] for i=1:k+1]
    rate = Vector{Float64}()
    cur_h_idx = 1
    cur_h = hs[cur_h_idx]
    next_cp_idx = 1
    upper = (next_cp_idx == k + 1) ? T : cps[next_cp_idx]
    for x in test_points
        while x > upper
            next_cp_idx += 1
            upper = (next_cp_idx == k + 1) ? T : cps[next_cp_idx]
            cur_h_idx += 1
            cur_h = hs[cur_h_idx]
        end
        push!(rate, cur_h)
    end
    rate
end

# compute posterior mean rate curve

function plot_posterior_mean_rate()
    test_points = collect(1.0:10.0:T)
    rates = Vector{Vector{Float64}}()
    num_samples = 0
    num_steps = 8000
    for reps=1:20
        (trace, _) = generate(model, (T,), observations)
        for iter=1:num_steps
            if iter % 1000 == 0
                println("iter $iter of $num_steps, k: $(trace[K])")
            end
            trace = mcmc_step(trace)
            if iter > 4000
                num_samples += 1
                rate_vector = get_rate_vector(trace, test_points)
                @assert length(rate_vector) == length(test_points)
                push!(rates, rate_vector)
            end
        end
    end
    posterior_mean_rate = zeros(length(test_points))
    for rate in rates
        posterior_mean_rate += rate / Float64(num_samples)
    end
    ymax = 0.010
    figure()
    plot(test_points, posterior_mean_rate, color="black")
    scatter(points, -rand(length(points)) * (ymax/6.), color="black", s=5)
    ax = gca()
    xlim = [0., T]
    plot(xlim, [0., 0.], "--")
    ax[:set_xlim](xlim)
    ax[:set_ylim](-ymax/5., ymax)
    savefig("posterior_mean_rate.pdf")
end

function plot_trace_plot()
    figure(figsize=(8, 4))

    # reversible jump
    (trace, _) = generate(model, (T,), observations)
    rate1 = Float64[]
    num_clusters_vec = Int[]
    burn_in = 0
    for iter=1:burn_in + 1000
        trace = mcmc_step(trace)
        if iter > burn_in
            push!(num_clusters_vec, trace[K])
        end
    end
    subplot(2, 1, 1)
    plot(num_clusters_vec, "b")

    # simple MCMC
    (trace, _) = generate(model, (T,), observations)
    rate1 = Float64[]
    num_clusters_vec = Int[]
    burn_in = 0
    for iter=1:burn_in + 1000
        trace = simple_mcmc_step(trace)
        if iter > burn_in
            push!(num_clusters_vec, trace[K])
        end
    end
    subplot(2, 1, 2)
    plot(num_clusters_vec, "b")

    savefig("trace_plot.pdf")
end

println("showing prior samples...")
show_prior_samples()

println("showing posterior samples...")
show_posterior_samples()

println("estimating posterior mean rate...")
plot_posterior_mean_rate()

println("making trace plot...")
plot_trace_plot()


text:
using ModelingToolkit
MTK = ModelingToolkit

# function MTK.to_diffeq(eq::Equation)
#     if isintermediate(eq)
#
#     end
#     (x, t, n) = flatten_differential(eq.lhs)
#     (isa(t, Operation) && isa(t.op, Variable) && isempty(t.args)) ||
#         throw(ArgumentError("invalid independent variable $t"))
#     (isa(x, Operation) && isa(x.op, Variable) && length(x.args) == 1 && isequal(first(x.args), t)) ||
#         throw(ArgumentError("invalid dependent variable $x"))
#     return t.op, DiffEq(x.op, n, eq.rhs)
# end
#
# function DAESystem(eqs)
#     reformatted = to_dae.(eqs)
# end
isnamedin(varname, varlist) = reduce((x1,x2)-> x1 || x2.name==varname, varlist; init=false)

abstract type Reference <: Expression end

struct VarRef <: Reference
    sys::MTK.AbstractSystem
    ref
end

struct ParamRef <: Reference
    sys::MTK.AbstractSystem
    ref
end

struct Component
    name::Symbol
    subcomponents::Array{Component,1}
    connections::Array{Equation,1}
    sys::Union{ODESystem, Nothing}
    function Component(name, eqs::Array{Equation,1})
        sc, con, sys = parse(eqs)
        return new(Symbol(name), sc, con, sys)
    end
end

macro component(sys::ODESystem)

end

function Base.parse(eqs::Array{Equation,1})
    sc = get_systems(eqs)

end

struct SuperSystem <: MTK.AbstractSystem
    systems::Array{MTK.AbstractSystem,1}
    con_eqs::Array{Equation,1}
    eqs::Array{MTK.DiffEq,1}
    iv::Variable
    dvs::Array{Variable,1}
    ps::Array{Variable,1}
    jac::Base.RefValue{Array{Expression,2}}
    Wfact::Base.RefValue{Array{Expression,2}}
    Wfact_t::Base.RefValue{Array{Expression,2}}
end
function SuperSystem(eqs::Array{Equation,1})
end

get_systems(x) = []
get_systems(ode::ODESystem) = [ode]
get_systems(ref::Reference) = [ref.sys]
get_systems(op::Operation) = get_systems(op.args)
get_systems(arr::Array) = vcat(get_systems.(arr)...) |> unique
function get_systems(eq::Equation)
    @assert isa(eq.lhs, VarRef) "lhs must be a VarRef. DAEs aren't supported yet."
    return get_systems([eq.lhs, eq.rhs]) |> Array{MTK.AbstractSystem}
end

Base.convert(::Type{Expr}, ref::Reference) =
    MTK.build_expr(:call, Any[Symbol])

function Base.getindex(sys::ODESystem, key)
    if isnamedin(key, sys.ps)
        return ParamRef(sys, key)
    elseif isnamedin(key, sys.dvs)
        return VarRef(sys, key)
    else
        throw(ArgumentError("invalid index: $key"))
    end
end

function Base.setindex!(sys::ODESystem, val::Reference, ref)
end

function Base.collect(sys::SuperSystem)

end

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
using OrdinaryDiffEq, FractionalDiffEq, Plots

function fun(du, u, p, x)
    du=u+2*x^2-x-3
end

prob = SecondOrderODEProblem(fun, 0, 0, (0, 5))
sol = OrdinaryDiffEq.solve(prob, Tsit5())

fodeprob = SingleTermFODEProblem((x, u)->u+2*x^2-x-3, 2, 0, (0, 5))
fodesol = FractionalDiffEq.solve(fodeprob, 0.001, PECE())

analytical(x)=-2*x^2+x+exp(-x)-1

plot(sol, vars=(2))
plot!(fodesol, ls=:dash)

rightfun(t, x) = 0
multiprob = MultiTermsFODEProblem([1, -9, 20], [2, 1, 0], rightfun, [1, 2], (0, 1))
multisol = FractionalDiffEq.solve(multiprob, 0.01, PIPECE())
anal(x)=exp(4*x)*(3-2*exp(x))

analytics=anal.(collect(0:0.01:1))
plot(multisol)
plot!(collect(0:0.01:1), analytics, ls=:dash)