text:
# f1 = 3*x1-cos(x2*x3)-0.5
# f2 = x1^2-81*(x2+0.1)^2+sin(x3)+1.06
# f3 = e^(-x1*x2)+20*x3+(10*π-3)/3

x1 = 0.01
x2 = 0.01
x3 = -0.01
f1 = 3*x1-cos(x2*x3)-0.5
f2 = x1^2-81*(x2+0.1)^2+sin(x3)+1.06
f3 = e^(-x1*x2)+20*x3+(10*π-3)/3

x10 = x1
x20 = x2
x30 = x3

dx1 = 0.01
dx2 = 0.01
dx3 = 0.01
x1 = x1 + dx1
x2 = x2 + dx2
x3 = x3 + dx3

while (abs(f1)>1e-6)|(abs(f2)>1e-6)|(abs(f3)>1e-6)
  f10 = f1
  f20 = f2
  f30 = f3

  f1 = 3*x1-cos(x2*x3)-0.5
  f2 = x1^2-81*(x2+0.1)^2+sin(x3)+1.06
  f3 = e^(-x1*x2)+20*x3+(10*π-3)/3

  x11 = x1
  x22 = x2
  x33 = x3

  x1 = x1 + -f1/(f1-f10)
  x2 = x2 + -f1/(f1-f10)
  x3 = x3 + -f1/(f1-f10)

  x10 = x11
  x20 = x22
  x30 = x33
end

print("x1 is : $(x1)\n")
print("x2 is : $(x2)\n")
print("x3 is : $(x3)\n")

ext: jl

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
function RTBPEnergy(u, mu)
"""
     RTBP-Nbody problem
"""

 @inbounds begin
    x = u[1]   # x
    y = u[2]   # y
    px = u[3]  # px
    py = u[4]  # py

    umu=1-mu

    r1 = ((x+mu)^2 + y^2)^(1/2)
    r2 = ((x-umu)^2 + y^2)^(1/2)

    Energy=(px*px+py*py)/2+px*y-py*x-umu/r1-mu/r2-mu*umu/2

    return(Energy)
    end

end

f = (du,u,p,t) -> begin

  μ=p[1]
  uμ=1-μ
  @inbounds begin
  # 1 = y₁
  # 2 = y₂
  # 3 = y₁'
  # 4 = y₂'
  D₁ = ((u[1]+μ)^2 + u[2]^2)^(3/2)
  D₂ = ((u[1]-uμ)^2 + u[2]^2)^(3/2)
  du[1] = u[3]+u[2]
  du[2] = u[4]-u[1]
  du[3] = u[4] - uμ*(u[1]+μ)/D₁ - μ*(u[1]-uμ)/D₂
  du[4] = -u[3] - uμ*u[2]/D₁ - μ*u[2]/D₂
  end
end


text:
# benchmarks inner loop acquisition functions
using LinearAlgebra
using GaussianDistributions
using SARA
using SARA: stripe_uncertainty_sampling, TemperatureProfile, stripe_to_global, get_relevant_T
using ForwardDiff
using Base.Threads: @threads
using SARA: oSARA

# TODO: throw out sampled data!
############### generate sampling paths with different policies
function learning_path(policy, G::Gaussian, x::AbstractMatrix, y::AbstractMatrix,
                       σ²::AbstractMatrix, stripe_conditions, nsample::Int)
    n = 512
    T = range(200, 1200, length = n + 1) # plotting ranges
    τ = range(2.5, 4, length = n)
    gradients = zeros(n+1, n, nsample)
    C = G
    indices = zeros(Int, nsample)
    chosen_conditions = zeros(2, nsample)
    helper(x, i) = @views vec(x[:, i]) # selects the columns indexed by i, and vectorizes the result
    for i in 1:nsample
        println("stripe $i")
        xj = policy(C, stripe_conditions)
        chosen_conditions[:, i] .= xj
        j = findfirst(==(xj), stripe_conditions)
        indices[i] = j
        indices_i = @view indices[1:i]
        sort!(indices_i) # technically not necessary
        xi, yi, σ²i = helper(x, indices_i), helper(y, indices_i), helper(σ², indices_i)
        C = conditional(G, xi, Gaussian(yi, σ²i), tol = 1e-6)
        μ = mean(C, vec(x))
        δ = vec(y) - μ
        gradients[:, :, i] .= mean(C).(tuple.(T, τ'))
    end
    return T, τ, gradients, chosen_conditions
end

function outer_optical_learning(x::AbstractMatrix, y::AbstractMatrix, σ²::AbstractMatrix,
                                T_offset, stripe_conditions)
    # subsample observations per stripe
    stride = 2 # lower stride important to not get "waves"
    samplesperstripe = 16
    subsample(x::AbstractMatrix) = @view x[end+1-stride*samplesperstripe:stride:end, :]
    subsample(x::AbstractVector) = @view x[end+1-stride*samplesperstripe:stride:end]
    x = subsample(x)
    y = subsample(y)
    σ² = subsample(σ²)

    # standardize global data
    y_normalization = std(y)
    y /= y_normalization
    σ² /= y_normalization^2 # noise variance

    # experimental parameters
    nsample = 32 # number of active learning samples

    # set up relevant T function
    constant_offset_bool = true
    nout = 32
    constant_offset = Val(constant_offset_bool)
    T_range = range(T_offset[1], T_offset[2], length = nout)
    T_range = subsample(T_range)
    T_offset = [maximum(T_range), minimum(T_range)]
    c_min, c_max = T_offset
    nout = length(T_range)
    relevant_T = get_relevant_T(constant_offset, c_min, c_max, nout)

    # stripe uncertainty policy
    policy = stripe_uncertainty_sampling(relevant_T)

    T_length = 30.
    τ_length = .4
    k_scale = 1/4
    σ² ./= k_scale^2
    k = oSARA([T_length, τ_length])
    G = Gaussian(k)
    temperatures, dwelltimes, gradients, chosen_conditions =
                learning_path(policy, G, x, y, σ², stripe_conditions, nsample)
    if false
        f = h5open("outer_optical_learning.h5", "w")
        f["gradients"] = gradients
        f["nsample"] = nsample
        f["kernel"] = ["oSARA"]
        f["temperatures"] = collect(temperatures)
        f["dwelltimes"] = collect(dwelltimes)
        f["temperature lengthscale"] = T_length
        f["dwelltime lengthscale"] = τ_length
        f["y_normalization"] = y_normalization
        f["conditions"] = conditions
        f["dimensions"] = ["temperature", "dwelltime", "nstripes"]
        f["stride"] = stride
        f["samples per stripe"] = samplesperstripe
        close(f)
    end
    return temperatures, dwelltimes, gradients, chosen_conditions, stripe_conditions
end

# load outer loop data
using HDF5
path = "SARA/NatCom2020/outer/data/"
# path = "NatCom2020/outer/"
# file = "Bi2O3_19F44_01_outer_loop_data.h5" # without input noise
# file = "Bi2O3_19F44_01_outer_loop_input_noise_data.h5"
file = "Bi2O3_19F44_01_outer_optical_gradients.h5"
f = h5open(path * file, "r")
dwelltimes = read(f, "dwelltimes")
temperatures = read(f, "temperatures")
gradients = read(f, "gradients")
var_gradients = read(f, "var(gradients)")
T_offset = read(f, "T_offset")
close(f)

# IDEA: throw out data below τ = 2.5?
# ind = dwelltimes .≥ 2.5
# dwelltimes = dwelltimes[ind]
# temperatures = temperatures[ind]

# normalize global data
y_normalization = std(gradients)
gradients ./= y_normalization
var_gradients ./= y_normalization^2 # noise variance

x = tuple.(temperatures, dwelltimes)
y = gradients
σ² = var_gradients

temperature_correction(x::NTuple{2}) = (x[1] + T_offset[2], x[2])
stripe_conditions = temperature_correction.(x[end, :])

temperatures, dwelltimes, gradient_evolution, chosen_conditions = outer_optical_learning(x, y, σ², T_offset, stripe_conditions)

# using Plots
# plotly()
# heatmap(temperatures, dwelltimes, gradient_evolution[:, :, end]')
# scatter!(conditions[1, :], conditions[2, :])


text:
function dpDirect(shot::Shot, v::Float64, src::Source, dx::Float64, dt::Float64)
    nr = length(shot.irz)
    nt = size(shot.d, 1 )
    isz = src.iz; isx = src.ix;
    lsrc= length(src.p)
    damp= hanning(11)[2:6]
    if shot.irz[1] != isz
       error("source and receivers are not on same level")
    end
    dp  = ones(nt, nr)
    for ir = 1 : nr
        l  = abs((shot.irx[ir]-isx)*dx)
        it = ceil(Int64,l/(v*dt))+1+lsrc*2
        dp[1:it, ir] = zeros(it)
        dp[it+1:it+5, ir] = damp
    end
    return dp
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
using Plots

# hyperbola.png --------------------------------------------------------------------------

paramlen=1.5
t0=0.6*paramlen
ran = LinRange(-paramlen, paramlen, 100)

# Hyperbola
plot(framestyle=:origin, aspect_ratio=:equal)
plot!(cosh, sinh, ran, color=:black, lab="")
plot!(t -> -cosh(t), sinh, ran, color=:black, lab="")

# Guides
plot!(t -> t, t -> t, -sinh(last(ran)):0.01:sinh(last(ran)), style=:dash, color=:red, lab="")
plot!(t -> t, t -> -t, -sinh(last(ran)):0.01:sinh(last(ran)), style=:dash, color=:red, lab="")
plot!([(0, sinh(t0)), (cosh(t0), sinh(t0))], style=:dot, color=:black, lab="")
plot!([(cosh(t0), 0), (cosh(t0), sinh(t0))], style=:dot, color=:black, lab="")

# Annotations
annotate!(cosh(t0)+0.1, sinh(t0)/2, Plots.text("sinh t", :left, 8))
annotate!(cosh(t0)/2, sinh(t0)+0.1, Plots.text("cosh t", :up, 8))
savefig("hyperbola.png")

# hyperbolic-paraboloid.png -------------------------------------------------------------

x = LinRange(-10, 10, 100)
y = LinRange(-10, 10, 100)
z(x, y) = x^2-y^2
plot(x, y, z, st=:surface)
savefig("hyperbolic-paraboloid.png")

# image-parametrization-Y.png ------------------------------------------------------------

x = LinRange(-10, 10, 100)
y = LinRange(-10, 10, 100)
z(x, y) = abs(y) <= abs(x) ? x^2-y^2 : NaN
plot(x, y, z, st=:surface)
savefig("image-parametrization-Y.png")

max_stars_repo_name: useredsa/exercises-surfaces-geometry


text:
##
##!------------------------------------------------------------------------------------------------!
##  digraph.jl
##
##  Tests for digraph.jl
##
##  Creation date: 18/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##


text:
module TestQuery
    using TablesDemo
    using Base.Test

    a = [1, 2, 3]
    b = [4, 5, 6]
    c = ["a", "b", "c"]
    d = [:a, :b, :c]
    tbl1 = Table(
        a = NullableArray(a),
        b = NullableArray(b),
        c = NullableArray(c),
        d = NullableArray(d)
    )
    _tbl1 = copy(tbl1)

    # test select
    for (fld, col) in eachcol(tbl1)
        @eval qrya = @query select(tbl1, $fld)
        @eval qryb = @query tbl1 |> select($fld)
        tbl2a = collect(qrya)
        tbl2b = collect(qryb)
        @test isequal(tbl1, _tbl1)
        @test isequal(
            tbl2a,
            @eval Table($fld = $col)
        )
        @test isequal(tbl2a, tbl2b)
    end

    f = [ a[i] * b[i] for i in eachindex(a) ]
    tbl2 = Table(f = NullableArray(f))
    qrya = @query select(tbl1, f = a * b)
    qryb = @query tbl1 |> select(f = a * b)
    tbl2a = collect(qrya)
    tbl2b = collect(qryb)
    @test isequal(tbl1, _tbl1)
    @test isequal(tbl2, tbl2a)
    @test isequal(tbl2a, tbl2b)

    # test filter
    tbl2 = Table(
        a = NullableArray([3]),
        b = NullableArray([6]),
        c = NullableArray(["c"]),
        d = NullableArray([:c])
    )
    qrya = @query filter(tbl1, a > 2)
    qryb = @query tbl1 |> filter(a > 2)
    tbl2a = collect(qrya)
    tbl2b = collect(qryb)
    @test isequal(tbl1, _tbl1)
    @test isequal(tbl2, tbl2a)
    @test isequal(tbl2a, tbl2b)

    # test summarize
    tbl2 = Table(
        b_avg = NullableArray([mean([4, 5, 6])])
    )
    qrya = @query summarize(tbl1, b_avg = mean(b))
    qryb = @query tbl1 |> summarize(b_avg = mean(b))
    tbl2a = collect(qrya)
    tbl2b = collect(qryb)
    @test isequal(tbl1, _tbl1)
    @test isequal(tbl2, tbl2a)
    @test isequal(tbl2a, tbl2b)

    # test combinations
    tbl2 = Table(
        b = NullableArray([4])
    )
    qrya = @query select(filter(tbl1, a == 1), b)
    qryb = @query tbl1 |>
        filter(a == 1) |>
        select(b)
    tbl2a = collect(qrya)
    tbl2b = collect(qryb)
    @test isequal(tbl1, _tbl1)
    @test isequal(tbl2a, tbl2)
    @test isequal(tbl2a, tbl2b)
end


text:
const new_package_title_regex = r"^New package: (\w*?) v(\S*?)$"

const new_version_title_regex = r"^New version: (\w*?) v(\S*?)$"

const commit_regex = r"(?:^|\n|\r\n)(?:\-|\*) Commit: (\w*?)(?:$|\n|\r\n)"

is_new_package(pull_request::GitHub.PullRequest) = occursin(new_package_title_regex, title(pull_request))

is_new_version(pull_request::GitHub.PullRequest) = occursin(new_version_title_regex, title(pull_request))

function check_authorization(pkg, pr_author_login, authorized_authors,
                             authorized_authors_special_jll_exceptions,
                             error_exit_if_automerge_not_applicable)
    if pr_author_login ∉ vcat(authorized_authors,
                              authorized_authors_special_jll_exceptions)
        throw_not_automerge_applicable(
            AutoMergeAuthorNotAuthorized,
            "Author $(pr_author_login) is not authorized to automerge. Exiting...";
            error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable
        )
        return :not_authorized
    end

    # A JLL-only author (e.g. `jlbuild`) is not allowed to register
    # non-JLL packages.
    this_is_jll_package = is_jll_name(pkg)
    if (!this_is_jll_package && pr_author_login ∉ authorized_authors)
        throw_not_automerge_applicable(
            AutoMergeAuthorNotAuthorized,
            "This package is not a JLL package. Author $(pr_author_login) is not authorized to register non-JLL packages. Exiting...";
            error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable
        )
        return :not_authorized
    end

    if pr_author_login ∈ authorized_authors_special_jll_exceptions
        return :jll
    end

    return :normal
end

function parse_pull_request_title(::NewVersion,
                                  pull_request::GitHub.PullRequest)
    m = match(new_version_title_regex, title(pull_request))
    pkg = convert(String, m.captures[1])::String
    version = VersionNumber(m.captures[2])
    return pkg, version
end

function commit_from_pull_request_body(pull_request::GitHub.PullRequest)
    pr_body = body(pull_request)
    m = match(commit_regex, pr_body)
    return convert(String, m.captures[1])::String
end

function parse_pull_request_title(::NewPackage,
                                  pull_request::GitHub.PullRequest)
    m = match(new_package_title_regex, title(pull_request))
    pkg = convert(String, m.captures[1])::String
    version = VersionNumber(m.captures[2])
    return pkg, version
end

function pull_request_build(api::GitHub.GitHubAPI,
                            pr_number::Integer,
                            current_pr_head_commit_sha::String,
                            registry::GitHub.Repo,
                            registry_head::String;
                            whoami::String,
                            auth::GitHub.Authorization,
                            authorized_authors::Vector{String},
                            authorized_authors_special_jll_exceptions::Vector{String},
                            error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable,
                            master_branch::String,
                            master_branch_is_default_branch::Bool,
                            suggest_onepointzero::Bool,
                            registry_deps::Vector{<:AbstractString} = String[],
                            check_license::Bool,
                            public_registries::Vector{<:AbstractString} = String[],
                            read_only::Bool)::Nothing
    pr = my_retry(() -> GitHub.pull_request(api, registry, pr_number; auth=auth))
    _github_api_pr_head_commit_sha = pull_request_head_sha(pr)
    if current_pr_head_commit_sha != _github_api_pr_head_commit_sha
        throw(AutoMergeShaMismatch("Current commit sha (\"$(current_pr_head_commit_sha)\") does not match what the GitHub API tells us (\"$(_github_api_pr_head_commit_sha)\")"))
    end
    result = pull_request_build(api,
                                pr,
                                current_pr_head_commit_sha,
                                registry,
                                registry_head;
                                auth=auth,
                                authorized_authors=authorized_authors,
                                authorized_authors_special_jll_exceptions=authorized_authors_special_jll_exceptions,
                                error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable,
                                master_branch=master_branch,
                                master_branch_is_default_branch=master_branch_is_default_branch,
                                suggest_onepointzero=suggest_onepointzero,
                                whoami=whoami,
                                registry_deps=registry_deps,
                                check_license=check_license,
                                public_registries=public_registries,
                                read_only=read_only)
    return result
end

function pull_request_build(api::GitHub.GitHubAPI,
                            pr::GitHub.PullRequest,
                            current_pr_head_commit_sha::String,
                            registry::GitHub.Repo,
                            registry_head::String;
                            auth::GitHub.Authorization,
                            authorized_authors::Vector{String},
                            authorized_authors_special_jll_exceptions::Vector{String},
                            error_exit_if_automerge_not_applicable::Bool,
                            master_branch::String,
                            master_branch_is_default_branch::Bool,
                            suggest_onepointzero::Bool,
                            whoami::String,
                            registry_deps::Vector{<:AbstractString} = String[],
                            check_license::Bool,
                            public_registries::Vector{<:AbstractString} = String[],
                            read_only::Bool)::Nothing
    # 1. Check if the PR is open, if not quit.
    # 2. Determine if it is a new package or new version of an
    #    existing package, if neither quit.
    # 3. Check if the author is authorized, if not quit.
    # 4. Call the appropriate method for new package or new version.
    if !is_open(pr)
        throw_not_automerge_applicable(
            AutoMergePullRequestNotOpen,
            "The pull request is not open. Exiting...";
            error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable
        )
        return nothing
    end

    if is_new_package(pr)
        registration_type = NewPackage()
    elseif is_new_version(pr)
        registration_type = NewVersion()
    else
        throw_not_automerge_applicable(
            AutoMergeNeitherNewPackageNorNewVersion,
            "Neither a new package nor a new version. Exiting...";
            error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable,
        )
        return nothing
    end

    pkg, version = parse_pull_request_title(registration_type, pr)
    pr_author_login = author_login(pr)
    authorization = check_authorization(pkg, pr_author_login,
                                        authorized_authors,
                                        authorized_authors_special_jll_exceptions,
                                        error_exit_if_automerge_not_applicable)

    if authorization == :not_authorized
        return nothing
    end

    registry_master = clone_repo(registry)
    if !master_branch_is_default_branch
        checkout_branch(registry_master, master_branch)
    end
    data = GitHubAutoMergeData(;api = api,
                               registration_type = registration_type,
                               pr = pr,
                               pkg = pkg,
                               version = version,
                               current_pr_head_commit_sha = current_pr_head_commit_sha,
                               registry = registry,
                               auth = auth,
                               authorization = authorization,
                               registry_head = registry_head,
                               registry_master = registry_master,
                               suggest_onepointzero = suggest_onepointzero,
                               whoami = whoami,
                               registry_deps = registry_deps,
                               public_registries = public_registries,
                               read_only = read_only)
    pull_request_build(data, registration_type; check_license=check_license)
    rm(registry_master; force = true, recursive = true)
    return nothing
end


text:
using Revise
using SMLMData
using PlotlyJS

filepath="Y:\\Projects\\BaGoL\\21_11_02_DNAPAINT_EGFR\\activated_EGF\\Results\\Cell_05\\Label_01\\Data_2021-11-4-6-11-39"
filename="Data_2021-11-4-6-11-39_Results.mat"
smd=SMLMData.SMITEsmd(filepath,filename)

smld=SMLMData.SMLD2D(smd)

plot(scattergl(x=smld.x, y=smld.y, mode="markers",markersize = .1, linewidth=.1))


text:
function get_rop_sections()
    return ["MODIFICATION CODE", "BUS VOLTAGE", "ADJUSTABLE BUS SHUNTS",
            "BUS LOADS", "ADJUSTABLE BUS LOAD", "GENERATOR DISPATCH",
            "ACTIVE POWER DISPATCH", "GENERATOR RESERVE",
            "GENERATION REACTIVE CAP", "ADJUSTABLE BRANCH",
            "PIECEWISE LINEAR COST", "PIECEWISE QUADRATIC COST",
            "POLYNOMIAL EXPONENTIAL COST", "PERIOD RESERVES",
            "BRANCH FLOWS", "INTERFACE FLOWS", "LINEAR CONSTRAINT DEP"]
end

function get_rop_section_info()
    generator_dispatch = [["BUS", 1, Int], ["GENID", 2, String],
                          ["DSPTBL", 4, Int]]
    active_power_dispatch = [["TBL", 1, Int], ["CTBL", 7, Int]]
    piecewise_linear_cost = [["LTBL", 1, Int], ["NPARIS", 3, Int]]

    section_info = Dict{String,Array}(
        "MODIFICATION CODE" => [],
        "BUS VOLTAGE" => [],
        "ADJUSTABLE BUS SHUNTS" => [],
        "BUS LOADS" => [],
        "ADJUSTABLE BUS LOAD" => [],
        "GENERATOR DISPATCH" => generator_dispatch,
        "ACTIVE POWER DISPATCH" => active_power_dispatch,
        "GENERATOR RESERVE" => [],
        "GENERATION REACTIVE CAP" => [],
        "ADJUSTABLE BRANCH" => [],
        "PIECEWISE LINEAR COST" => piecewise_linear_cost,
        "PIECEWISE QUADRATIC COST" => [],
        "POLYNOMIAL EXPONENTIAL COST" => [],
        "PERIOD RESERVES" => [],
        "BRANCH FLOWS" => [],
        "INTERFACE FLOWS" => [],
        "LINEAR CONSTRAINT DEP" => [])

    return section_info
end

function idx_rop()
    ROP_BUS = 1
    ROP_ID = 2
    ROP_NPAIRS = 3
    ROP_PAIRSTART = 4

    return ROP_BUS, ROP_ID, ROP_NPAIRS, ROP_PAIRSTART
end

function rop_join_tables(ropdata::Dict{String,Array})
    gendspt = ropdata["GENERATOR DISPATCH"]
    acpdspt = ropdata["ACTIVE POWER DISPATCH"]
    pwlcost = ropdata["PIECEWISE LINEAR COST"]

    ndata = size(gendspt, 1)
    data = Array{Any}(undef, ndata, 4)

    for i=1:ndata
        dsptbl = gendspt[i, 3]

        if dsptbl <= ndata && dsptbl == acpdspt[dsptbl, 1]
            acp_row = dsptbl
        else
            acp_row = findfirst(dsptbl .== acpdspt[:, 1])
        end

        if acpdspt[acp_row, 2] == pwlcost[acp_row, 1]
            pwl_row = acp_row
        else
            pwl_row = findfirst(acpdspt[acp_row, 2] .== pwlcost[:, 1])
        end

        data[i, 1] = gendspt[i, 1]
        data[i, 2] = gendspt[i, 2]
        data[i, 3] = pwlcost[pwl_row, 2]
        data[i, 4] = pwlcost[pwl_row, 3]
    end

    return data
end

function rop_pwl(lines::Array{String}, pos::Int, info::Array;
                 delim=",")
    NPAIRS_LOC = 2
    start = last = pos
    ndata = npairs = tot_npairs = 0

    line = strip(lines[last])
    while line[1] ∉ [ '0', 'Q' ]
        fields = split(line, delim)
        npairs = parse(info[2][3], fields[info[2][2]])
        tot_npairs += npairs
        ndata += 1
        last += 1 + npairs
        line = strip(lines[last])
    end

    data = Array{Any}(undef, ndata, length(info)+1)
    pairs = Array{Float64}(undef, tot_npairs, 2)

    offset = start
    npairs = tot_npairs = 0
    for i=1:ndata
        fields = split(lines[offset], delim)

        for j=1:length(info)
            data[i,j] = parse(info[j][3], fields[info[j][2]])
        end

        offset += 1
        npairs = data[i,NPAIRS_LOC]
        data[i,end] = tot_npairs
        for j=1:npairs
            tot_npairs += 1
            fields = split(lines[offset], delim)
            pairs[tot_npairs,1] = parse(Float64, fields[1])
            pairs[tot_npairs,2] = parse(Float64, fields[2])
            offset += 1
        end
    end

    pos = last + 1
    return pos, data, pairs
end

function parse_rop(filename::AbstractString; delim=",")
    f = open(filename, "r")
    lines = readlines(f)
    close(f)

    sections = get_rop_sections()
    section_info = get_rop_section_info()
    pos = 1

    ropdata = Dict{String,Array}()
    for s in sections
        info = section_info[s]

        if isempty(info)
            pos = parse_skipsection(lines, pos) + 1
            continue
        end

        if s == "PIECEWISE LINEAR COST"
            pos, data, pairs = rop_pwl(lines, pos, info; delim=delim)
            ropdata["PAIRS"] = pairs
        else
            pos, data = parse_data(lines, pos, info; delim=delim)
        end

        ropdata[s] = data
    end

    # Perform join operation
    if haskey(ropdata, "PIECEWISE LINEAR COST")
        data = rop_join_tables(ropdata)
        ropdata["JOINED GENERATOR DISPATCH"] = data
    end

    return ropdata
end


text:
# Benchmark (Performance) a Solver

using BenchmarkTools;

include("GenerateQuadraticProgram.jl");

# const BENCHMARK_NUM_SAMPLES = 10;
# const BENCHMARK_NUM_EVALS   = 1;
# const BENCHMARK_RUN_TIME    = 90;

function GenerateElementsVector(numElementsMin = 20, numElementsMax = 2000; logSpace = true)
    if (logSpace)
        # Log Scale: More data points in the high values
        vNumElements = round.(Int, 10 .^ (LinRange(log10(numElementsMin), log10(numElementsMax), numDims)));
    else
        vNumElements = round.(Int, LinRange(numElementsMin, numElementsMax, numDims));
    end

    return vNumElements;
end

function BenchmarkSolver(solverFun, problemClass::ProblemClass, vNumElements, vNumConstraints; benchMarkSamples = 10, benchMarkEvals = 1, benchMarkSeconds = 120)

    # mR = Matrix{Float64}(undef, numDims, 5); #<! Min Time [Nano Sec], Median Time [Nano Sec], Max Time [Nano Sec], Number of Allocations, Allocations Size
    cBenchMark = Vector{Any}(undef, numDims);

    for ii in 1:numDims
        mP, vQ, mA, vL, vU = GenerateRandomQP(problemClass, vNumElements[ii], numConstraints = vNumConstraints[ii]);
        # Some problems might change the actual data
        numElements     = size(mP, 1);
        numConstraints  = size(mA, 1);

        vXX = zeros(numElements);

        sBenchMark = @benchmarkable $solverFun(vX, $mP, $vQ, $mA, $vL, $vU) setup = (vX = copy($vXX)); #<! We must interpolate the function as it is not in the global scope
        cBenchMark[ii] = run(sBenchMark, samples = benchMarkSamples, evals = benchMarkEvals, seconds = benchMarkSeconds);
    end

    return cBenchMark;

end

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
using RobotZoo
using StaticArrays, LinearAlgebra

model = RobotZoo.Cartpole()
n,m = RD.dims(model)
N = 11
x,u = rand(model)
t,dt = 1.1,0.1
z = KnotPoint(x,u,t,dt)

#--- Generate some constraints
# Circle Constraint
xc = SA[1,1,1]
yc = SA[1,2,3]
r  = SA[1,1,1]
cir = CircleConstraint(n, xc, yc, r)

# Goal Constraint
xf = @SVector rand(n)
goal = GoalConstraint(xf)

# Linear Constraint
p = 5
A = @SMatrix rand(p,n+m)
b = @SVector rand(p)
lin = LinearConstraint(n,m,A,b, Inequality())

# Bound Constraint
xmin = -@SVector rand(n)
xmax = +@SVector rand(n)
umin = -@SVector rand(m)
umax = +@SVector rand(m)
bnd = BoundConstraint(n,m, x_min=xmin, x_max=xmax, u_min=umin, u_max=umax)

# Dynamics Constraint
dmodel = RD.DiscretizedDynamics{RD.RK4}(model)

##--- Create the Constraint List
cons = ConstraintList(n,m,N)
add_constraint!(cons, cir, 1:N)
add_constraint!(cons, goal, N)
add_constraint!(cons, lin, 2:N-1)
add_constraint!(cons, bnd, 1:N-1)

prob, = Problems.Cartpole(N=N)

##--- Create an ALConstraintSet
ilqr = Altro.iLQRSolver(prob, use_static=Val(false))
conset = Altro.ALConstraintSet{Float64}()
Altro.initialize!(conset, cons, ilqr.Z, ilqr.opts, ilqr.obj.J, ilqr.Efull)
Altro.reset!(conset)

cval = conset[1]
@test cval isa Altro.ALConstraint{Float64, typeof(cir)}
@test cval.inds == cons.inds[1]
@test size(cval.jac[1]) == (RD.output_dim(cir), n)
@test size(cval.jac) == (N,)
cval.jac[1] .= 1
@test cval.jac[1] == ones(RD.output_dim(cir), n)

cval = conset[3]
@test cval isa Altro.ALConstraint{Float64,typeof(lin)}
@test cval.inds == cons.inds[3]
@test size(cval.jac[1]) == (RD.output_dim(lin), n+m)
@test size(cval.jac) == (N-2,)
cval.jac[2] .= 1
cval.jac[2][:,n+1:end] .= 2
@test cval.jac[2] ≈ [ones(RD.output_dim(lin), n) 2*ones(RD.output_dim(lin),m)]

# Test iteration
@test all([con.con for con in conset] .=== [con for con in cons])
@test length(conset) == length(cons)

##--- Test evaluation
Z = SampledTrajectory{n,m}([rand(n) for k = 1:N], [rand(m) for k = 1:N-1], dt=fill(dt,N-1))
Altro.evaluate_constraints!(conset, Z)
@test conset[3].vals[2] ≈ RD.evaluate(lin, Z[3])
@test conset[2].vals[1] ≈ RD.evaluate(goal, Z[end])

Altro.constraint_jacobians!(conset, Z)
@test conset[2].jac[1] ≈ I(n)
∇c = TO.gen_jacobian(lin)
c = zeros(RD.output_dim(lin))
RD.jacobian!(lin, ∇c, c, Z[5])
@test conset[3].jac[5] ≈ ∇c ≈ A

##--- Stats functions
TO.evaluate_constraints!(conset)
vals = [conval.vals for conval in conset]
viols = map(enumerate(vals)) do (i,val)
    pos(x) = x > 0 ? x : zero(x)
    if TO.sense(cons[i]) == Inequality()
        [pos.(v) for v in val]
    else
        [abs.(v) for v in val]
    end
end
pos(x) = max(0,x)
c_max = map(enumerate(vals)) do (i,val)
    if TO.sense(cons[i]) == Inequality()
        maximum(maximum.([pos.(c) for c in val]))
    else
        maximum(norm.(val,Inf))
    end
end

@test maximum(c_max) ≈ max_violation(conset)
@test maximum(maximum.([maximum.(viol) for viol in viols])) ≈ max_violation(conset)

p = 2
@test norm([norm(norm.(val,p),p) for val in viols],p) ≈ Altro.normviolation!(conset, p)
p = 1
@test norm([norm(norm.(val,p),p) for val in viols],p) ≈ Altro.normviolation!(conset, p)
p = Inf
@test norm([norm(norm.(val,p),p) for val in viols],p) ≈ Altro.normviolation!(conset, p)
@test !(Altro.normviolation!(conset, Inf) ≈ Altro.normviolation!(conset, 2))

# TODO: test penalty reset
conset[1].opts.usedefault[:penalty_initial] = true
conset[1].opts.penalty_initial = 10
Altro.resetparams!(conset[1])
Altro.reset_penalties!(conset)
@test conset[1].opts.penalty_initial == ilqr.opts.penalty_initial
@test Altro.max_penalty(conset) ≈ 1.0

Altro.setparams!(conset[1], penalty_initial = 10.0)
Altro.resetparams!(conset[1])
@test conset[1].opts.penalty_initial == 10.0

Altro.reset_penalties!(conset)
@test Altro.max_penalty(conset) ≈ 10.0

conset[3].μ[6] = @SVector fill(30, RD.output_dim(cons[3]))
@test Altro.max_penalty(conset) == 30
Altro.reset_penalties!(conset)
@test Altro.max_penalty(conset) ≈ 10.0

Altro.evaluate_constraints!(conset, Z)
conset[2].vals[1][2] = 2*max_violation(conset)

@test Altro.findmax_violation(conset) == "GoalConstraint at time step 11 at index 2"
Altro.findmax_violation(conset)
conset[4].vals[2][3] = 2*max_violation(conset)
@test Altro.findmax_violation(conset) == "BoundConstraint at time step 2 at x max 3"

#--- AL Updates
Altro.reset_duals!(conset)
Altro.dualupdate!(conset)
@test conset[2].λ[1] ≈ conset[2].vals[1]
@test conset[1].λ[3] ≈ clamp.(10*conset[1].vals[3], 0, Inf)
λ0 = copy(conset[1].λ[3])

Altro.dualupdate!(conset)
@test conset[1].λ[3] ≈ clamp.(λ0 .+ 10*conset[1].vals[3], 0, Inf)
Altro.reset_duals!(conset)
@test conset[1].λ[3] == zeros(RD.output_dim(cons[1]))

Altro.reset_penalties!(conset)
Altro.penaltyupdate!(conset)
@test Altro.max_penalty(conset) ≈ 100
@test conset[2].μ[1] ≈ fill(10, RD.output_dim(cons[2]))


text:
using PureEOS
using Base.Test

# write your own tests here
println("Viscosity test")
@test waterVisc(298.15,998/18.015268) ≈ 889.735100 atol=0.001
@test waterVisc(298.15,1200/18.015268) ≈ 1437.649467 atol=0.001
@test waterVisc(373.15,1000/18.015268) ≈ 307.883622 atol=0.001
@test waterVisc(433.15,1/18.015268) ≈ 14.538324 atol=0.001
@test waterVisc(433.15,1000/18.015268) ≈ 217.685358 atol=0.001
@test waterVisc(873.15,1/18.015268) ≈ 32.619287 atol=0.001
@test waterVisc(873.15,100/18.015268) ≈ 35.802262 atol=0.001
@test waterVisc(873.15,600/18.015268) ≈ 77.430195 atol=0.001
@test waterVisc(1173.15,1/18.015268) ≈ 44.217245 atol=0.001
@test waterVisc(1173.15,100/18.015268) ≈ 47.640433 atol=0.001
@test waterVisc(1173.15,400/18.015268) ≈ 64.154608 atol=0.001
println("thermalconductivity test")
@test waterTherm(647.35,1/18.015268) ≈ 51.9298924 atol=0.001
@test waterTherm(647.35,122/18.015268) ≈ 130.922885 atol=0.1
@test waterTherm(647.35,222/18.015268) ≈ 367.787459 atol=20
@test waterTherm(647.35,272/18.015268) ≈ 757.959776 atol=40
@test waterTherm(647.35,322/18.015268) ≈ 1443.75556 atol=200
@test waterTherm(647.35,372/18.015268) ≈ 650.319402 atol=40
@test waterTherm(647.35,422/18.015268) ≈ 448.883487 atol=20
@test waterTherm(647.35,750/18.015268) ≈ 600.961346 atol=0.001

println("ethanol thermalconductivity test")
#@test ethanolTherm(300.0,850.0/46.0) ≈ 209.68 atol=2
#@test ethanolTherm(400.0,2.0/46.0) ≈ 26.108 atol=1
#@test ethanolTherm(400.0,690.0/46.0) ≈ 149.21 atol=1.5
#@test ethanolTherm(500.0,10.0/46.0) ≈ 39.594 atol=1
#@test ethanolTherm(500.0,10.0/46.0) ≈ 40.755 atol=1

#Data from NIST


text:
#*************************************************************************#
#
# Predictor type
#
#*************************************************************************#
"""
    Predictor

Type used to make a prediction about a temporal series.
"""
struct Predictor
    # Range of values of interest.
    start::Int
    length::Int
    # Internal Dynamic Model of the predictor.
    dynmdl::DifferenceEquation
    # Input to `dynmdl`.
    u::Matrix{Float64}
    # Initial conditions used by the `dynmdl`.
    y0::Matrix{Float64}
end

#*************************************************************************#
#
# Usefull constructors
#
#*************************************************************************#
"""
    one_step_ahead(mdl, yterms, uterms, y, u) -> pred::Predictor

One-step-ahead predictor for the given data set.
"""
function one_step_ahead{R<:Real,
                        S<:Real}(mdl::ParametricModel,
                                 yterms::Vector{Vector{Int}},
                                 uterms::Vector{Vector{Int}},
                                 y::Matrix{R},
                                 u::Matrix{S})
    # Get dynmdl
    dynmdl = DifferenceEquation(mdl, uterms=[yterms; uterms])

    # Get info about dynmdl
    iodelay = dynmdl.iodelay
    m = dynmdl.m
    n = dynmdl.n
    ny, nd = size(y)

    # Data for Simulation
    u_pred = [y[:, 1:nd-iodelay]; u[:, 1:nd-iodelay]]
    y0 = Matrix{Float64}(0, 0)

    Predictor(iodelay+m, nd-iodelay-m+1, dynmdl, u_pred, y0)
end

"""
    one_step_ahead(mdl, yterms, uterms, iddata) -> pred::Predictor

One-step-ahead predictor for the given data set.
"""
function one_step_ahead(mdl::ParametricModel,
                        yterms::Vector{Vector{Int}},
                        uterms::Vector{Vector{Int}},
                        iddata::IdData)
    one_step_ahead(mdl, yterms, uterms, iddata.y, iddata.u)
end

"""
    free_run_simulation(mdl, yterms, uterms, y, u) -> pred::Predictor

Free-run simulation for the given data set.
"""
function free_run_simulation(mdl::ParametricModel,
                             yterms::Vector{Vector{Int}},
                             uterms::Vector{Vector{Int}},
                             y::Matrix{<:Real},
                             u::Matrix{<:Real})
    # Get optimal predictor
    dynmdl = DifferenceEquation(mdl, yterms=yterms, uterms=uterms)

    # Get info about optimal predictor
    iodelay = dynmdl.iodelay
    m = dynmdl.m
    n = dynmdl.n
    n0 = max(iodelay+m, n+1)
    ny, nd = size(y)

    # Data for computing error
    u_pred = u[:, n0-iodelay-m+1:nd-iodelay]
    y0 = y[:, n0-n:n0-1]

    Predictor(n0, nd-n0+1, dynmdl, u_pred, y0)
end

"""
    free_run_simulation(mdl, yterms, uterms, iddata) -> pred::Predictor

Free-run simulation for the given data set.
"""
function free_run_simulation(mdl::ParametricModel,
                             yterms::Vector{Vector{Int}},
                             uterms::Vector{Vector{Int}},
                             iddata::IdData)
    free_run_simulation(mdl, yterms, uterms, iddata.y, iddata.u)
end

#*************************************************************************#
#
# Generate Buffers
#
#*************************************************************************#
"""
    y_buffer(pred::Predictor) -> y

Return the array to be used as buffer for the predicted output sequence `y`.
"""
function y_buffer(pred::Predictor)
    return Matrix{Float64}(pred.dynmdl.ny, pred.length)
end

"""
    dΘ_buffer(pred::Predictor) -> dΘ

Return the array to be used as buffer for the array of derivatives `dΘ`.
"""
function dΘ_buffer(pred::Predictor)
    return Array{Float64}(pred.dynmdl.ny, pred.length,
                          length_params(pred.dynmdl.mdl))
end

"""
    dy0_buffer(pred::Predictor) -> dy0

Return the array to be used as buffer for the array of derivatives `dy0`.
"""
function dy0_buffer(pred::Predictor)
    return Array{Float64}(pred.dynmdl.ny, pred.length,
                          pred.dynmdl.ny, pred.dynmdl.n)
end

#*************************************************************************#
#
# predict!
#
#*************************************************************************#
function predict!(pred::Predictor,
                  y::AbstractMatrix{<:Real},
                  Θ::AbstractVector{<:Real},
                  y0::AbstractMatrix{<:Real}=pred.y0;
                  dΘ::AbstractArray{<:Real, 3}=Array{Float64}(0, 0, 0),
                  dy0::AbstractArray{<:Real, 4}=Array{Float64}(0, 0, 0, 0))

    simulate!(pred.dynmdl, pred.u, Θ, y0, y,  dΘ=dΘ, dy0=dy0)

end

#*************************************************************************#
#
# predict
#
#*************************************************************************#
function predict(pred::Predictor,
                 Θ::AbstractVector{<:Real},
                 y0::AbstractMatrix{<:Real}=pred.y0)

    y = y_buffer(pred)

    predict!(pred, y, Θ, y0)

    return y
end

#*************************************************************************#
#
# get_slice
#
#*************************************************************************#
function get_slice(pred::Predictor,
                   v::AbstractVector{<:Real})
    return v[pred.start:(pred.start+pred.length-1)]
end

function get_slice(pred::Predictor,
                   v::AbstractMatrix{<:Real})
    return v[:, pred.start:(pred.start+pred.length-1)]
end


text:
#
# Self - Water
#

using ComplexMixtures ; const CM = ComplexMixtures
using PDBTools
using Plots
ENV["GKSwstype"] = "nul"

atoms = readPDB("../../test/data/NAMD/structure.pdb")

# Self correlation, thus the solute and solvent indexes are identical
water_atoms = select(atoms,"water")
water = CM.Selection( water_atoms, natomspermol=3)

# Input options for the calcualtion
options = CM.Options(binstep=0.2)

# Run MDDF calculation, and get the resutls in the R structure
nlabel="lcP"
trajectory = CM.Trajectory("../../test/data/NAMD/trajectory.dcd",water)
@time R = CM.mddf(trajectory,options)

plot(layout=(4,1))

sp=1
plot!(ylabel="MDDF or RDF",subplot=sp)
plot!(R.d,R.mddf,subplot=sp,label="mddf")
plot!(R.d,R.rdf,subplot=sp,label="rdf")
plot!(legend=:topright,subplot=sp)

sp=2
plot!(ylabel="KB",subplot=sp)
plot!(R.d,R.kb,subplot=sp,label="mddf")
plot!(R.d,R.kb_rdf,subplot=sp,label="rdf")
plot!(legend=:topright,subplot=sp)

sp=3
plot!(ylabel="atom contrib",subplot=sp)
for i in 1:water.natomspermol
  plot!(R.d,R.solute_atom[:,i],subplot=sp,label="")
end
y1 = similar(R.d)
for i in 1:size(R.d,1)
  y1[i] = sum(R.solvent_atom[i,:])
end
plot!(R.d,y1,subplot=sp,label="sum")
plot!(legend=:topright,subplot=sp)

sp=4
plot!(ylabel="count",subplot=sp)
plot!(R.d,R.md_count,subplot=sp,label="R md")
plot!(legend=:topleft,subplot=sp)

plot!(size=(600,800))
savefig("./example3.pdf")


text:
@doc raw"""
    Vector2

An object `v` of the type [`Vector2`](@ref) is a vector in the two-dimensional
vector space ``ℝ^2``.

Geometrically spoken, a vector is the difference of two points ``p_2, p_1`` and
and denotes the direction and the distance from ``p_1`` to ``p_2``.

CGAL defines a symbolic constant [`NULL_VECTOR`](@ref). We will explicitly
state where you can pass this constant as an argument instead of a vector
initialized with zeros.
"""
Vector2

"""
    Vector2(a::Point2, b::Point2)

Introduces the vector `b - a`.
"""
Vector2(a::Point2, b::Point2)

"""
    Vector2(s::Segment2)

Introduces the vector `target(s) - source(s)`.
"""
Vector2(s::Segment2)

"""
    Vector2(r::Ray2)

Introduces the vector having the same direction as `r`.
"""
Vector2(r::Ray2)

"""
    Vector2(l::Line2)

Introduces the vector having the same direction as `t`.
"""
Vector2(l::Line2)

"""
    Vector2(::NullVector = NULL_VECTOR)

Introduces a null vector `v`.
"""
Vector2(::NullVector = NULL_VECTOR)

"""
    Vector2(x::Real, y::Real)

Introduces a vector `v` initialized to ``(x, y)``.
"""
Vector2(x::Real, y::Real) = Vector2(convert(FT, x), convert(FT, y))

@doc raw"""
    Vector2(hx::Real, hy::Real, hw::Real = 1)

Introduces a vector `v` initialized to `(hx/hw,hy/hw)`.

!!! info "Precondition"

    `hw` ``≠`` `0`
"""
Vector2(hx::Real, hy::Real, hw::Real) =
    Vector2(convert(RT, hx), convert(RT, hy), convert(RT, hw))

"""
    hx(p::Vector2)

Returns the homogeneous ``x`` coordinate.
"""
hx(p::Vector2)

"""
    hy(p::Vector2)

Returns the homogeneous ``y`` coordinate.
"""
hy(p::Vector2)

"""
    hw(p::Vector2)

Returns the homogenizing coordinate.
"""
hw(p::Vector2)

"""
    x(p::Vector2)

Returns the Cartesian ``x`` coordinate, that is [`hx()`](@ref)/[`hw()`](@ref).
"""
x(p::Vector2)

"""
    y(p::Vector2)

Returns the Cartesian ``y`` coordinate, that is [`hy()`](@ref)/[`hw()`](@ref).
"""
y(p::Vector2)

@doc raw"""
    homogeneous(v::Vector2, i::Integer)

Returns the `i`ᵗʰ homogeneous coordinate of `v`.

!!! info "Precondition"

    ``0 ≤ i ≤ 2``.
"""
homogeneous(v::Vector2, i::Integer)

@doc raw"""
    cartesian(v::Vector2, i::Integer)

Returns the `i`ᵗʰ Cartesian coordinate of `v`.

!!! info "Precondition"

    ``0 ≤ i ≤ 1``.
"""
cartesian(v::Vector2, i::Integer)

"""
    dimension(v::Vector2)

Returns the dimension (the constant 2).
"""
dimension(v::Vector2)

"""
    direction(v::Vector2)

Returns the direction which passes through `v`.
"""
direction(v::Vector2)

"""
    transform(v::Vector2, t::AffTransformation2)

Returns the vector obtained by applying `t` on `v`.
"""
transform(v::Vector2, t::AffTransformation2)

"""
    perpendicular(v::Vector2, o::Orientation)

Returns the vector perpendicular to `v` in clockwise or counterclockwise
orientation.
"""
perpendicular(v::Vector2, o::Orientation)

"""
    squared_length(v::Vector2)

Returns the squared length of `v`.
"""
squared_length(v::Vector2)

"""
    ==(v::Vector2, w::Vector2)

Test for equality: two vectors are equal, iff their ``x`` and ``y`` coordinates
are equal.

You can compare a vector with the `NULL_VECTOR`.
"""
==(v::Vector2, w::Vector2)

"""
    +(v::Vector2, w::Vector2)

Addition.
"""
+(v::Vector2, w::Vector2)

"""
    -(v::Vector2, w::Vector2)

Subtraction.
"""
-(v::Vector2, w::Vector2)

"""
    -(v::Vector2)

Returns the opposite vector.
"""
-(v::Vector2)

"""
    *(v::Vector2, w::Vector2)

Returns the scalar product (= inner product) of the two vectors.
"""
*(v::Vector2, w::Vector2)

"""
    /(v::Vector2, s::Real)

Division by a scalar.
"""
/(v::Vector2, s::Real)

@cxxdereference Base.:/(v::Vector2, s::Real) = v / convert(RT, s)


text:
using CUDA.APIUtils

@testset "call" begin
    # ccall throws if the lib doesn't exist, even if not called
    foo(x) = (x && ccall((:whatever, "nonexisting"), Cvoid, ()); 42)
    if VERSION < v"1.4.0-DEV.653"
        @test_throws ErrorException foo(false)
    else
        foo(false)
    end

    # @runtime_ccall prevents that
    bar(x) = (x && @runtime_ccall((:whatever, "nonexisting"), Cvoid, ()); 42)
    @test bar(false) == 42
    # but should still error nicely if actually calling the library
    @test_throws ErrorException bar(true)

    # ccall also doesn't support non-constant arguments
    lib = Ref("libjulia")
    baz() = ccall((:jl_getpid, lib[]), Cint, ())
    @test_throws TypeError baz()

    # @runtime_ccall supports that
    qux() = @runtime_ccall((:jl_getpid, lib[]), Cint, ())
    @test qux() == getpid()

    # decoding ccall/@runtime_ccall
    @test decode_ccall_function(:(ccall((:fun, :lib)))) == "fun"
    @test decode_ccall_function(:(@runtime_ccall((:fun, :lib)))) == "fun"
end

@testset "@enum_without_prefix" begin
    mod = @eval module $(gensym())
        using CUDA.APIUtils
        @enum MY_ENUM MY_ENUM_VALUE
        @enum_without_prefix MY_ENUM MY_
    end

    @test mod.ENUM_VALUE == mod.MY_ENUM_VALUE
end
@testset "@checked" begin
    mod = @eval module $(gensym())
        using CUDA.APIUtils

        const checks = Ref(0)
        macro check(ex)
            esc(quote
                $checks[] += 1
                $ex
            end)
        end

        @checked function foo()
            ccall(:jl_getpid, Cint, ())
        end
    end

    @test mod.checks[] == 0
    @test mod.foo() == getpid()
    @test mod.checks[] == 1
    @test mod.unsafe_foo() == getpid()
    @test mod.checks[] == 1
end

@testset "@argout" begin
  f() = 1
  f(a) = 2
  f(a,b) = 3

  @test CUDA.@argout(f()) == nothing
  @test CUDA.@argout(f(out(4))) == 4
  @test CUDA.@argout(f(out(5), out(6))) == (5,6)
end


text:
module _XDistEn
export XDistEn
using StatsBase: fit, Histogram, skewness
using Statistics: mean, std
    """
        XDist, Ppi = XDistEn(Sig)

    Returns the cross-distribution entropy estimate (`XDist`) and the
    corresponding distribution probabilities (`Ppi`) estimated between the data
    sequences contained in `Sig` using the default parameters:
    embedding dimension = 2, time delay = 1, binning method = 'Sturges',
    logarithm = base 2, normalisation = w.r.t # of histogram bins

        XDist, Ppi = XDistEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1, Bins::Union{Int,String}="Sturges", Logx::Real=2, Norm::Bool=true)

    Returns the cross-distribution entropy estimate (`XDist`) estimated between the
    data sequences contained in `Sig` using the specified 'keyword' = arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer   [default: 2]    \n
    `tau`   - Time Delay, a positive integer            [default: 1]    \n
    `Bins`  - Histogram bin selection method for distance distribution,
              an integer > 1 indicating the number of bins, or one of the
              following strings {'sturges','sqrt','rice','doanes'}
              [default: 'sturges'] \n
    `Logx`  - Logarithm base, a positive scalar         [default: 2]
              ** enter 0 for natural log**\n
    `Norm`  - Normalisation of DistEn value:
              [false]  no normalisation.
              [true]   normalises w.r.t # of histogram bins [default] \n

    # See also `XSampEn`, `XApEn`, `XPermEn`, `XCondEn`, `DistEn`, `DistEn2D`, `XMSEn`

    # References:
        [1] Yuanyuan Wang and Pengjian Shang,
            "Analysis of financial stock markets through the multiscale
            cross-distribution entropy based on the Tsallis entropy."
            Nonlinear Dynamics
            94.2 (2018): 1361-1376.

    """
    function XDistEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1,
        Bins::Union{Int,String}="Sturges", Logx::Real=2, Norm::Bool=true)

    (size(Sig,2) > size(Sig,1)) ? Sig = transpose(Sig) : nothing
    Logx == 0  ? Logx = exp(1) : nothing

    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 0) ? nothing : error("m:     must be an integer > 0")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx:     must be a positive number > 0")
    if typeof(Bins)<:Int
        (Bins>1) ? nothing : error("Bins:  must be an integer > 1 (or name of binning method)")
    elseif typeof(Bins)<:String
        (lowercase(Bins) in ["sturges","sqrt","rice","doanes"]) ? nothing :
        error("Bins:    must be one of the following strings
                'sturges', 'sqrt', 'rice', 'doanes' (or an integer >1)")
    end

    S1 = Sig[:,1]; S2 = Sig[:,2];
    Nx = size(S1,1) - ((m-1)*tau)
    Zm1 = zeros(Nx,m)
    Zm2 = zeros(Nx,m)
    for n = 1:m
        Zm1[:,n] = S1[(n-1)*tau + 1:Nx+(n-1)*tau]
        Zm2[:,n] = S2[(n-1)*tau + 1:Nx+(n-1)*tau]
    end

    DistMat = zeros(Nx,Nx);
    for k = 1:Nx
        DistMat[k,:] = maximum(abs.(repeat(transpose(Zm1[k,:]),outer=Nx) - Zm2),dims=2)
    end

    Ny = Nx*Nx
    DistMat = reshape(DistMat,1,Ny)
    if eltype(Bins)<:Char
        if lowercase(Bins) == "sturges"
            Bx = ceil(log2(Ny) + 1)
        elseif lowercase(Bins) == "rice"
            Bx = ceil(2*(Ny^(1/3)));
        elseif lowercase(Bins) == "sqrt"
            Bx = ceil(sqrt(Ny));
        elseif lowercase(Bins) == "doanes"
            sigma = sqrt(6*(Ny-2)/((Ny+1)*(Ny+3)));
            Bx = ceil(1+log2(Ny)+log2(1+abs(skewness(DistMat)/sigma)));
        else
            error("Please enter a valid binning method")
        end
    else
        Bx = Bins
    end

    By = collect(range(minimum(DistMat),maximum(DistMat),length=Int(Bx+1)))
    By[end] += 1; By[1] -= 1
#    Ppi = fit(Histogram, transpose(DistMat), By).weights/Ny
    Ppi = fit(Histogram, DistMat[:], By).weights/Ny

    if round(sum(Ppi),digits=6) != 1
        @warn("Potential error estimating probabilities (p = $(sum(Ppi))")
        Ppi = Ppi[Ppi.!=0]
    elseif any(Ppi.==0)
        print("Note: $(sum(Ppi.==0))/$(length(Ppi)) bins were empty")
        Ppi = Ppi[Ppi.!=0]
    end
    XDist = -sum(Ppi.*log.(Logx, Ppi))
    Norm ? XDist = XDist/log(Logx, Bx) : nothing

    return XDist, Ppi
    end
end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
function init_grid_bot(client, config, bot_path)
    @info(now())
    @info(config["name"])
    @info("Preparation")
    @info("-------------")

    #
    # Add min_qty, precision, sell_precision, min_step, tick_size
    if !haskey(config, "precision") || !haskey(config, "sell_precision")
        @info("Adding precisions")
        buy_precision, sell_precision, min_qty, min_step, tickSize = Pinguu.get_precision(client, config["symbol"])
        config["precision"] = buy_precision
        config["sell_precision"] = sell_precision
        config["min_qty"] = min_qty
        config["min_step"] = min_step
        config["tick_size"] = tickSize
        open(joinpath(bot_path, "config.json"), "w") do f
            JSON.print(f, config, 4)
        end
    end

    #
    ## Place initial order
    if config["init_state"]["status"] == "not placed"
        @info("Placing initial order")
        order = try Pinguu.place_init_grid_order(client, config)
            catch
            nothing
        end
        if order !== nothing
            config["init_state"]["status"] = "placed"
            config["init_state"]["order"] = order
            open(joinpath(bot_path, "config.json"), "w") do f
                JSON.print(f, config, 4)
            end
        end
        return
    end

    #
    ## Check initial order. Change to "FILLED" if order has been fileed
    if config["init_state"]["status"] == "placed"
        @info("Checking initial order")
        tmp = try Pinguu.check_init_grid_order(client, config)
            catch
            nothing
        end
        @info(tmp["status"])
        if tmp !== nothing
            if tmp["status"] == client.ORDER_STATUS_FILLED
                config["init_state"]["status"] = client.ORDER_STATUS_FILLED
                config["init_state"]["order"] = [tmp]
                open(joinpath(bot_path, "config.json"), "w") do f
                    JSON.print(f, config, 4)
                end
            end
        end
        return
    end

    #
    ## Calculate average price, if several orders have been executed
    if config["init_state"]["status"] == client.ORDER_STATUS_FILLED
        @info("Evaluating filled initial order")
        vals = []
        prices = []
        executedQty = []
        for order in config["init_state"]["order"]
            cQty   = typeof(order["cummulativeQuoteQty"]) != Float64 ? parse(Float64, order["cummulativeQuoteQty"]) : order["cummulativeQuoteQty"]
            oPrice = typeof(order["price"]) != Float64 ? parse(Float64, order["price"]) : order["price"]
            exQty  = typeof(order["executedQty"]) != Float64 ? parse(Float64, order["executedQty"]) : order["executedQty"]
            push!(vals, cQty)
            push!(prices, oPrice)
            push!(executedQty, exQty)
        end
        avg_price = sum(vals) / sum(vals ./ prices)
        avg_price = Pinguu.price2string(avg_price, config["sell_precision"])
        config["init_state"]["price"] = avg_price
        config["init_state"]["cummulativeQuoteQty"] = sum(vals) * (1-0.00075) # Estimate Fees
        config["init_state"]["executedQty"] = sum(executedQty)
        config["init_state"]["status"] = "order completed"

        #
        ## Add fees of inital orders
        i = 1
        for i in eachindex(config["init_state"]["order"])
            order = config["init_state"]["order"][i]
            if !haskey(order, "trades")
                trades = client.get_my_trades(symbol=config["symbol"], orderId=order["orderId"])
                config["init_state"]["order"][i]["trades"] = trades
                open(joinpath(bot_path, "config.json"), "w") do f
                    JSON.print(f, config, 4)
                end
            end
            if !haskey(order, "fee")
                order["fee"] = calc_order_fee(client, order["trades"], config)
                open(joinpath(bot_path, "config.json"), "w") do f
                    JSON.print(f, config, 4)
                end
            end
        end
        if !haskey(config["init_state"], "fee")
            fee = 0
            for order in config["init_state"]["order"]
                fee += order["fee"]
            end
            config["init_state"]["fee"] = fee
        end
        config["init_state"]["sold"] = 0
        open(joinpath(bot_path, "config.json"), "w") do f
            JSON.print(f, config, 4)
        end
        return
    end

    #
    ## Create folders
    !isdir(joinpath(bot_path, "closed_orders")) ? mkpath(joinpath(bot_path, "closed_orders")) : "dir exists"
    !isdir(joinpath(bot_path, "closed_init_orders")) ? mkpath(joinpath(bot_path, "closed_init_orders")) : "dir exists"
    !isdir(joinpath(bot_path, "history")) ? mkpath(joinpath(bot_path, "history")) : "dir exists"
    !isdir(joinpath(bot_path, "temp")) ? mkpath(joinpath(bot_path, "temp")) : "dir exists"

    #
    ## Create Grid
    if !isfile(joinpath(bot_path, "trade_grid.json"))
        trade_grid = try Pinguu.create_trade_grid(client, config, bot_path)
        catch
            @info("Trade grid creation failed")
            return
        end
        open(joinpath(bot_path, "trade_grid.json"), "w") do f
            JSON.print(f, trade_grid, 4)
        end
        return
    end

    #
    ## Start trading! :-9
    config["status"]       = "active"
    config["start"]        = string(now())
    config["thes_factor"]  = 1.0
    config["thes_volume"]  = 0.0
    config["extra_volume"] = Dict()
    config["bot_path"]     = bot_path
    open(joinpath(bot_path, "config.json"), "w") do f
        JSON.print(f, config, 4)
    end
end
export init_grid_bot


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_subsidy_revenue(path::AbstractString, inputs::Dict, setup::Dict, dfCap::DataFrame, EP::Model)

Function for reporting subsidy revenue earned if a generator specified `Min_Cap` is provided in the input file. GenX will print this file only the shadow price can be obtained form the solver. Do not confuse this with the Minimum Capacity Carveout constraint, which is for a subset of generators, and a separate revenue term will be calculated in other files. The unit is \$.
"""
function write_subsidy_revenue(path::AbstractString, inputs::Dict, setup::Dict, dfCap::DataFrame, EP::Model)
	dfGen = inputs["dfGen"]
	#NumberOfMinCapReqs = inputs["NumberOfMinCapReqs"]

	dfSubRevenue = DataFrame(region = dfGen[!,:region], Resource = inputs["RESOURCES"], zone = dfGen[!,:Zone], Cluster = dfGen[!,:cluster], R_ID = dfGen[!,:R_ID])
	#dfSubRevenue.SubsidyRevenue .= 0.0

	dfSubRevenue.SubsidyRevenue = zeros(size(dfSubRevenue, 1))
	#dfSubRevenue[:,:SubsidyRevenue] = zeros(size(dfSubRevenue, 1))

	#dfSubRevenue[!,:SubsidyRevenue] .= 0.0
	for y in (dfGen[(dfGen[!,:Min_Cap_MW].>0) ,:][!,:R_ID])
		dfSubRevenue[y,:SubsidyRevenue] = (value.(EP[:eTotalCap])[y]) * (dual.(EP[:cMinCap])[y])
	end

	if setup["ParameterScale"] == 1
		dfSubRevenue.SubsidyRevenue = dfSubRevenue.SubsidyRevenue*(ModelScalingFactor^2) #convert from Million US$ to US$
	end
	### calculating tech specific subsidy revenue

	dfRegSubRevenue = DataFrame(region = dfGen[!,:region], Resource = inputs["RESOURCES"], zone = dfGen[!,:Zone], Cluster = dfGen[!,:cluster], R_ID = dfGen[!,:R_ID])
	dfRegSubRevenue.SubsidyRevenue = zeros(size(dfRegSubRevenue, 1))
	if (setup["MinCapReq"] >= 1)
		for mincap in 1:inputs["NumberOfMinCapReqs"] # This key only exists if MinCapReq >= 1, so we can't get it at the top outside of this condition.
			for y in dfGen[(dfGen[!,Symbol("MinCapTag_$mincap")].== 1) ,:][!,:R_ID]
			   dfRegSubRevenue[y,:SubsidyRevenue] = (value.(EP[:eTotalCap])[y]) * (dual.(EP[:cZoneMinCapReq])[mincap])
			end
		end
	end

	if setup["ParameterScale"] == 1
		dfRegSubRevenue.SubsidyRevenue = dfRegSubRevenue.SubsidyRevenue*(ModelScalingFactor^2) #convert from Million US$ to US$
	end

	CSV.write(joinpath(path, "SubsidyRevenue.csv"), dfSubRevenue)
	CSV.write(joinpath(path, "RegSubsidyRevenue.csv"), dfRegSubRevenue)
	return dfSubRevenue, dfRegSubRevenue
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
"""
    LCMType

Supertype of concrete Julia `struct`s that represent LCM message types.

Subtypes must be `mutable struct`s and may use the following field types:

* `Bool`
* numeric types: `Int8`, `Int16`, `Int32`, `Int64`, `Float32`, `Float64`
* bytes (encoded in the same was as `Int8`): `UInt8`;
* `String`;
* another `LCMType`;
* `Vector` or a subtype of `StaticVector`, for which the element type must also be
one of the previously specified types or another `Vector` or `StaticVector`.

Any size fields must come **before** the `Vector` fields to which they correspond.

The following methods must be defined for a concrete subtype of `LCMType` (say `MyType`):

* `dimensions(x::MyType, ::Val{fieldsym})`
* `fingerprint(::Type{MyType})`

Note that the `@lcmtypesetup` macro can be used to generate these methods automatically.

Also note that ideally, all of these methods would be generated from the LCM message type
definition, but that is currently not the case.
"""
abstract type LCMType end

"""
    dimensions(x::LCMType, ::Val{fieldsym})

Return a tuple of LCMDimensions describing the size of `getfield(x, fieldsym)`.
"""
function dimensions end

"""
    fingerprint(::Type{T}) where T<:LCMType

Return the fingerprint of LCM type `T` as an Int64 (exactly 8 bytes).
"""
function fingerprint end

# Types that are encoded in network byte order, as dictated by the LCM type specification.
const NetworkByteOrderEncoded = Union{Int8, Int16, Int32, Int64, Float32, Float64, UInt8}

# Julia types that correspond to LCM primitive types
const LCMPrimitive = Union{Int8, Int16, Int32, Int64, Float32, Float64, String, Bool, UInt8}

# Default values for all of the possible field types of an LCM type:
defaultval(::Type{T}) where {T<:LCMType} = T()
defaultval(::Type{Bool}) = false
defaultval(::Type{T}) where {T<:NetworkByteOrderEncoded} = zero(T)
defaultval(::Type{String}) = ""
defaultval(::Type{Array{T, N}}) where {T, N} = Array{T, N}(undef, ntuple(i -> 0, Val(N))...)
defaultval(::Type{SA}) where {T, SA<:StaticArray{<:Any, T}} = _defaultval(SA, T, Length(SA))
_defaultval(::Type{SA}, ::Type{T}, ::Length{L}) where {SA, T, L} = SA(ntuple(i -> defaultval(T), Val(L)))

# Generated default constructor for LCMType subtypes
@generated function (::Type{T})() where T<:LCMType
    constructor_arg_exprs = [:(defaultval(fieldtype(T, $i))) for i = 1 : fieldcount(T)]
    :(T($(constructor_arg_exprs...)))
end

# Field dimension information
abstract type LCMDimension end

struct LCMDimensionVar{T} <: LCMDimension end
LCMDimensionVar(fieldname::Symbol) = LCMDimensionVar{fieldname}()
sizestring(::LCMDimensionVar{T}) where {T} = string(T)

struct LCMDimensionConst{T} <: LCMDimension end
LCMDimensionConst(constsize::Integer) = LCMDimensionConst{Int(constsize)}()
sizestring(::LCMDimensionConst{T}) where {T} = string(T)

@enum DimensionMode LCM_CONST LCM_VAR
dimensionmode(::Type{<:LCMDimensionVar}) = LCM_VAR
dimensionmode(::Type{<:LCMDimensionConst}) = LCM_CONST

makedim(fieldname::Symbol) = LCMDimensionVar(fieldname)
makedim(constsize::Int) = LCMDimensionConst(constsize)

@inline evaldims(x::LCMType) = ()
@inline function evaldims(x::LCMType, dimhead::LCMDimensionConst{constsize}, dimtail::LCMDimension...) where {constsize}
    (Int(constsize), evaldims(x, dimtail...)...)
end
@inline function evaldims(x::LCMType, dimhead::LCMDimensionVar{fieldname}, dimtail::LCMDimension...) where {fieldname}
    (Int(getfield(x, fieldname)), evaldims(x, dimtail...)...)
end

# Hash computation
# Port of https://github.com/ZeroCM/zcm/blob/e9c7bfc401ea15aa64291507da37d1163e5506c0/gen/ZCMGen.cpp#L73-L94
# NOTE: actual LCM implementation: https://github.com/lcm-proj/lcm/blob/992959adfbda78a13a858a514636e7929f27ed16/lcmgen/lcmgen.c#L114-L132
# uses an int64_t for v and thus relies on undefined behavior; see https://stackoverflow.com/a/4009922.
function sign_extended_right_shift(val::UInt64, nshift::UInt64)
    (val >> 63) == 0 ? val >> nshift :
        (val >> nshift) | ~((1 << (UInt64(64) - nshift)) - UInt64(1))
end

hashupdate(v::UInt64, c::Cchar) = ((v << 8) ⊻ sign_extended_right_shift(v, UInt64(55))) + c
hashupdate(v::UInt64, c::Char) = hashupdate(v, Cchar(c))
hashupdate(v::UInt64, n::Integer) = hashupdate(v, Cchar(n))

function hashupdate(v::UInt64, s::AbstractString)
    v = hashupdate(v, Cchar(length(s)))
    for char in s
        v = hashupdate(v, char)
    end
    v
end

hashupdate(v::UInt64, sym::Symbol) = hashupdate(v, string(sym))
hashupdate(v::UInt64, mode::DimensionMode) = hashupdate(v, Cchar(mode))
function hashupdate(v::UInt64, dim::T) where {T<:LCMDimension}
    v = hashupdate(v, dimensionmode(T))
    hashupdate(v, sizestring(dim))
end

isprimitive(::Type{<:LCMType}) = false
isprimitive(::Type{T}) where {T<:LCMPrimitive} = true
isprimitive(::Type{<:AbstractArray{T}}) where {T} = isprimitive(T)

lcmtypename(::Type{Int8}) = "int8_t"
lcmtypename(::Type{Int16}) = "int16_t"
lcmtypename(::Type{Int32}) = "int32_t"
lcmtypename(::Type{Int64}) = "int64_t"
lcmtypename(::Type{Float32}) = "float"
lcmtypename(::Type{Float64}) = "double"
lcmtypename(::Type{String}) = "string"
lcmtypename(::Type{Bool}) = "boolean"
lcmtypename(::Type{UInt8}) = "byte"
lcmtypename(::Type{<:AbstractArray{T}}) where {T} = lcmtypename(T)
lcmtypename(::Type{T}) where {T<:LCMType} = string(T)

# port of https://github.com/lcm-proj/lcm/blob/992959adfbda78a13a858a514636e7929f27ed16/lcmgen/lcmgen.c#L248-L282
function basehash(::Type{T}) where T<:LCMType
    v = UInt64(0x12345678)
    for field in fieldnames(T)
        v = hashupdate(v, field)
        F = fieldtype(T, field)
        if isprimitive(F)
            v = hashupdate(v, lcmtypename(F))
        end
        dims = dimensions(T, Val(field))
        v = hashupdate(v, length(dims))
        for dim in dims
            v = hashupdate(v, dim)
        end
    end
    v
end

computehash(::Type{T}, parents::Vector{DataType}) where {T<:LCMPrimitive} = zero(UInt64)
computehash(::Type{<:AbstractArray{T}}, parents::Vector{DataType}) where {T} = computehash(T, parents)
function computehash(::Type{T}, parents::Vector{DataType}) where T<:LCMType
    T in parents && return zero(UInt64)
    hash = basehash(T)
    for field in fieldnames(T)
        F = fieldtype(T, field)
        hash += computehash(F, vcat(parents, T))
    end
    (hash << 1) + ((hash >> 63) & 1)
end

struct FingerprintException <: Exception
    T::Type
end

@noinline function Base.showerror(io::IO, e::FingerprintException)
    print(io, "LCM message fingerprint did not match type ", e.T, ". ")
    print(io, "This means that you are trying to decode the wrong message type, or a different version of the message type.")
end

function checkfingerprint(io::IO, ::Type{T}) where T<:LCMType
    fp = fingerprint(T)
    fpint = fp isa Int64 ? fp : ntoh(reinterpret(Int64, Vector(fp))[1]) # TODO: remove for next release; fingerprints are now Int64s
    decodefield(io, Int64) == fpint || throw(FingerprintException(T))
end

# Resizing
@generated function Base.resize!(x::T) where T<:LCMType
    exprs = Expr[]
    for fieldname in fieldnames(T)
        F = fieldtype(T, fieldname)
        if F <: Array
            push!(exprs, quote
                resizearrayfield!(x, $(Val(fieldname)), $F)
            end)
        end
    end
    quote
        $(exprs...)
        return nothing
    end
end

@generated function resizearrayfield!(x::T, ::Val{fieldname}, ::Type{Array{S, N}}) where {T<:LCMType, fieldname, S, N}
    quote
        Base.@_inline_meta
        newsize = evaldims(x, dimensions(T, $(Val(fieldname)))...)
        if newsize !== size(x.$fieldname)
            x.$fieldname = Array{S, N}(undef, newsize...)
            @inbounds for i in eachindex(x.$fieldname)
                x.$fieldname[i] = defaultval(S)
            end
        end
    end
end

# check_valid
"""
    check_valid(x::LCMType)

Check that the array sizes of `x` match their corresponding size fields.
"""
@generated function check_valid(x::T) where T<:LCMType
    exprs = Expr[]
    for fieldname in fieldnames(T)
        F = fieldtype(T, fieldname)
        if F <: Array
            push!(exprs, quote
                if size(x.$fieldname) != evaldims(x, LCMCore.dimensions(T, $(Val(fieldname)))...)
                    throw(DimensionMismatch())
                end
            end)
        end
    end
    quote
        $(exprs...)
        return nothing
    end
end

# Decoding
function decode!(x::LCMType, io::IO)
    checkfingerprint(io, typeof(x))
    decodefield!(x, io)
end

"""
    decode_in_place(T)

Specify whether type `T` should be decoded in place, i.e. whether to use a
`decodefield!` method instead of a `decodefield` method.
"""
function decode_in_place end

decode_in_place(::Type{<:LCMType}) = true
@generated function decodefield!(x::T, io::IO) where T<:LCMType
    field_assignments = Vector{Expr}(undef, fieldcount(x))
    for (i, fieldname) in enumerate(fieldnames(x))
        F = fieldtype(x, fieldname)
        field_assignments[i] = quote
            $F <: Array && resizearrayfield!(x, $(Val(fieldname)), $F)
            if decode_in_place($F)
                decodefield!(x.$fieldname, io)
            else
                x.$fieldname = decodefield(io, $F)
            end
        end
    end
    return quote
        $(field_assignments...)
        return x
    end
end

decode_in_place(::Type{Bool}) = false
decodefield(io::IO, ::Type{Bool}) = read(io, UInt8) == 0x01

decode_in_place(::Type{<:NetworkByteOrderEncoded}) = false
decodefield(io::IO, ::Type{T}) where {T<:NetworkByteOrderEncoded} = ntoh(read(io, T))

decode_in_place(::Type{String}) = false
function decodefield(io::IO, ::Type{String})
    len = ntoh(read(io, UInt32))
    ret = String(read(io, len - 1))
    read(io, UInt8) # strip off null
    ret
end

decode_in_place(::Type{<:Array}) = true
function decodefield!(x::Array{T}, io::IO) where T
    if decode_in_place(T)
        @inbounds for i in reversedimindices(x)
            decodefield!(x[i], io)
        end
    else
        @inbounds for i in reversedimindices(x)
            x[i] = decodefield(io, T)
        end
    end
    x
end

decode_in_place(::Type{SA}) where {T, SA<:StaticArray{<:Any, T}} = decode_in_place(T)
function decodefield!(x::StaticArray{<:Any, T}, io::IO) where T
    decode_in_place(T) || error()
    @inbounds for i in reversedimindices(x)
        decodefield!(x[i], io)
    end
    x
end
@inline decodefield(io::IO, ::Type{SA}) where {SA<:StaticArray} = _decodefield(io, Size(SA), SA)
@generated function _decodefield(io::IO, ::Size{s}, ::Type{SA}) where {s, T, SA<:StaticArray{<:Any, T}}
    decode_in_place(T) && error()
    exprs = [:(decodefield(io, T)) for i = 1 : prod(s)]
    quote
        Base.@_inline_meta
        reversedims($SA(tuple($(exprs...))))
    end
end

# Encoding
"""
    encode(io::IO, x::LCMType)

Write an LCM byte representation of `x` to `io`.
"""
function encode(io::IO, x::LCMType)
    encodefield(io, fingerprint(typeof(x)))
    encodefield(io, x)
end

@generated function encodefield(io::IO, x::LCMType)
    encode_exprs = Vector{Expr}(undef, fieldcount(x))
    for (i, fieldname) in enumerate(fieldnames(x))
        encode_exprs[i] = :(encodefield(io, x.$fieldname))
    end
    quote
        check_valid(x)
        $(encode_exprs...)
        io
    end
end

encodefield(io::IO, x::Bool) = write(io, ifelse(x, 0x01, 0x00))
encodefield(io::IO, x::NetworkByteOrderEncoded) = write(io, hton(x))
function encodefield(io::IO, x::String)
    write(io, hton(UInt32(length(x) + 1)))
    write(io, x)
    write(io, UInt8(0))
end
function encodefield(io::IO, A::AbstractArray)
    for i in reversedimindices(A)
        encodefield(io, A[i])
    end
end

# Sugar
encode(data::AbstractVector{UInt8}, x::LCMType) = encode(FastWriteBuffer(data), x)
encode(x::LCMType) = (stream = FastWriteBuffer(); encode(stream, x); flush(stream); take!(stream))

decode!(x::LCMType, data::AbstractVector{UInt8}) = decode!(x, FastReadBuffer(data))
decode(data::AbstractVector{UInt8}, ::Type{T}) where {T<:LCMType} = decode!(T(), data)

# @lcmtypesetup macro and related functions
"""
    @lcmtypesetup(lcmtype, dimensioninfos...)

Generate the following methods for a concrete LCMType subtype (say `MyType`):

* `dimensions(x::MyType, ::Val{fieldsym})`, for all fields
* `fingerprint(::Type{MyType})`

The `lcmtype` argument should be the name of a concrete LCMType subtype.
The `dimensioninfos` arguments can be used to define which fields determine
the size of which variable-size array fields. Each `dimensioninfos` element
should have the form

```julia
arrayfieldname => (size1, size2, ...)
```

where `arrayfieldname` is the name of an `Array` field and `size1`, `size2`
etc. are the dimensions of the `Array`, where both integers and field names
that come *before* `arrayfieldname` in the type definition may be used.

# Examples
```julia
mutable struct MyType <: LCMType
    alength::Int32
    c_inner_length::Int32
    a::Vector{Float64}
    b::SVector{3, Float32}
    c::Matrix{Int64}
end

@lcmtypesetup(MyType,
    a => (alength, ),
    c => (3, c_inner_length)
)
```
"""
macro lcmtypesetup(lcmt, dimensioninfos...)
    # LCMCore.dimensions methods for variable dimensions
    vardimmethods = map(dimensioninfos) do dimensioninfo
        @assert dimensioninfo.head == :call
        @assert dimensioninfo.args[1] == :(=>)
        vecfieldname = dimensioninfo.args[2]::Symbol
        @assert dimensioninfo.args[3].head == :tuple
        dims = dimensioninfo.args[3].args
        quote
            let dimtuple = tuple($(LCMCore.makedim.(dims)...))
                LCMCore.dimensions(::Type{$(esc(lcmt))}, ::Val{$(QuoteNode(vecfieldname))}) = dimtuple
            end
        end
    end

    # LCMCore.dimensions methods for constant dimensions
    makeconstdimmethods = quote
        let T = $(esc(lcmt))
            for field in fieldnames(T)
                F = fieldtype(T, field)
                if F <: Array
                    # skip
                elseif F <: StaticArray
                    let dimtuple = LCMCore.makedim.(size(F))
                        LCMCore.dimensions(::Type{T}, ::Val{field}) = dimtuple
                    end
                else
                    LCMCore.dimensions(::Type{T}, ::Val{field}) = ()
                end
            end
        end
    end

    # LCMCore.fingerprint method
    fingerprint = quote
        let hash = reinterpret(Int64, LCMCore.computehash($(esc(lcmt)), DataType[]))
            LCMCore.fingerprint(::Type{$(esc(lcmt))}) = hash
        end
    end

    quote
        $(vardimmethods...)
        $makeconstdimmethods
        $fingerprint
    end
end


text:
# generates a lcov.info file in the format generated by `geninfo`. This format
# can be parsed by a variety of useful utilities to display coverage info

export LCOV
"""
Coverage.LCOV Module

This module provides functionality to generate LCOV info format files from
Julia coverage data. It exports the `writefile` function.
"""
module LCOV

using Compat
using Coverage

export writefile

"""
    writefile(outfile::AbstractString, fcs)

Write the given coverage data to a file in LCOV info format. The data can either
be a `FileCoverage` instance or a vector of `FileCoverage` instances.
"""
function writefile(outfile::AbstractString, fcs)
    open(outfile, "w") do f
        write(f, fcs)
    end
end

"""
    write(io::IO, fcs)

Write the given coverage data to an `IO` stream in LCOV info format. The data
can either be a `FileCoverage` instance or a vector of `FileCoverage` instances.
"""
function write end

function write(io::IO, fcs::Vector{FileCoverage})
    for fc in fcs
        write(io, fc)
    end
end

function write(io::IO, fc::FileCoverage)
    instrumented = 0
    covered = 0
    println(io, "SF:$(fc.filename)")
    for (line, cov) in enumerate(fc.coverage)
        (lineinst, linecov) = writeline(io, line, cov)
        instrumented += lineinst
        covered += linecov > 0 ? 1 : 0
    end
    println(io, "LH:$covered")
    println(io, "LF:$instrumented")
    println(io, "end_of_record")
end

# document the writeline function instead of individual methods
"""
    writeline(io::IO, line::Int, count)

Write LCOV data for a single line to the given `IO` stream. Returns a 2-tuple
of the number of lines instrumented (0 or 1) and the count for how many times
the line was executed during testing.
"""
function writeline end

function writeline(io::IO, line::Int, count::Int)
    println(io, "DA:$line,$count")
    (1, count)
end
function writeline(io::IO, line::Int, count::Nothing)
    # skipped line, nothing to do here
    (0, 0)
end

end


text:
using Dates
using BenchmarkTools

include("../Fastback.jl");

# w/ delayed orders min     = 14.4 ms
# direct placement          = 5.4 ms
@benchmark begin
    AAPL = Instrument("AAPL");
    dt = DateTime(2018, 1, 2, 9, 30, 0)
    acc = Account(100_000.0)
    for i = 1:100_000
        nbbo = NBBO(dt, 100.0 + i / 1000, 102 + i / 1000)
        execute_order!(acc, OpenOrder(AAPL, 100), nbbo)
        update_account!(acc, AAPL, NBBO(dt, 100, 102))
        execute_order!(acc, CloseOrder(acc.open_positions[1]), nbbo)
        update_account!(acc, AAPL, NBBO(dt, 100, 102))
    end
end evals=3 samples=40

# w/ delayed orders min     = 2.1 ms
# direct placement          = 1.1 ms
@benchmark begin
    AAPL = Instrument("AAPL");
    dt = DateTime(2018, 1, 2, 9, 30, 0)
    acc = Account(100_000.0)
    for i = 1:100_000
        nbbo = NBBO(dt, 100.0 + i / 1000, 102 + i / 1000)
        if i == 1
            execute_order!(acc, OpenOrder(AAPL, 100), nbbo)
        end
        update_account!(acc, AAPL, nbbo)
    end
end evals=3 samples=40


text:
canvas = BrailleCanvas(40, 10, origin_x = 0., origin_y = 0., width = 1., height = 1.)
p = @inferred Plot(canvas)
@test_throws ArgumentError Plot(canvas, margin=-1)
test_ref(
    "references/plot/empty.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)
canvas = BrailleCanvas(30, 5, origin_x = 0., origin_y = 0., width = 1., height = 1.)
p = @inferred Plot(canvas)
test_ref(
    "references/plot/empty_small.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)

canvas = BrailleCanvas(40, 10, origin_x = 0., origin_y = 0., width = 1., height = 1.)
lines!(canvas, 0., 0., 1., 1., :blue)
lines!(canvas, .2, .7, 1., 0., :red)
lines!(canvas, 0., 2., .5, 0., :green)
points!(canvas, .5, .9)

p = @inferred Plot(canvas)
test_ref(
    "references/plot/canvas_only.txt",
    @io2str(print(IOContext(::IO, :color=>true), p))
)
test_ref(
    "references/plot/canvas_only.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)

for border in (:solid, :corners, :barplot, :bold, :ascii, :none, :dashed, :dotted)
    local p = @inferred Plot(canvas, border = border)
    test_ref(
        "references/plot/border_$(string(border)).txt",
        @io2str(show(IOContext(::IO, :color=>true), p))
    )
end

p = @inferred Plot(canvas, xlabel = "x", ylabel= "y")
test_ref(
    "references/plot/xylabel.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)

p = @inferred Plot(canvas, title = "testtitle", xlabel = "x", ylabel= "y", margin = 4, padding = 5)
@test @inferred(annotate!(p, :l, ":l auto 1")) === p
@test @inferred(annotate!(p, :r, ":r auto 1", :red)) === p
test_ref(
    "references/plot/padding.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)

p = @inferred Plot(canvas, title = "testtitle", xlabel = "x", ylabel= "y", margin = 4, padding = 5, labels = false)
@test @inferred(annotate!(p, :l, ":l auto 1")) === p
@test @inferred(annotate!(p, :r, ":r auto 1", :red)) === p
test_ref(
    "references/plot/padding_nolabels.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)

p = @inferred Plot(canvas, title = "testtitle")
test_ref(
    "references/plot/title.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)
@test @inferred(annotate!(p, :l, ":l auto 1")) === p
@test @inferred(annotate!(p, :r, ":r auto 1", :red)) === p
@test @inferred(annotate!(p, :l, ":l auto 2", :green)) === p
@test @inferred(annotate!(p, :r, ":r auto 2", :blue)) === p
test_ref(
    "references/plot/title_auto.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)

for i in 1:10
    @test @inferred(annotate!(p, :l, i, "$i", color=:green)) === p
    @test @inferred(annotate!(p, :r, i, "$i", :yellow)) === p
end
@test_throws ArgumentError annotate!(p, :bl, 5, ":l  5", :red)
@test_throws ArgumentError annotate!(p, :d, ":l auto", :red)
@test @inferred(annotate!(p, :l, 5, ":l  5", :red)) === p
@test @inferred(annotate!(p, :r, 5, "5  :r", color=:blue)) === p
@test @inferred(annotate!(p, :t, ":t", color=:yellow)) === p
@test @inferred(annotate!(p, :tl, ":tl", :cyan)) === p
@test @inferred(annotate!(p, :tr, ":tr")) === p
annotate!(p, :bl, ":bl", :blue)
annotate!(p, :b, ":b", :green)
annotate!(p, :br, ":br", :white)

@test @inferred(lines!(p, 0., 1., 1., 0., color=(0,0,255))) === p
@test @inferred(pixel!(p, 10, 1, color = (255,255,0))) === p
@test @inferred(points!(p, 0.05, .75, color = (0,255,0))) === p

ttl = "title!(plot, text)"
@test @inferred(title!(p, ttl)) === p
@test @inferred(title(p)) == ttl
xlab = "xlabel!(plot, text)"
@test @inferred(xlabel!(p, xlab)) === p
@test @inferred(xlabel(p)) == xlab
ylab = "ylabel!(plot, text)"
@test @inferred(ylabel!(p, ylab)) === p
@test @inferred(ylabel(p)) == ylab

test_ref(
    "references/plot/full_deco.txt",
    @io2str(show(IOContext(::IO, :color=>true), p))
)
test_ref(
    "references/plot/full_deco_nocolor.txt",
    @io2str(show(::IO, p))
)


text:
using ThArrays
using Test

@testset "Scalar Creation" begin

    @testset "Scalar creation" begin
        s = Scalar(1.0)
        @test s[] == 1.0
    end

    @testset "Scalar data type" begin
        for typ in [Float16, Float32, Float64,
                    Bool,
                    Int8, Int16, Int32, Int64]
            s = Scalar(typ(1))
            @test s.type == typ
        end
    end
end


text:
using PIMDtools
using Test

function test1()
    itrj_start = 4
    itrj_end = 100
    atominfo = [("Pd",20),("Ru",16),("Al",92)]
    a = 12.6238002777
    unitcell = Array{Array{Float64,1},1}(undef,3)
    for i=1:3
        unitcell[i] = zeros(Float64,3)
        unitcell[i][i] = a/a0
    end
    trj = Trajectory("trj.out",itrj_start,itrj_end,atominfo,unitcell=unitcell)
    positions = get_relative_positions(trj)

    dirname = "./"
    period = 10
    start_itrj = itrj_start
    maxsteps = 10
    calc_MSD(trj,dirname,positions,period,start_itrj,maxsteps)
end

function testbinary()
    itrj_start = 4
    itrj_end = 100
    atominfo = [("Pd",20),("Ru",16),("Al",92)]
    a = 12.6238002777
    mkpath("testdir")
    filename = "trj.out"
    dirname = "./testdir"
    num_of_structures = 5

    inputdirname = "./testdir"
    period = 10
    start_itrj = itrj_start
    maxsteps = 10

    unitcell = Array{Array{Float64,1},1}(undef,3)
    for i=1:3
        unitcell[i] = zeros(Float64,3)
        unitcell[i][i] = a/a0
    end

    toBinary_r(filename,dirname,itrj_start,itrj_end,num_of_structures,atominfo,unitcell)

    calc_MSD_fromfile(inputdirname,"./",itrj_start,
        period,start_itrj,maxsteps)
end

@testset "PIMDtools.jl" begin

    test1()
    testbinary()
    # Write your tests here.
end


text:
function init_constraint!(
    com::CS.CoM,
    constraint::XorConstraint,
    fct,
    set::XorSet;
)
    !init_lhs_and_rhs!(com, constraint, fct, set) && return false

    if constraint.complement_lhs !== nothing
        init_constraint!(com, constraint.complement_lhs, constraint.complement_lhs.fct, constraint.complement_lhs.set)
    end
    if constraint.complement_rhs !== nothing
        init_constraint!(com, constraint.complement_rhs, constraint.complement_rhs.fct, constraint.complement_rhs.set)
    end
    return true
end

"""
    function is_constraint_violated(
        com::CoM,
        constraint::BoolConstraint,
        fct,
        set::XorSet,
    )

Check if both of the inner constraints are violated or whether both are solved
"""
function is_constraint_violated(
    com::CoM,
    constraint::BoolConstraint,
    fct,
    set::XorSet,
)
    both_violated = is_lhs_constraint_violated(com, constraint) && is_rhs_constraint_violated(com, constraint)
    both_violated && return true
    lhs_solved = is_constraint_solved(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_solved = is_constraint_solved(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    return lhs_solved && rhs_solved
end

"""
    still_feasible(com::CoM, constraint::XorConstraint, fct, set::XorSet, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `vidx` to `value`.
**Attention:** This assumes that it isn't violated before.
"""
function still_feasible(
    com::CoM,
    constraint::XorConstraint,
    fct,
    set::XorSet,
    vidx::Int,
    value::Int,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    lhs_feasible = !lhs_violated
    lhs_solved = false
    if lhs_feasible
        lhs_solved = is_constraint_solved_when_fixed(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value)
    end
    if lhs_feasible && !lhs_solved
        lhs_indices = constraint.lhs.indices
        for i in 1:length(lhs_indices)
            if lhs_indices[i] == vidx
                lhs_feasible = still_feasible(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value)
                break
            end
        end
    end
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    rhs_feasible = !rhs_violated
    rhs_solved = false
    if rhs_feasible
        rhs_solved = is_constraint_solved_when_fixed(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value)
    end
    if rhs_feasible && !rhs_solved
        rhs_indices = constraint.rhs.indices
        for i in 1:length(rhs_indices)
            if rhs_indices[i] == vidx
                rhs_feasible = still_feasible(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value)
                break
            end
        end
    end

    # at least one must be feasible
    if !lhs_feasible
        return rhs_feasible
    elseif !rhs_feasible
        return lhs_feasible
    end
    # not allowed that both are already solved
    return !(lhs_solved && rhs_solved)
end

"""
    prune_constraint!(com::CS.CoM, constraint::XorConstraint, fct, set::XorSet; logs = true)

Reduce the number of possibilities given the `XorConstraint` by pruning both parts
Return whether still feasible
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::XorConstraint,
    fct,
    set::XorSet;
    logs = true,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    if lhs_violated && rhs_violated
        return false
    end
    # check if one is already solved
    lhs_solved = is_constraint_solved(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_solved = is_constraint_solved(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    if lhs_solved && rhs_solved
        return false
    end

    # if one is solved => complement prune the other
    # Todo implement for activated complement constraints
    if lhs_solved && constraint.complement_rhs !== nothing
        return prune_constraint!(com, constraint.complement_rhs, constraint.complement_rhs.fct, constraint.complement_rhs.set; logs=logs)
    end
    if rhs_solved && constraint.complement_lhs !== nothing
        return prune_constraint!(com, constraint.complement_lhs, constraint.complement_lhs.fct, constraint.complement_lhs.set; logs=logs)
    end

    # if both aren't solved yet we only prune if one is violated
    if !lhs_solved && !rhs_solved
        if lhs_violated
            activate_rhs!(com, constraint)
            return prune_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set; logs=logs)
        end
        if rhs_violated
            activate_lhs!(com, constraint)
            return prune_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set; logs=logs)
        end
    end

    return true
end


text:
import Base: (+), (-), (*), (/), (//), (^), divrem

export monic

# Addition

function (+)(p::SimplePolynomial, q::SimplePolynomial)
    n = max(degree(p), degree(q))
    T = typeof(p(1) + q(1))
    c = zeros(T, n + 1)
    for j = 0:n
        c[j+1] = p[j] + q[j]
    end
    return SimplePolynomial(c)
end

# polynomial + number
function (+)(p::SimplePolynomial, a::T) where {T<:CoefX}
    return p + SimplePolynomial(a)
end

function (+)(a::T, p::SimplePolynomial) where {T<:CoefX}
    return p + SimplePolynomial(a)
end

# Subtraction

function (-)(p::SimplePolynomial)
    coefs = [-x for x in p.data]
    return SimplePolynomial(coefs)
end

(-)(p::SimplePolynomial, q::SimplePolynomial) = p + (-q)
(-)(p::SimplePolynomial, a::T) where {T<:CoefX} = p + (-a)
(-)(a::T, p::SimplePolynomial) where {T<:CoefX} = a + (-p)

# Multiplication

function (*)(a::T, p::SimplePolynomial) where {T<:CoefX}
    coefs = [a * x for x in p.data]
    return SimplePolynomial(coefs)
end

(*)(p::SimplePolynomial, a::T) where {T<:CoefX} = a * p

function (*)(p::SimplePolynomial, q::SimplePolynomial)
    n = degree(p) + degree(q)

    if p == 0 || q == 0
        return SimplePolynomial(0)
    end

    coefs = [p[0] * q[0] for i = 0:n]

    for k = 1:n
        coefs[k+1] = 0 * coefs[k+1]  # clear it out
        for j = 0:k
            coefs[k+1] += p[j] * q[k-j]
        end
    end
    return SimplePolynomial(coefs)
end

# Integer exponentiation

function (^)(p::SimplePolynomial, k::S) where {S<:Integer}

    if k < 0
        f = inv(p)
        return f^(-k)
    end

    if k == 0
        return SimplePolynomial(Int(1))
    end

    if k == 1
        return p
    end

    if k == 2
        return p * p
    end

    if k % 2 == 0  # even exponent
        j = k ÷ 2
        q = p^j
        return q * q
    end

    j = k ÷ 2
    q = p^j
    return q * q * p
end

# Division

function (/)(p::SimplePolynomial, a::T) where {T<:CoefX}
    coefs = p.data .// a
    return SimplePolynomial(coefs)
end
(//)(p::SimplePolynomial, a::T) where {T<:CoefX} = p / a

"""
`monic(p::SimplePolynomial)` returns a new polynomial formed by
dividing all coefficients by the leading coefficient, thereby
yielding a monic polynomial. If `p` is zero, a zero polynomial
is returned.
"""
function monic(p::SimplePolynomial)::SimplePolynomial
    if p == 0
        return p
    end
    return p / lead(p)
end

function divrem(a::SimplePolynomial, b::SimplePolynomial)
    if b == 0
        DivideError()
    end

    #ST = typeof((a[0]+b[0])//1)
    ST = Int
    if eltype(a) <: AbstractMod || eltype(b) <: AbstractMod
        ST = eltype(a)
    else
        ST = typeof((lead(a) + lead(b)) // 1)
    end

    da = degree(a)
    db = degree(b)

    if da < db
        return SimplePolynomial(), a
    end

    leader = a.data[end] // b.data[end]

    coefs = zeros(ST, da - db + 1)
    coefs[end] = leader

    q = SimplePolynomial(coefs)
    aa = a - b * q
    qq, r = divrem(aa, b)

    return q + qq, r
end


text:
module thismorning

export interest

export rate

function interest(amount, rate)
 return amount * (1 + rate)
end

function rate(amount, interest)
 return interest / amount
end

end

max_stars_repo_name: jlvvlj/thismorning

max_issues_repo_name: jlvvlj/thismorning

max_forks_repo_name: jlvvlj/thismorning


text:
using Kahuna
using Test

@testset "kahuna_read" begin

    @testset ".dm3 files" begin
        @test 2 + 2 == 4
    end

    @testset ".dm4 files" begin
        @test 2 + 2 == 4
    end

    @testset ".hdf5/.h5 files" begin
        @test 2 + 2 == 4
    end

    @testset ".mat files" begin
        # matfile = "test/sample_files/test_fileio_mat.mat";
        matfile = "sample_files/test_fileio_mat.mat";
        @test typeof(kahuna_read(matfile, "mat0d")) == Float64
        @test typeof(kahuna_read(matfile, "mat1d")) == Array{Float64,2} && size(kahuna_read(matfile, "mat1d")) == (1,10)
        @test typeof(kahuna_read(matfile, "mat2d")) == Array{Float64,2} && size(kahuna_read(matfile, "mat2d")) == (10,10)
        @test typeof(kahuna_read(matfile, "mat3d")) == Array{Float64,3} && size(kahuna_read(matfile, "mat3d")) == (10,10,10)
        @test typeof(kahuna_read(matfile, "mat4d")) == Array{Float64,4} && size(kahuna_read(matfile, "mat4d")) == (10,10,10,10)

        @test kahuna_read(matfile; mode="list") == Set(["mat0d", "mat1d", "mat2d", "mat4d", "mat3d"])

        @test kahuna_read(matfile) == Dict(map(x -> x => kahuna_read(matfile, x), collect(kahuna_read(matfile; mode="list"))))
    end

    @testset ".mib files" begin

        mibfile512_12bit = "sample_files/test_512_12bit_single.mib";

        # mibfiles = [mibfile256_1bit, mibfile256_6bit, mibfile256_12bit,
        #             mibfile256_1bit_raw, mibfile256_6bit_raw, mibfile256_12bit_raw,
        #             mibfile512_1bit, mibfile512_6bit, mibfile512_12bit,
        #             mibfile512_1bit_raw, mibfile512_6bit_raw, mibfile512_12bit_raw];
        mibfiles = [mibfile512_12bit]

        for mibfile in mibfiles
            mib_images, mib_headers = kahuna_read(mibfile)
            @test typeof(mib_images) == Array{Array{UInt16,2},1}
            @test typeof(mib_headers) == Array{MIBHeader,1}
            # @test typeof(kahuna_read(mibfile, [1, 10])) == Array{Float64,2} && size(kahuna_read(mibfile, [1, 10])) == (1,10)
            # @test typeof(kahuna_read(mibfile, [1, 10])) == Array{Float64,2} && size(kahuna_read(mibfile, [1, 10])) == (10,10)
            # @test typeof(kahuna_read(mibfile, [1, 10])) == Array{Float64,2} && size(kahuna_read(mibfile, [1, 10])) == (10,10)
        end

    end

    @testset ".toml files" begin
        @test 2 + 2 == 4
    end

    @testset ".jld files" begin
        @test 2 + 2 == 4
    end

end

@testset "kahuna_write" begin

    @testset ".hdf5/.h5 files" begin
        @test 2 + 2 == 4
    end

    @testset ".toml files" begin
        @test 2 + 2 == 4
    end

    @testset ".jld files" begin
        @test 2 + 2 == 4
    end

end