text:
function process_bouts(df0::AbstractDataFrame; observe = false)
    df1 = filter!(r -> !ismissing(r.Poke) && !r.Incorrect, df0)
    bdf = combine(groupby(df1,:Patch)) do dd
        count_bout_by_patch(dd)
    end
    sort!(df1,[:Patch])
    df1[!,:Bout] = bdf.bout

    # transform!(groupby(df1,:Patch), :RewardAvailable => count_bout_bypatch => :Bout)
    df1 = df1[:,[:Richness, :Travel,:State,:Poke,:ActivePort,
        :Port, :PokeIn, :PokeOut, :Duration,
        :Patch, :Bout,:RewardAvailable, :RewardConsumption,
        :SubjectID, :Taskname, :Experimentname, :Startdate]]
    observe && open_html_table(df1)
    bouts = combine(groupby(df1,[:Bout, :Patch])) do dd
        if dd[1,:State] == "Forage"
            forage_idx = findall(ismatch.(r"^Poke",dd.Port))
        elseif dd[1,:State] == "Travel"
            forage_idx = findall(ismatch.(r"^Trav",dd.Port))
        end
        if isempty(forage_idx)
            println("Bout #$(dd[1,:Bout]), Patch #$(dd[1,:Patch]) has no pokes matching state $(dd[1,:State])")
            df2 = DataFrame(
                In = missing,
                Out = missing,
                ForageTime_total = missing,
                ForageTime_Sum = missing,
                Pokes = missing,
                Rewarded = missing,
                )
                for x in [:Patch,:State,:Richness,:Travel, :ActivePort,:SubjectID,
                        :Taskname, :Experimentname, :Startdate]
                    df2[!,x] .= dd[1,x]
                end
        else
            df2 = DataFrame(
                In = dd[forage_idx[1],:PokeIn],
                Out = dd[forage_idx[end],:PokeOut],
                ForageTime_total = dd[forage_idx[end],:PokeOut] - dd[forage_idx[1],:PokeIn],
                ForageTime_Sum = sum(dd[forage_idx,:Duration]),
                Pokes = length(forage_idx),
                Rewarded = any(.!ismissing.(dd.RewardAvailable)),
                # Giveup = any(ismatch.(r"^TravPoke",dd.Port))
            )
            for x in [:Patch,:State,:Richness,:Travel, :ActivePort,:SubjectID,
                    :Taskname, :Experimentname, :Startdate]
                df2[!,x] .= dd[forage_idx,x][1]
            end
        end
        return df2
    end
    sort!(bouts,[:Patch,:Bout,:In])
    bouts[!,:GiveUp] = vcat(ismatch.(r"^Travel",bouts[2:end,:State]),[false])
    bouts = bouts[:,[:In,:Out, :ForageTime_total, :ForageTime_Sum,
            :Pokes, :Rewarded, :GiveUp, :Bout,
            :Patch,:State,:Richness,:Travel, :ActivePort,
            :SubjectID,:Taskname, :Experimentname, :Startdate
        ]]

    observe && open_html_table(bouts)
    return bouts
end
##
function count_bout_by_patch(df1)
    bdf = df1[:,[:Port,:RewardAvailable,:State]]
    bdf[!,:bout] .= 1
    c = 1
    bdf[!,:rewarding] .= false
    rewarding = false
    bdf[!,:rewarded] .= false
    rewarded = false
    bdf[!,:traveling] .= false
    traveling = false
    for r in eachrow(bdf)
        !ismissing(r.RewardAvailable) && (rewarding = true)
        rewarding && ismatch(r"^Rew",r.Port) && (rewarded = true)
        if rewarded && ismatch(r"^Poke",r.Port)
            c += 1
            rewarded = false
            ismissing(r.RewardAvailable) ? rewarding = false : rewarding = true
        end
        if !traveling && ismatch(r"^Travel$", r.State)
            c += 1
            traveling = true
        end
        if traveling && ismatch(r"^Forage$", r.State)
            traveling = false
        end
        r.rewarding = rewarding
        r.rewarded = rewarded
        r.traveling = traveling
        r.bout = c
    end
    return bdf
end
# function count_bout(df::AbstractDataFrame)
#     res = []
#     c = 1
#     # travel = false
#     rewarded = false
#     consumed = false
#     state = df[1,:State]
#     ## add consumed condition to change bout
#     for r in eachrow(df)
#         !ismissing(r.RewardAvailable) && (rewarded = true)
#         (ismatch(r"^Rew",r.Port) && (rewarded == true)) && (consumed = true)
#         if consumed && ismatch(r"^Poke",r.Port) && r.State == "Forage"
#             rewarded = false
#             consumed = false
#             c += 1
#             # elseif travel && ismatch(r"^P",r.Port) && r.State == "Forage"
#             #     travel = false
#             #     c += 1
#         elseif r.State != state && ismatch(r"(^Poke)|(^Trav)",r.Port)
#             state = r.State
#             c += 1
#         end
#         # r.State == "Travel" ? (travel = true) : (travel = false)
#         push!(res,c)
#     end
#     return res
# end

# function count_bout_bypatch(RewardAvailable)#df::AbstractDataFrame)
#     c = 1
#     res = []
#     for r in RewardAvailable
#         push!(res,c)
#         !ismissing(r) && (c += 1)
#     end
#     return res
# end

function count_bout_bypatch(df)
    res = []
    c = 1
    rewarding = false
    rewarded = false
    traveling = false
    for r in eachrow(df)
        !ismissing(r.RewardAvailable) && (rewarding = true)
        rewarding && ismatch(r"^Rew",r.Port) && (rewarded = true)
        if rewarded && ismatch(r"^Poke",r.Port)
            c += 1
            rewarded = false
            rewarding = false
        end
        if !traveling && ismatch(r"^Travel$", r.State)
            c += 1
            traveling = true
        end
        if traveling && ismatch(r"^Forage$", r.State)
            traveling = false
        end
        push!(res,c)
    end
    df[!,:Bout] = res
    return df[:,[:Patch, :Bout,:Poke, :State, :ActivePort, :Port, :Incorrect, :PokeIn, :PokeOut, :Duration,
        :Richness, :Travel, :P, :TravelComplete, :RewardAvailable, :RewardDelivery, :RewardConsumption,
        :Rsync_count, :Rsync_time, :StateIn, :StateOut, :T, :AFT, :IFT, :SubjectID, :Taskname,
        :Experimentname, :Startdate]]
end

# function count_bout_df(df)
#     bdf = DataFrame(C = repmat([1], nrow(df)))
#     bdf[!,]
# end


text:
# Autogenerated wrapper script for Arpack_jll for aarch64-apple-darwin
export libarpack

using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Arpack")
JLLWrappers.@declare_library_product(libarpack, "@rpath/libarpack.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libarpack,
        "lib/libarpack.2.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Shoco"
version = v"0.1.0" # NOTE: Shoco itself has no tagged versions

version_sha = "4dee0fc850cdec2bdb911093fe0a6a56e3623b71"

# Collection of sources required to build Shoco
sources = [
    "https://github.com/Ed-von-Schleck/shoco/archive/$version_sha.tar.gz" =>
        "20caed48304d4dcafd04dff6db2f4adaba21943811400da5b8008f14edab7667",
]

# Bash recipe for building across all platforms
script = raw"""
    cd ${WORKSPACE}/srcdir/shoco-*
    if [[ ${target} == *apple* ]]; then
        flag=-dynamiclib
    else
        flag=-shared
    fi
    ${CC} shoco.c -o libshoco.${dlext} -fPIC -std=c99 ${flag}
    if [ ! -d ${prefix}/lib ]; then
        mkdir -p ${prefix}/lib
    fi
    cp libshoco.* ${prefix}/lib
    """

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    FreeBSD(:x86_64),
    Linux(:aarch64, :glibc),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:armv7l, :musl, :eabihf),
    Linux(:i686, :glibc),
    Linux(:i686, :musl),
    Linux(:powerpc64le, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:x86_64, :musl),
    MacOS(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [LibraryProduct(prefix, "libshoco", :libshoco)]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# Performance analytics

## Returns
function returns_annualized(rets::Returns)::Float64
    return returns_annualized(values(rets))
end

function returns_annualized(rets::AbstractArray{Float64})::Float64
    n = 252 / length(rets)
    overall_return = last(cumprod(1 .+ rets))
    return round((overall_return ^ (n))-1, digits = 3)
end

## Volatility

function std_annualized(rets::TimeArray)::Float64
    return std_annualized(values(rets))
end

function std_annualized(rets::AbstractArray{Float64})::Float64
    return sqrt(252) * std(rets)
end

function semid_annualized(rets::TimeArray)::Float64
    return semid_annualized(values(rets))
end

function semid_annualized(rets::AbstractArray{Float64})::Float64
    return sqrt(252) * semideviation(rets)
end

## Frequently used quantiative measurements

function sharpe_annualized(rets::Union{Returns, AbstractArray{Float64}})::Float64
    return round(returns_annualized(rets) / std_annualized(rets), digits = 3)
end

function sortino_annualized(rets::Union{Returns, AbstractArray{Float64}})::Float64
    return round(returns_annualized(rets) / semid_annualized(rets), digits = 3)
end

## Drawdowns, etc.

# function drawdown(ta::TimeArray)
#     ta_cumprod = returns_cumprod(ta)
#     return upto(maximum, ta_cumprod)
# end

## Rolling metrics

function rolling_returns(t::TimeArray)::TimeArray
    return moving(return_annualized, t, 252)
end

function rolling_std(t::TimeArray)::TimeArray
    return moving(std_annualized, t, 252)
end

function rolling_sharpe(t::TimeArray)::TimeArray
    return moving(sharpe_annualized, t, 252)
end


text:
import Pkg
import UUIDs

@inline is_jll_name(name::AbstractString) = endswith(name, "_jll")

"""
    compress_versions(pool::Vector{VersionNumber}, subset::Vector{VersionNumber})

Given `pool` as the pool of available versions (of some package) and `subset` as some
subset of the pool of available versions, this function computes a `VersionSpec` which
includes all versions in `subset` and none of the versions in its complement.
"""
@inline function compress_versions(pool::Vector{VersionNumber},
                                   subset::Vector{VersionNumber})
    # Explicitly drop prerelease/build numbers, as those can confuse this.
    # TODO: Rewrite all this to use VersionNumbers instead of VersionBounds
    drop_build_prerelease(v::VersionNumber) = VersionNumber(v.major, v.minor, v.patch)
    pool = drop_build_prerelease.(pool)
    subset = sort!(drop_build_prerelease.(subset))

    complement = sort!(setdiff(pool, subset))
    ranges = Pkg.Types.VersionRange[]
    @label again
    isempty(subset) && return Pkg.Types.VersionSpec(ranges)
    a = first(subset)
    for b in reverse(subset)
        a.major == b.major || continue
        for m = 1:3
            lo = Pkg.Types.VersionBound((a.major, a.minor, a.patch)[1:m]...)
            for n = 1:3
                hi = Pkg.Types.VersionBound((b.major, b.minor, b.patch)[1:n]...)
                r = Pkg.Types.VersionRange(lo, hi)
                if !any(v in r for v in complement)
                    filter!(!in(r), subset)
                    push!(ranges, r)
                    @goto again
                end
            end
        end
    end
end

# @inline function compress_versions(pool::Vector{VersionNumber}, subset)
#     compress_versions(pool, filter(in(subset), pool))
# end

@inline function load_versions(path::String)
    versions_file = joinpath(dirname(path), "Versions.toml")
    versions_dict = Pkg.TOML.parsefile(versions_file)
    sort!([VersionNumber(v) for v in keys(versions_dict)])
end

@inline function load(path::String,
                      versions::Vector{VersionNumber} = load_versions(path))
    compressed = Pkg.TOML.parsefile(path)
    uncompressed = Dict{VersionNumber,Dict{Any,Any}}()
    for (vers, data) in compressed
        vs = Pkg.Types.VersionSpec(vers)
        for v in versions
            v in vs || continue
            merge!(get!(uncompressed, v, Dict()), deepcopy(data))
        end
    end
    return uncompressed
end

@inline function compress(path::String,
                          uncompressed::Dict,
                          versions::Vector{VersionNumber} = load_versions(path))
    inverted = Dict()
    for (ver, data) in uncompressed, (key, val) in data
        val isa Pkg.TOML.TYPE || (val = string(val))
        push!(get!(inverted, key => val, VersionNumber[]), ver)
    end
    compressed = Dict()
    for ((k, v), vers) in inverted
        for r in compress_versions(versions, sort!(vers)).ranges
            get!(compressed, string(r), Dict{String,Any}())[k] = v
        end
    end
    return compressed
end

@inline function save(path::String,
                      uncompressed::Dict,
                      versions::Vector{VersionNumber} = load_versions(path))
    compressed = compress(path, uncompressed)
    open(path, write=true) do io
        Pkg.TOML.print(io, compressed, sorted=true)
    end
end

@inline function _save_uncompressed(path::String,
                                    uncompressed::Dict,
                                    versions::Vector{VersionNumber} = load_versions(path))
    uncompressed_str = _make_keys_strings(uncompressed)
    open(path, write=true) do io
        Pkg.TOML.print(io, uncompressed_str, sorted=true)
    end
end

@inline function _make_keys_strings(uncompressed::Dict)
    result = Dict{Any, Any}()
    for k in keys(uncompressed)
        v = uncompressed[k]
        new_k = _my_string(k)
        result[new_k] = v
    end
    return result
end

@inline function _my_string(x)::String
    return String(x)
end

@inline function _my_string(v::VersionNumber)::String
    return string(v)
end


text:
using SnoopCompile, Pkg
project_path = dirname(@__FILE__)

### Log the compiles
# This only needs to be run once (to generate "/tmp/comonicon_compiles.log")

SnoopCompile.@snoopc ["--project=$(project_path)"] "/tmp/ion_compiles.log" begin
    using IonCLI, Comonicon
    include(Comonicon.PATH.project(IonCLI, "deps", "precompile.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts

data = SnoopCompile.read("/tmp/ion_compiles.log")
pc = SnoopCompile.parcel(reverse!(data[2]))

open(joinpath(project_path, "deps", "statements.jl"), "w+") do io
    for (k, v) in pc
        for ln in sort(v)
            println(io, ln)
        end
    end
end


text:
############################
# Specialized factorizations
############################

export idfact

@doc doc"""
An interpolative decomposition.

For a matrix `A`, the interpolative decomposition `F` contains the matrices `B`
and `P` computed by `idfact()`. See the documentation of `idfact()` for details.

References:

    \cite{Cheng2005, Liberty2007}
""" ->
immutable Interpolative{T} <: Factorization{T}
    B :: AbstractMatrix{T}
    P :: AbstractMatrix{T}
end

@doc doc"""
Compute the interpolative decomposition of `A`

    $$
    A ≈ B * P
    $$

    where:
    - `B`'s columns are a subset of the columns of `A`
    - some subset of `P`'s columns are the `k x k` identity,
    no entry of `P` exceeds magnitude 2, and
    - $$||B * P - A|| ≲ σ(A, k+1)$$, the (`k+1`)st singular value of `A`.

Inputs:

    `A`: Matrix to factorize
    `k`: Number of columns of A to return in B
    `l`: Length of random vectors to project onto

Implementation note:

    This is a hacky version of the algorithms described in \cite{Liberty2007}
    and \cite{Cheng2005}. The former refers to the factorization (3.1) of the
    latter.  However, it is not actually necessary to compute this
    factorization in its entirely to compute an interpolative decomposition.
    Instead, it suffices to find some permutation of the first k columns of Y =
    R * A, extract the subset of A into B, then compute the P matrix as B\A
    which will automatically compute P using a suitable least-squares
    algorithm.

    The approximation we use here is to compute the column pivots of Y,
    rather then use the true column pivots as would be computed by a column-
    pivoted QR process.

Reference:

    \cite[Algorithm I]{Liberty2007}

    @article{Cheng2005,
        author = {Cheng, H and Gimbutas, Z and Martinsson, P G and Rokhlin, V},
        doi = {10.1137/030602678},
        issn = {1064-8275},
        journal = {SIAM Journal on Scientific Computing},
        month = jan,
        number = {4},
        pages = {1389--1404},
        title = {On the Compression of Low Rank Matrices},
        volume = {26},
        year = {2005}
    }
""" ->
function idfact(A, k::Int, l::Int)
    m, n = size(A)
    R = randn(l, m)
    Y = R * A #size l x n

    #Compute column pivots of first k columns of Y
    maxvals = map(j->maximum(abs(view(Y, :, j))), 1:n)
    piv = sortperm(maxvals, rev=true)[1:k]

    B = A[:, piv]
    Interpolative(B, B\A)
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RDatasets
using Random
using Statistics
using LinearAlgebra
using Plots

iris = dataset("datasets", "iris")
randn(10)
plot(1:10)

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
"""
Authors: Yaroslav Salii, 2020+
         Kara Ignatenko, 2021

Main file for Oboe. Putting all the processing together,
debugging, etc.

Eventually, there'll be a couple debug scenarios and some command-line args parsing

Run me by typing `include("oboe-main.jl").processOboe(...)` at the julia REPL

oboe-main.jl
2020-12-10  v.0.0 A Hello, World!
2021-02-02  v.0.1 Just the AP-AP to node-node processing
2021-02-19  v.0.5 Full processing, commented out; still no IO
2021-07-08  v.0.6 Implemented IO via CLI arguments
2021-07-12  v.0.7 Moved CLI parsing to a separate file
2021-04-14  v.0.8 Added support for tract datasets other than NW and checking
                  if the output already exists
2021-07-19  v.1.0 fixed issue where unmatched tract IDs led to a crash
2021-08-17  v.1.1 updated to support the new API of Oboe.jl and FromFile module import
2021-09-08  v.1.2 isolate aggregation/partition function pair dispatch
2021-10-13  v.1.3 automatically normalize single-digit FIPS codes without a leading zero
2021-11-23  v.1.4 move dispatchAggPart to aggregation.jl, closer to the funcs it dispatches
"""

module OboeMain

"The FIPS codes for each of the contiguous US states (and DC)"
fipsAll =  ["01", "04", "05", "06", "08", "09", "10", "11", "12", "13", "16", "17", "18",
            "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31",
            "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "44", "45",
            "46", "47", "48", "49", "50", "51", "53", "54", "55", "56"]

"Wash. + Ore, 2072 tracts"
fipsNW = ["41","53"]

"NW's first neighbors except Cal., 2834 tracts"
fipsNW1 = ["16","32","41","53"]

"NW's 1st and 2nd neighbors except Cal., 4830 tracts, 259 ctys"
fipsNW2 = ["04","16","30","32","41","49","53","56"]

"Cal., 7038 tracts, 58 ctys"
fipsCA = ["06"]

"West Coast, Cal. + Ore. + Wash., 9910 tracts, 133 ctys"
fipsWCT= ["06","41","53"]

using FromFile
using ArgParse
@from "./Oboe/Oboe.jl" import Oboe
using DataFrames

export processOboe

println(Oboe.callsign)

#lame logging thing
myshow = obj -> println(first(obj,5))

"blueprint for getProcessedAPs() method "
#Oboe.getProcessedAPs(Oboe.rdBTS() |> Oboe.grpBTS,Oboe.rdAPs())

function processOboe(name, agg; fips=fipsAll, useNW=false, force=false)
    #construct the output file name and make sure the files don't already exist
    iname = name * agg
    if !force && Oboe.checkIfFilesExist(iname)
        error("Output data for $iname already exists, use the force flag to overwrite")
    end
    println("Preparing the tract, air travel, and commute data...")
    fips = map(Oboe.normalizeOneDigitFIPS, fips) # convert FIPS codes such as "1" to "01"
    if useNW
        iFluteFile="a~NW" * Oboe.fn.sep * Oboe.fn.fltInitSuff
        #read the census tracts corresponding to the name provided in the args
        nsRaw = Oboe.rdFluteTract(iFluteFile)
    else
        wholeUS = Oboe.rdWholeUS()
        nsRaw = Oboe.censorFluteTractByFIPS(wholeUS, fips)
    end
    #nsRaw |> myshow

    #read flight data from BTS
    rawBTS = Oboe.rdBTS()
    #read Openflights.org airport data: [:IATA_Code,:LAT,:LNG,:Name]
    rawAPs = Oboe.rdAPs()

    #all AP-AP travel, as list o'pairs [:ORG,:DST,:PSG], :ORG and :DST are :IATA_Code
    pBTS = Oboe.grpBTS(rawBTS)
    #cache the smallest reasonable APs, [:IATA_Code,:LAT,:LNG,:IN.+:OUT ≥ 2500]
    APs = Oboe.getProcessedAPs(pBTS, rawAPs)
    #to each node, assign a designated AP from `APs` -> +[:IATA_Code]
    ns = Oboe.assignDsgAPs(nsRaw,APs)
    #now find the :Pop of each APs' catchment area, and chuck that into a `Dict`
    d = Oboe.mkAP_pop_dict(ns)
    #go on, compute the nodes' passenger shares (add the :shr col), with `d` in mind
    ns2 = Oboe.assignPsgShares(ns,d)
    #ns2 |> myshow

    #generate matrix of AP-to-AP flights with aux. arrays mapping indices to IATA codes
    fmx = Oboe.mkFlightMx2(pBTS, ns2.IATA_Code |> unique; daily=true)
    #generate commute table
    cmt = Oboe.rdTidyWfsByFIPS(fips, ns2)
    #cmt |> myshow

    #do the processing, aggregating by tract/state/county/ap as selected in the args

    println("Processing started. Aggregation mode [$agg] selected.")
    aggnpart = Oboe.dispatchAggPart(agg) #dispatch the aggregate/partition function pair
    println("Aggregating")
    @time aggregated = ns2 |> aggnpart[1] #aggregate as selected by dispatcher
    println("Partitioning")
    @time partitioned = aggnpart[2](ns2,aggregated) #partition as selected by dispatcher

    skipagg = (agg == "tra") ? true : false #no aggregation necessary if "tra" requested
    println("Computing aggregated commute matrix")
    @time cmtMx = skipagg ? Oboe.mkCmtMx(ns2, cmt) : Oboe.mkCmtMx(ns2, aggregated, partitioned, cmt)
    println("Computing aggregated air travel matrix")
    @time psgMx = skipagg ? Oboe.mkPsgMx(ns2, fmx) : Oboe.mkPsgMx(ns2, fmx, aggregated, partitioned)

    #infect2! is aware of all aggregation types
    iv  = Oboe.infect2!(Oboe.ns2iv_sterile(aggregated),aggType=agg)

    println("Generated instance name: ",iname,"_",nrow(iv))
    println("Writing -init.csv and -trav.dat")
    Oboe.writeMe(iname,iv,psgMx + cmtMx)
    println("All done. Terminating.")
end

end #end module OboeMain


text:
function betterstrings_strjoin(line::Ref{UInt}, ns::Namespace, arr::AetherlangObject, sep::AetherlangObject=AetherlangObject(""))::AetherlangObject
    AetherlangObject(join([arr.dataref[i] for i in 1:length(arr.dataref)], sep.dataref))
end
function betterstrings_strsplit(line::Ref{UInt}, ns::Namespace, str::AetherlangObject, sep::AetherlangObject=AetherlangObject(" "))::AetherlangObject
    AetherlangObject(split(str.dataref, sep.dataref))
end
function betterstrings_charemplace(line::Ref{UInt}, ns::Namespace, str::AetherlangObject, indx::AetherlangObject, char::AetherlangObject)::AetherlangObject
    if length(char.dataref) > 1
        throw(AetherlangError("$(repr(char.dataref)) is not a single character"))
    end
    if indx.dataref == 1
        return AetherlangObject(char.dataref*str.dataref[2:end])
    elseif indx.dataref == length(str.dataref)
        return AetherlangObject(str.dataref[1:end-1]*char.dataref)
    end
    AetherlangObject(str.dataref[1:indx.dataref-1]*char.dataref*str.dataref[indx.dataref+1:end])
end
betterstrings_tolowercase(line::Ref{UInt}, ns::Namespace, str::AetherlangObject)::AetherlangObject = AetherlangObject(lowercase(str.dataref))
betterstrings_touppercase(line::Ref{UInt}, ns::Namespace, str::AetherlangObject)::AetherlangObject = AetherlangObject(uppercase(str.dataref))
function betterstrings_digit(line::Ref{UInt}, ns::Namespace, str::AetherlangObject)::AetherlangObject
    for c in str.dataref
        if !isdigit(c)
            return AetherlangObject()
        end
    end
    AetherlangObject(1)
end
function betterstrings_punct(line::Ref{UInt}, ns::Namespace, str::AetherlangObject)::AetherlangObject
    for c in str.dataref
        if !ispunct(c)
            return AetherlangObject()
        end
    end
    AetherlangObject(1)
end
function betterstrings_printable(line::Ref{UInt}, ns::Namespace, str::AetherlangObject)::AetherlangObject
    for c in str.dataref
        if !isprint(c)
            return AetherlangObject()
        end
    end
    AetherlangObject(1)
end
function betterstrings_whitespace(line::Ref{UInt}, ns::Namespace, str::AetherlangObject)::AetherlangObject
    for c in str.dataref
        if !isspace(c)
            return AetherlangObject()
        end
    end
    AetherlangObject(1)
end

betterstrings_namespace_modify = Namespace(
    "strjoin" => AetherlangObject(betterstrings_strjoin),
    "strsplit" => AetherlangObject(betterstrings_strsplit),
    "charemplace" => AetherlangObject(betterstrings_charemplace),
    "tolowercase" => AetherlangObject(betterstrings_tolowercase),
    "touppercase" => AetherlangObject(betterstrings_touppercase),
    "digit?" => AetherlangObject(betterstrings_digit),
    "punct?" => AetherlangObject(betterstrings_punct),
    "printable?" => AetherlangObject(betterstrings_printable),
    "whitespace?" => AetherlangObject(betterstrings_whitespace)
)
betterstrings_namespace_modify["charemp"] = betterstrings_namespace_modify["charemplace"]

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
module BPFnative

using LLVM, LLVM.Interop
using InteractiveUtils
using Libdl
using GPUCompiler
using CBinding
using Preferences

# Kernel Internal API
const use_vmlinux = parse(Bool, @load_preference("use_vmlinux", "false"))
if use_vmlinux
    const has_vmlinux = try
        using VMLinuxBindings
        c"VMLinuxBindings.struct pt_regs"
        true
    catch err
        @warn "Failed to load Linux Kernel definitions: $err"
        false
    end
else
    const has_vmlinux = false
end
function enable_vmlinux(ans::Bool=true)
    ans_str = repr(ans)
    if ans_str != repr(use_vmlinux)
        @set_preferences!("use_vmlinux"=>ans_str)
        @info "$(ans ? "Enabled" : "Disabled") loading vmlinux; please restart Julia for changes to take effect"
    end
end

# Common API
module API
if !parse(Bool, get(ENV, "JULIA_BPFNATIVE_DISABLE_ARTIFACTS", "0"))
    using Libbpf_jll
else
    const libbpf = "libbpf"
end
import ..BPFnative: has_vmlinux
import ..CBinding
import ..CBinding: @c_str, Cptr
if has_vmlinux
    import ..VMLinuxBindings
    const VMLinux = VMLinuxBindings
end
include("utils.jl")
include("common.jl")
include("libbpf.jl")
if Sys.islinux()
include("network.jl")
end
end

# Runtime API
module RT
import ..API
using ..LLVM
using ..LLVM.Interop
import Core: LLVMPtr
import ..CBinding
import ..CBinding: @c_str, Cptr
include("runtime/maps_core.jl")
include("runtime/bpfcall.jl")
include("runtime/maps.jl")
include("runtime/buffers.jl")
include("runtime/helpers.jl")
include("runtime/intrinsics.jl")
include("runtime/constants.jl")
include("runtime/utils.jl")
end

# Host API
module Host
import ..API
include("host/syscall.jl")
include("host/maps.jl")
include("host/socket.jl")
include("host/kallsyms.jl")
end

# Compiler
include("compiler.jl")
include("reflection.jl")

# Easy-to-use probes
include("probes.jl")

# Useful utilities
include("extra.jl")
include("xdp.jl")

end # module


text:
# The order here is designed to avoid an ambiguity warning in convert,
# see the top of ImageAxes
using ImageAxes
using ImageCore
using Images
import Images
using Printf
import AxisArrays
using OffsetArrays

import Base: +, -, *, /
import Base: permutedims

const ViewIndex = Union{Base.ViewIndex, Colon}

export
    # types
    DirectImage,

    # functions
    copyheaders,
    arraydata,
    headers,
    shareheaders,
    spatialheaders

#### types and constructors ####

using OrderedCollections

# Valid types for FITS headers. Restrict to only creating these values
# instead of erroring when trying to write them out.
ValTypes = Union{String, Bool, Integer, AbstractFloat, Nothing}

const FLEN_COMMENT = 73  # max length of a keyword comment string

struct CommentedValue
    value::ValTypes
    comment::String
    function CommentedValue(value, comment)
        if sizeof(comment) > FLEN_COMMENT
                @warn "Truncating comment to 73 bytes"
                comment = comment[begin:begin+FLEN_COMMENT]
        end
        new(value,comment)
    end
end

# Concrete types
"""
`DirectImage` is an AbstractArray that can have metadata, stored in a dictionary.
Construct an image with `DirectImage(A, props)` (for a headers dictionary
`props`), or with `DirectImage(A, prop1=val1, prop2=val2, ...)`.
"""
struct DirectImage{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    data::A
    headers::OrderedDict{Symbol,CommentedValue}
end

# TODO: friendly constructors

function DirectImage(data::AbstractArray{T,N}) where {T,N}
    return DirectImage{T,N,typeof(data)}(data, OrderedDict{Symbol,CommentedValue}())
end

function DirectImage(data::AbstractArray{T,N}, headers::OrderedDict{Symbol,CommentedValue}) where {T,N}
    return DirectImage{T,N,typeof(data)}(data, headers)
end

function DirectImage(data::AbstractArray{T,N}, headers::OrderedDict{Symbol,CommentedValue}, axes::AbstractRange...) where {T,N}
    data_offset = OffsetArray(data, axes...)
    return DirectImage{T,N,typeof(data_offset)}(data_offset, headers)
end

const DirectImageArray{T,N,A<:Array} = DirectImage{T,N,A}
const DirectImageAxis{T,N,A<:AxisArray} = DirectImage{T,N,A}

Base.size(A::DirectImage) = size(arraydata(A))
Base.size(A::DirectImageAxis, Ax::Axis) = size(arraydata(A), Ax)
Base.size(A::DirectImageAxis, ::Type{Ax}) where {Ax<:Axis} = size(arraydata(A), Ax)
Base.axes(A::DirectImage) = axes(arraydata(A))
Base.axes(A::DirectImageAxis, Ax::Axis) = axes(arraydata(A), Ax)
Base.axes(A::DirectImageAxis, ::Type{Ax}) where {Ax<:Axis} = axes(arraydata(A), Ax)

"""
    origin(image_or_array)

Return the index corresponding to the "origin" of the data. If a regular array
or DirectImage, this is the same as `Images.center`.

This function diverges from `Images.center` if the data has offset indices. In that
case, it returns the index closest to zero.

Example:
julia> rand(3,3) |> origin

"""
function origin(A::AbstractArray)
    Ifirst = first(CartesianIndices(A))
    if Ifirst == CartesianIndex((1 for _ in 1:length(Ifirst))...)
        return Images.center(A)
    else
        # We should be able to find this without iterating through all indices

        # How to find the index closest to zero?
        # CartesianIndex(Images.center.(axes(A))...)
        indices = Int[]
        for ax in 1:length(size(A))
            push!(indices, argmin(abs.(UnitRange(axes(A,ax)))))
        end
        return indices

        # closest_to_zero = firstindex(CartesianIndices(a))
        # closest_dist = zero(eltype(closest_to_zero))
        # for I in CartesianIndices(a)
        #     squares = zero(eltype(I))
        #     for i in eachindex(I)
        #         squares += i^2
        #     end

    end
end
function origin(A::AbstractArray, axis)
    return origin(A)[axis]
end
export origin

datatype(::Type{DirectImage{T,N,A}}) where {T,N,A<:AbstractArray} = A

Base.IndexStyle(::Type{M}) where {M<:DirectImage} = IndexStyle(datatype(M))

AxisArrays.HasAxes(A::DirectImageAxis) = AxisArrays.HasAxes{true}()

# getindex and setindex!
for AType in (DirectImage, DirectImageAxis)
    @eval begin
        @inline function Base.getindex(img::$AType{T,1}, i::Int) where T
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds ret = arraydata(img)[i]
            ret
        end
        @inline function Base.getindex(img::$AType, i::Int)
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds ret = arraydata(img)[i]
            ret
        end
        @inline function Base.getindex(img::$AType{T,N}, I::Vararg{Int,N}) where {T,N}
            @boundscheck checkbounds(arraydata(img), I...)
            @inbounds ret = arraydata(img)[I...]
            ret
        end

        @inline function Base.setindex!(img::$AType{T,1}, val, i::Int) where T
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds arraydata(img)[i] = val
            val
        end
        @inline function Base.setindex!(img::$AType, val, i::Int)
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds arraydata(img)[i] = val
            val
        end
        @inline function Base.setindex!(img::$AType{T,N}, val, I::Vararg{Int,N}) where {T,N}
            @boundscheck checkbounds(arraydata(img), I...)
            @inbounds arraydata(img)[I...] = val
            val
        end
    end
end

@inline function Base.getindex(img::DirectImageAxis, ax::Axis, I...)
    result = arraydata(img)[ax, I...]
    maybe_wrap(img, result)
end
@inline function Base.getindex(img::DirectImageAxis, i::Union{Integer,AbstractVector,Colon}, I...)
    result = arraydata(img)[i, I...]
    maybe_wrap(img, result)
end
maybe_wrap(img::DirectImage{T}, result::T) where T = result
maybe_wrap(img::DirectImage{T}, result::AbstractArray{T}) where T = copyheaders(img, result)

@inline function Base.setindex!(img::DirectImageAxis, val, ax::Axis, I...)
    setindex!(arraydata(img), val, ax, I...)
end
@inline function Base.setindex!(img::DirectImageAxis, val, i::Union{Integer,AbstractVector,Colon}, I...)
    setindex!(arraydata(img), val, i, I...)
end

Base.view(img::DirectImage, ax::Axis, I...) = shareheaders(img, view(arraydata(img), ax, I...))
Base.view(img::DirectImage{T,N}, I::Vararg{ViewIndex,N}) where {T,N} = shareheaders(img, view(arraydata(img), I...))
Base.view(img::DirectImage, i::ViewIndex) = shareheaders(img, view(arraydata(img), i))
Base.view(img::DirectImage, I::Vararg{ViewIndex,N}) where {N} = shareheaders(img, view(arraydata(img), I...))

@inline function Base.getproperty(img::DirectImage, propname::Symbol)::ValTypes
    if !haskey(headers(img), propname)
        propname = Symbol(uppercase(string(propname)))
    end
    return headers(img)[propname].value
end
@inline function Base.getindex(img::DirectImage, propname::Union{Symbol,String})::ValTypes
    if !haskey(headers(img), propname)
        propname = Symbol(uppercase(string(propname)))
    end
    return headers(img)[propname].value
end
@inline function Base.getindex(img::DirectImage, propname::Union{Symbol,String}, ::typeof(/))::String
    if !haskey(headers(img), propname)
        propname = Symbol(uppercase(string(propname)))
    end
    return headers(img)[propname].comment
end
@inline function Base.setproperty!(img::DirectImage, propname::Symbol, val)
    if !haskey(headers(img), propname)
        propname = Symbol(uppercase(string(propname)))
    end
    if haskey(headers(img), propname)
        headers(img)[propname] = CommentedValue(val, headers(img)[propname].comment)
    else
        headers(img)[propname] = CommentedValue(val, "")
    end
    return img
end
@inline function Base.setindex!(img::DirectImage, val, propname::Union{Symbol,String})
    if !haskey(headers(img), propname)
        propname = Symbol(uppercase(string(propname)))
    end
    if haskey(headers(img), propname)
        headers(img)[propname] = CommentedValue(val, headers(img)[propname].comment)
    else
        headers(img)[propname] = CommentedValue(val, "")
    end
    return img
end
@inline function Base.setindex!(img::DirectImage, comment, propname::Union{Symbol,String}, ::typeof(/))
    propname = Symbol(uppercase(string(propname)))
    if haskey(headers(img), propname)
        headers(img)[propname] = CommentedValue(img[propname], string(comment))
    else
        headers(img)[propname] = CommentedValue(nothing, string(comment))
    end
    comment
end

Base.propertynames(img::DirectImage) = (keys(headers(img))...,)

Base.copy(img::DirectImage) = DirectImage(copy(arraydata(img)), deepcopy(headers(img)))

Base.convert(::Type{DirectImage}, A::DirectImage) = A
Base.convert(::Type{DirectImage}, A::AbstractArray) = DirectImage(A)
Base.convert(::Type{DirectImage{T}}, A::DirectImage{T}) where {T} = A
Base.convert(::Type{DirectImage{T}}, A::DirectImage) where {T} = shareheaders(A, convert(AbstractArray{T}, arraydata(A)))
Base.convert(::Type{DirectImage{T}}, A::AbstractArray{T}) where {T} = DirectImage(A)
Base.convert(::Type{DirectImage{T}}, A::AbstractArray) where {T} = DirectImage(convert(AbstractArray{T}, A))

# copy headers
function Base.copy!(imgdest::DirectImage, imgsrc::DirectImage, prop1::Symbol, props::Symbol...)
    setproperty!(imgdest, prop1, deepcopy(getproperty(imgsrc, prop1)))
    for p in props
        setproperty!(imgdest, p, deepcopy(getproperty(imgsrc, p)))
    end
    return imgdest
end

# similar
Base.similar(img::DirectImage, ::Type{T}, shape::Dims) where {T} = DirectImage(similar(arraydata(img), T, shape), copy(headers(img)))
Base.similar(img::DirectImageAxis, ::Type{T}) where {T} = DirectImage(similar(arraydata(img), T), copy(headers(img)))
Base.similar(img::DirectImage{T,N,<:OffsetArray}) where{T,N} = DirectImage(similar(arraydata(img), T), copy(headers(img)))

"""
    copyheaders(img::DirectImage, data) -> imgnew
Create a new "image," copying the headers dictionary of `img` but
using the data of the AbstractArray `data`. Note that changing the
headers of `imgnew` does not affect the headers of `img`.
See also: [`shareheaders`](@ref).
"""
copyheaders(img::DirectImage, data::AbstractArray) =
    DirectImage(data, copy(headers(img)))

"""
    shareheaders(img::DirectImage, data) -> imgnew
Create a new "image," reusing the headers dictionary of `img` but
using the data of the AbstractArray `data`. The two images have
synchronized headers; modifying one also affects the other.
See also: [`copyheaders`](@ref).
"""
shareheaders(img::DirectImage, data::AbstractArray) = DirectImage(data, headers(img))

# Delete a property!
Base.delete!(img::DirectImage, propname::Symbol) = delete!(headers(img), propname)

# Iteration
# Defer to the array object in case it has special iteration defined
Base.iterate(img::DirectImage) = Base.iterate(arraydata(img))
Base.iterate(img::DirectImage, s) = Base.iterate(arraydata(img), s)

# Show
const emptyset = Set()
function showim(io::IO, img::DirectImage)
    IT = typeof(img)
    print(io, eltype(img).name.name, " DirectImage with:\n  data: ", summary(arraydata(img)), "\n  headers")
    showdictlines(io, headers(img), get(img, :suppress, emptyset))
end
Base.show(io::IO, img::DirectImage) = showim(io, img)
Base.show(io::IO, ::MIME"text/plain", img::DirectImage) = showim(io, img)

function Base.reinterpret(::Type{T}, img::DirectImage) where {T}
    shareheaders(img, reinterpret(T, arraydata(img)))
end
if Base.VERSION >= v"1.6.0-DEV.1083"
    function Base.reinterpret(::typeof(reshape), ::Type{T}, img::DirectImage) where {T}
        shareheaders(img, reinterpret(reshape, T, arraydata(img)))
    end
end

"""
    arraydata(img::DirectImage) -> array
Extract the data from `img`, omitting the headers
dictionary. `array` shares storage with `img`, so changes to one
affect the other.
See also: [`headers`](@ref).
"""
ImageAxes.arraydata(img::DirectImage) = getfield(img, :data)

function Base.PermutedDimsArray(A::DirectImage, perm)
    ip = sortperm([perm...][[coords_spatial(A)...]])  # the inverse spatial permutation
    permutedims_props!(copyheaders(A, PermutedDimsArray(arraydata(A), perm)), ip)
end
ImageCore.channelview(A::DirectImage) = shareheaders(A, channelview(arraydata(A)))
ImageCore.rawview(A::DirectImage{T}) where {T<:Real} = shareheaders(A, rawview(arraydata(A)))
ImageCore.normedview(::Type{T}, A::DirectImage{S}) where {T<:FixedPoint,S<:Unsigned} = shareheaders(A, normedview(T, arraydata(A)))

# AxisArrays functions
AxisArrays.axes(img::DirectImageAxis) = AxisArrays.axes(arraydata(img))
AxisArrays.axes(img::DirectImageAxis, d::Int) = AxisArrays.axes(arraydata(img), d)
AxisArrays.axes(img::DirectImageAxis, Ax::Axis) = AxisArrays.axes(arraydata(img), Ax)
AxisArrays.axes(img::DirectImageAxis, ::Type{Ax}) where {Ax<:Axis} = AxisArrays.axes(arraydata(img), Ax)
AxisArrays.axisdim(img::DirectImageAxis, ax) = axisdim(arraydata(img), ax)
AxisArrays.axisnames(img::DirectImageAxis) = axisnames(arraydata(img))
AxisArrays.axisvalues(img::DirectImageAxis) = axisvalues(arraydata(img))

#### Properties ####

"""
    headers(img) -> props
Extract the headers dictionary `props` for `img`. `props`
shares storage with `img`, so changes to one affect the other.
See also: [`arraydata`](@ref).
"""
headers(img::DirectImage) = getfield(img, :headers)

import Base: hasproperty
hasproperty(img::DirectImage, k::Symbol) = haskey(headers(img), k)

Base.get(img::DirectImage, k::Symbol, default) = get(headers(img), k, default)

# So that defaults don't have to be evaluated unless they are needed,
# we also define a @get macro (thanks Toivo Hennington):
struct IMNothing end   # to avoid confusion in the case where dict[key] === nothing
macro get(img, k, default)
    quote
        img, k = $(esc(img)), $(esc(k))
        val = get(headers(img), k, IMNothing())
        return isa(val, IMNothing) ? $(esc(default)) : val
    end
end

ImageAxes.timeaxis(img::DirectImageAxis) = timeaxis(arraydata(img))
ImageAxes.timedim(img::DirectImageAxis) = timedim(arraydata(img))

ImageCore.pixelspacing(img::DirectImage) = pixelspacing(arraydata(img))

"""
    spacedirections(img)
Using DirectImagedata, you can set this property manually. For example, you
could indicate that a photograph was taken with the camera tilted
30-degree relative to vertical using
```
img["spacedirections"] = ((0.866025,-0.5),(0.5,0.866025))
```
If not specified, it will be computed from `pixelspacing(img)`, placing the
spacing along the "diagonal".  If desired, you can set this property in terms of
physical units, and each axis can have distinct units.
"""
ImageCore.spacedirections(img::DirectImage) = @get img :spacedirections spacedirections(arraydata(img))

ImageCore.sdims(img::DirectImageAxis) = sdims(arraydata(img))

ImageCore.coords_spatial(img::DirectImageAxis) = coords_spatial(arraydata(img))

ImageCore.spatialorder(img::DirectImageAxis) = spatialorder(arraydata(img))

ImageAxes.nimages(img::DirectImageAxis) = nimages(arraydata(img))

ImageCore.size_spatial(img::DirectImageAxis) = size_spatial(arraydata(img))

ImageCore.indices_spatial(img::DirectImageAxis) = indices_spatial(arraydata(img))

ImageCore.assert_timedim_last(img::DirectImageAxis) = assert_timedim_last(arraydata(img))

#### Permutations over dimensions ####

"""
    permutedims(img, perm, [spatialprops])
When permuting the dimensions of an DirectImage, you can optionally
specify that certain headers are spatial and they will also be
permuted. `spatialprops` defaults to `spatialproperties(img)`.
"""
permutedims

function permutedims_props!(ret::DirectImage, ip, spatialprops=spatialproperties(ret))
    if !isempty(spatialprops)
        for prop in spatialprops
            if hasproperty(ret, prop)
                a = getproperty(ret, prop)
                if isa(a, AbstractVector)
                    setproperty!(ret, prop, a[ip])
                elseif isa(a, Tuple)
                    setproperty!(ret, prop, a[ip])
                elseif isa(a, AbstractMatrix) && size(a,1) == size(a,2)
                    setproperty!(ret, prop, a[ip,ip])
                else
                    error("Do not know how to handle property ", prop)
                end
            end
        end
    end
    ret
end

function permutedims(img::DirectImageAxis, perm)
    p = AxisArrays.permutation(perm, axisnames(arraydata(img)))
    ip = sortperm([p...][[coords_spatial(img)...]])
    permutedims_props!(copyheaders(img, permutedims(arraydata(img), p)), ip)
end
function permutedims(img::DirectImage, perm)
    ip = sortperm([perm...][[coords_spatial(img)...]])
    permutedims_props!(copyheaders(img, permutedims(arraydata(img), perm)), ip)
end

# Note: `adjoint` does not recurse into DirectImage properties.
function Base.adjoint(img::DirectImage{T,2}) where {T<:Real}
    ip = sortperm([2,1][[coords_spatial(img)...]])
    permutedims_props!(copyheaders(img, adjoint(arraydata(img))), ip)
end

function Base.adjoint(img::DirectImage{T,1}) where T<:Real
    check_empty_spatialproperties(img)
    copyheaders(img, arraydata(img)')
end

"""
    spatialproperties(img)
Return a vector of strings, containing the names of properties that
have been declared "spatial" and hence should be permuted when calling
`permutedims`.  Declare such properties like this:
    img[:spatialproperties] = [:spacedirections]
"""
spatialproperties(img::DirectImage) = @get img :spatialproperties [:spacedirections]

function check_empty_spatialproperties(img)
    sp = spatialproperties(img)
    for prop in sp
        if hasproperty(img, prop)
            error("spatialproperties must be empty, have $prop")
        end
    end
    nothing
end

function Images.centered(A::DirectImage)
    new_ax =  map(i->axes(A,i).-floor(Int, origin(A,i)), eachindex(axes(A)))
    DirectImage(A, headers(A), new_ax...)
end
export centered

#### Low-level utilities ####
function showdictlines(io::IO, dict::AbstractDict, suppress::Set)
    i = 0
    for (k, v) in dict
        i += 1
        if k === :suppress
            continue
        end
        if !in(k, suppress)
            val = string(v.value)[1:min(20,end)]
            @printf(io, "\n    %-8s= %-20s / %s", k, val, v.comment)
        else
            print(io, "\n    ", k, ": <suppressed>")
        end
        if i > 10
            l = length(dict) - 10
            print(io, "\n    ...       $l additional headers redacted")
            break
        end
    end
end
showdictlines(io::IO, dict::AbstractDict, prop::Symbol) = showdictlines(io, dict, Set([prop]))

max_stars_repo_name: sefffal/SpaceImages


text:
module WagnerModel

include("simFunctions.jl")

export generate_config
export assess_stability
# Write your package code here.

end


text:
using ObservationSchemes
using Test
using StaticArrays, LinearAlgebra
const OBS = ObservationSchemes

@testset "OBS.jl" begin
    # Write your own tests here.
    @testset "default constructor" begin
        v = @SVector [1.0, 2.0, 3.0]
        tt = 1.0
        obs = LinearGsnObs(tt, v)
        @test obs.t == tt
        @test eltype(obs) == typeof(v)
        @test size(obs) == length(v)
        @test length(obs) == length(v)
        @test fpt_info(obs) == NoFirstPassageTimes
        @test obs.obs == v
        @test obs.L == SDiagonal{3}(I)
        @test obs.μ == zero(v)
        @test obs.Σ == SDiagonal{3}(I)*1e-11
        @test obs.full_obs == false
    end

    @testset "custom constructor 1" begin
        v = @SVector [1.0, 2.0, 3.0]
        tt = 2.0
        L = @SMatrix [1.0 2.0 0.0 4.0; 3.0 4.0 0.0 0.0; 0.0 1.0 2.0 1.0]
        Σ = SDiagonal{3}(I)
        fpt = FirstPassageTimeInfo(
            (1, 2),
            (1.0, 2.0),
            (true, false),
            (false, true),
            (1.0,),
        )
        obs = LinearGsnObs(tt, v; L = L, Σ = Σ, fpt = fpt)

        @test obs.t == tt
        @test eltype(obs) == typeof(v)
        @test size(obs) == length(v)
        @test length(obs) == length(v)
        @test fpt_info(obs) == typeof(fpt)
        @test obs.obs == v
        @test obs.L == L
        @test obs.μ == zero(v)
        @test obs.Σ == Σ
        @test obs.full_obs == false
    end
end

@testset "Container for all observations" begin
    struct LawA α; β; end
    struct LawB γ; β; end
    OBS.var_parameter_names(P::LawA) = [:α, :β]
    OBS.var_parameter_names(P::LawB) = [:γ, :β]

    all_obs = AllObservations()
    P, obs, t, x0_prior = LawA(10,20), [1.0, 2.0], 0.0, KnownStartingPt(0.0)
    P2, obs2, t2, x0_prior2 = LawB(30,40), [11.0, 12.0], 1.0, KnownStartingPt(2.0)

    @testset "recording observations" begin
        add_recording!(all_obs, (P=P, obs=obs, t0=t, x0_prior=x0_prior))
        @test all_obs.recordings[1].P == P
        @test all_obs.recordings[1].obs == obs
        @test all_obs.recordings[1].t0 == t
        @test all_obs.recordings[1].x0_prior == x0_prior

        add_recording!(all_obs, (P=P2, obs=obs2, t0=t2, x0_prior=x0_prior2))
        @test all_obs.recordings[1].P == P
        @test all_obs.recordings[1].obs == obs
        @test all_obs.recordings[1].t0 == t
        @test all_obs.recordings[1].x0_prior == x0_prior
        @test all_obs.recordings[2].P == P2
        @test all_obs.recordings[2].obs == obs2
        @test all_obs.recordings[2].t0 == t2
        @test all_obs.recordings[2].x0_prior == x0_prior2
    end

    @testset "adding dependency structure" begin
        all_obs = AllObservations()
        add_recording!(all_obs, (P=P, obs=obs, t0=t, x0_prior=x0_prior))
        add_recording!(all_obs, (P=P2, obs=obs2, t0=t2, x0_prior=x0_prior2))
        add_dependency!(all_obs, Dict(:β => [(1, :β), (2, :β)]))
        @test all_obs.param_depend[:β] == [1 => :β, 2 => :β]

        all_obs = AllObservations()
        add_recording!(all_obs, (P=P, obs=obs, t0=t, x0_prior=x0_prior))
        add_recording!(all_obs, (P=P2, obs=obs2, t0=t2, x0_prior=x0_prior2))
        add_dependency!(all_obs, Dict(:β => [(1, :β), (2, :β)]))
        @test all_obs.param_depend[:β] == [1 => :β, 2 => :β]
    end

    @testset "splitting at full observations" begin
        obs = [
            [
                LinearGsnObs(1.0, 1.0; full_obs=false),
                LinearGsnObs(2.0, 2.0; full_obs=true),
                LinearGsnObs(3.0, 3.0; full_obs=false),
            ],
            [
                LinearGsnObs(1.3, 1.0; full_obs=true),
                LinearGsnObs(2.3, 2.0; full_obs=true),
                LinearGsnObs(3.3, 3.0; full_obs=true),
            ],
            [
                LinearGsnObs(1.5, 1.0; full_obs=false),
                LinearGsnObs(2.5, 2.0; full_obs=false),
                LinearGsnObs(3.5, 3.0; full_obs=false),
                LinearGsnObs(4.5, 4.0; full_obs=false),
                LinearGsnObs(5.5, 5.0; full_obs=false),
            ],
        ]
        recordings = [
            (
                P=LawA(10,20),
                obs=obs[1],
                t0=0.0,
                x0_prior=KnownStartingPt(2.0),
            ),
            (
                P=LawA(10,20),
                obs=obs[2],
                t0=0.3,
                x0_prior=KnownStartingPt(-2.0),
            ),
            (
                P=LawB(30,40),
                obs=obs[3],
                t0=0.5,
                x0_prior=KnownStartingPt(10.0),
            ),
        ]
        all_obs = AllObservations()
        for recording in recordings
            add_recording!(all_obs, recording)
        end
        add_dependency!(
            all_obs,
            Dict(
                :α => [(1, :α), (2, :α)],
                :β => [(1, :β), (2, :β), (3, :β)],
            )
        )
        d = all_obs.param_depend
        @testset "pre-test α dependence $i" for i in 1:2
            @test d[:α][i] == (i => :α)
        end
        @testset "pre-test β dependence $i" for i in 1:3
            @test d[:β][i] == (i => :β)
        end
        @testset "pre-test presence of recordings $i" for i in 1:3
            @test all_obs.recordings[i] == recordings[i]
        end

        initialized_obs, old_to_new_dict = initialize(all_obs)
        d = initialized_obs.param_depend
        @testset "new α dependence $i" for i in 1:4
            @test d[:α][i] == (i => :α)
        end
        @testset "new β dependence $i" for i in 1:6
            @test d[:β][i] == (i => :β)
        end
        @test d[:REC3_γ][1] == (6 => :γ)

        @test length(initialized_obs.recordings) == 6

        @testset "initialized recording 1" begin
            r = initialized_obs.recordings[1]
            @test r.P == recordings[1].P
            #TODO figure out why wrong
            #@test r.obs == recordings[1].obs[1:2]
            @test r.t0 == recordings[1].t0
            @test r.x0_prior == recordings[1].x0_prior
        end
        @testset "initialized recording 2" begin
            r = initialized_obs.recordings[2]
            @test r.P == recordings[1].P
            #@test r.obs == recordings[1].obs[3:3]
            @test r.t0 == recordings[1].obs[2].t
            @test r.x0_prior == KnownStartingPt(recordings[1].obs[2].obs)
        end
        @testset "initialized recording 3" begin
            r = initialized_obs.recordings[3]
            @test r.P == recordings[2].P
            #@test r.obs == recordings[2].obs[1:1]
            @test r.t0 == recordings[2].t0
            @test r.x0_prior == recordings[2].x0_prior
        end
        @testset "initialized recording 4" begin
            r = initialized_obs.recordings[4]
            @test r.P == recordings[2].P
            #@test r.obs == recordings[2].obs[2:2]
            @test r.t0 == recordings[2].obs[1].t
            @test r.x0_prior == KnownStartingPt(recordings[2].obs[1].obs)
        end
        @testset "initialized recording 5" begin
            r = initialized_obs.recordings[5]
            @test r.P == recordings[2].P
            #@test r.obs == recordings[2].obs[3:3]
            @test r.t0 == recordings[2].obs[2].t
            @test r.x0_prior == KnownStartingPt(recordings[2].obs[2].obs)
        end
        @testset "initialized recording 6" begin
            r = initialized_obs.recordings[6]
            @test r.P == recordings[3].P
            #@test r.obs == recordings[3].obs
            @test r.t0 == recordings[3].t0
            @test r.x0_prior == recordings[3].x0_prior
        end

        @testset "packaging functions" begin
            struct Foo end
            struct Foo2 end
            struct Foo3 end
            packaged1 = package(Foo(), all_obs)
            packaged2 = package(Foo, initialized_obs)
            packaged3 = package([Foo(), Foo2(), Foo3()], old_to_new_dict, initialized_obs)
            packaged4 = package([Foo(), Foo2(), Foo3()], all_obs)
            @test packaged1 == [
                [Foo(), Foo(), Foo()],
                [Foo(), Foo(), Foo()],
                [Foo(), Foo(), Foo(), Foo(), Foo()],
            ]
            @test packaged2 == [
                [Foo, Foo],
                [Foo],
                [Foo],
                [Foo],
                [Foo],
                [Foo, Foo, Foo, Foo, Foo],
            ]
            @test packaged3 == [
                [Foo(), Foo()],
                [Foo()],
                [Foo2()],
                [Foo2()],
                [Foo2()],
                [Foo3(), Foo3(), Foo3(), Foo3(), Foo3()],
            ]
            @test packaged4 == [
                [Foo(), Foo(), Foo()],
                [Foo2(), Foo2(), Foo2()],
                [Foo3(), Foo3(), Foo3(), Foo3(), Foo3()],
            ]
        end

    end
end


text:
struct SampleClassEnumeration{D,M}
  class::Vector{Vector{Int}}
  composition::Vector{Vector{Int}}
  pairs::Vector{Tuple{Int,Int}}
  reactant_to_class::OrderedDict{Tuple{Int,Int},Int}
  pair_to_classes::OrderedDict{Tuple{Int,Int},Vector{Int}}
  isactive::Vector{Bool}
  dummy_composition::Vector{Int}
end

function SampleClassEnumeration{D,M}(composition, pairs, reactant_to_class, pair_to_classes, isactive) where {D,M}
  number_pairwise = count(x -> x[2] != 0, pairs)
  number_onsite   = count(x -> x[2] == 0, pairs)
  number_classes  = capacity(M(), D) * number_pairwise + number_onsite

  class = [Int[] for s in 1:number_classes]
  dummy_composition = zero(composition[1])

  SampleClassEnumeration{D,M}(
    class,
    composition,
    pairs,
    reactant_to_class,
    pair_to_classes,
    isactive,
    dummy_composition
  )
end

function build_reactant_pairs(reactions)
  # so, we need to build up the number of sample classes while forming the reaction channels.
  # this involves counting the number of unique reactants, singly for on-site, in pairs for pairwise.
  # ordering is important for the pairs too (though we can fix this later if we want...)
  pairs = Tuple{Int,Int}[]

  for reaction in reactions
    # gather up the pairs of reactants:
    push!(pairs, (reaction.input1, reaction.input2) )
  end

  # this returns the unique pairs of reactant species.  On-site reactions are denoted by the second reactant being `0`.
  pairs = unique(pairs)
end

function map_reactant_to_class(pairs, nbmax)
  # Now that we've got the pairs, we need to form the sample indices.
  # This is equivalent to enumerating the different
  # unique pairs of reactants, splitting the pairwise reactants into 2*d groups (1 for each possible number of
  # adjacent product particles) and leaving one group for each on-site reactant.  Each reaction (old reaction)
  # points to exactly one of these groups, the group that it samples from.

  # we need a way to turn a reactant pair (plus number of adjacent members of the second type) into a sample class
  # so we turn to Dictionaries.

  # probably ought to sizehint this too...
  # and add type restrictions...
  reactant_to_class = OrderedDict{Tuple{Int,Int},Int}()

  ### I'm going to do this the janky way, first thing that came to mind.  I'm sure there's a better way of doing this.

  counts = 1
  for pair in pairs
    # add the pair and the referenced sample class to the dictionary
    reactant_to_class[pair] = counts

    # check to see if the pair is pairwise or on-site
    # the second entry is `0` if the pair is on-site
    if pair[2] == 0
      # we only need to increase the count by 1, since only one
      # sampling index is needed
      counts += 1

      # if the pair is pairwise, need to increase the count by
      # nbmax, one for each possible number of neighbor reactants
    else
      counts += nbmax
    end
  end

  return reactant_to_class
end

function build_active_set(pairs, number_types)
  isactive = zeros(Bool, number_types + 1)

  for pair in pairs
    reactant = pair[1]
    isactive[reactant] = true
  end

  return isactive
end

function map_pair_to_classes(composition, reactant_to_class, isactive, number_types)
  pair_to_classes = OrderedDict{Tuple{Int,Int},Vector{Int}}()
  number_compositions = length(composition)

  for l in 1:(number_types + 1)
    # check if the species is active before doing anything
    if isactive[l]
      for k in 1:number_compositions
        ### get the composition associated with k, .
        comp = composition[k]

        # get an empty vector to collect the sample classes a particle of type l w/ neighborhood k takes part in
        sample_classes = Int[]

        # iterate over every reactant pair (l, j)
        iterate_reactant_pairs!(sample_classes, l, comp, reactant_to_class, number_types)

        # finally, put the sample_classes in the dictionary with key (l, k)
        pair_to_classes[(l, k)] = sample_classes
      end
    end
  end
  return pair_to_classes
end

function iterate_reactant_pairs!(sample_classes, l, comp, reactant_to_class, number_types)
  for j in 1:(number_types + 1)
    # get the number of neighbors of type j
    number_neighbors = comp[j]

    # if composition has type j...
    if number_neighbors != 0
      # get the sample class
      # return 0 if (l, j) is not a reacting pair
      s = get(reactant_to_class, (l, j), 0)

      # if (l,j) takes part in a reaction...
      if s != 0
        # push the appropriate sample class to the vector
        push!(sample_classes, s + number_neighbors - 1)
      end
    end
  end

  # finally, check to see if (l, 0) is a reactant pair
  # meaning l takes part in an on-site reaction
  s = get(reactant_to_class, (l, 0), 0)

  if s != 0
    # no need to correct s here
    push!(sample_classes, s)
  end

  return sample_classes
end

##### accessors

number_compositions(enum::SampleClassEnumeration) = length(enum.composition)

##### query API

# get the index assigned to a particular neighborhood class
function get_nbclass_index(enum::SampleClassEnumeration, comp)
  k = searchsortedfirst(enum.composition, comp, rev = true)

  if k > length(enum.composition)
    msg = """
    Composition $(comp) not found.
    """
    error(msg)
  end

  return k
end

# derive the neighbor class of a site with
# (i) nbmax - 1 open sites
# (ii) 1 neighbor of type i
function derive_nbclass(enum::SampleClassEnumeration{D,M}, i) where {D,M}
  # pull up the worker array
  comp = enum.dummy_composition

  # build the desired composition
  comp[1] = capacity(M(), D) - 1
  comp[i] = 1

  # search for the index; assumes composition is sorted somehow
  k = get_nbclass_index(enum, comp)

  # cleanup the mess
  fill!(comp, 0)

  return k
end

# get neighbor class after a i -> j type change,
# given that the site was in class k
function get_new_nbclass(enum::SampleClassEnumeration, k, i, j)
  # retrieve the list of compositions
  composition = enum.composition

  # retrieve the worker array
  comp = enum.dummy_composition

  # copy the composition we need to operate on
  copyto!(comp, composition[k])

  # generate the new composition
  comp[i] -= 1
  comp[j] += 1

  k_new = get_nbclass_index(enum, comp)

  # clean up
  fill!(comp, 0)

  return k_new
end

function add_class_member!(class, x::Site)
  i = label(x)
  idx = searchsortedfirst(class, i)
  insert!(class, idx, i)

  return nothing
end

function rmv_class_member!(class, x::Site)
  i = searchsortedfirst(class, label(x))
  deleteat!(class, i)

  return nothing
end

function sample_from_class(lattice, enum, s)
  return get_site(lattice, rand(enum.class[s]))
end

##### TODO Here be dragons.

# Update neighbor and sample classes for x, given that
# 	x changes from type i1 to j1
# 	y changes from type i2 to j2
function update_classes_particle!(enumeration, x, i1, j1, i2, j2)
  # unpack information
  k_old = get_neighbor_class(x)

  # 1. update neighborhood class if necessary
  if i2 != j2
    k_new = get_new_nbclass(enumeration, k_old, i2, j2)
    change_neighbor_class!(x, k_new)
  end

  # 2. update sample class
  # Case 1: i1 == j1 and i2 != j2; neighborhood structure changed, but not type
  if i1 == j1 && i2 != j2
    update_classes_neighbor_change!(enumeration, x, i1, k_old, i2, j2)
  end

  if i1 != j1
    # Case 2: i1 != j1 and i2 == j2; neighborhood structure for x did not change
    if i2 == j2
      # removals
      sample_class_removals!(enumeration, x, i1, k_old)

      # additions
      sample_class_additions!(enumeration, x, j1, k_old)
    else
    # Case 3: i1 != j1 and i2 != j2; neighborhood structure for x did change

      # removals
      sample_class_removals!(enumeration, x, i1, k_old)

      # additions
      sample_class_additions!(enumeration, x, j1, k_new)
    end
  end
  return nothing
end

function sample_class_removals!(enumeration, x, l, k)
  # unpack information
  isactive = enumeration.isactive
  pair_to_classes = enumeration.pair_to_classes
  class = enumeration.class

  if isactive[l]
    affected_classes = pair_to_classes[(l, k)]

    for s in affected_classes
      rmv_class_member!(class[s], x)
    end
  end

  return enumeration
end

function sample_class_additions!(enumeration, x, l, k)
  # unpack information
  isactive = enumeration.isactive
  pair_to_classes = enumeration.pair_to_classes
  class = enumeration.class

  if isactive[l]
    affected_classes = pair_to_classes[(l, k)]

    for s in affected_classes
      add_class_member!(class[s], x)
    end
  end

  return enumeration
end

function update_classes_neighbor_change!(enumeration, x, i1, k_old, i2, j2)
  isactive = enumeration.isactive
  reactant_to_class = enumeration.reactant_to_class
  class = enumeration.class
  composition = enumeration.composition

  if isactive[i1]
    # what class corresponds to the reactant pair (i1, i2)?
    s = get(reactant_to_class, (i1, i2), 0)

    # updates for losing particle of type i2
    if s != 0
      # get number of neighbors of type i2, the interacting particle
      number_neighbors = composition[k_old][i2]

      # index shifting
      s += number_neighbors - 1

      # remove from the sample class
      rmv_class_member!(class[s], x)

      if number_neighbors != 1
        # lost a particle of type i2, so add it to the previous class
        add_class_member!(class[s-1], x)
      end
    end

    # updates for gaining a particle of type j2
    s = get(reactant_to_class, (i1, j2), 0)

    if s != 0
      number_neighbors = composition[k_old][j2]

      s += number_neighbors

      add_class_member!(class[s], x)

      if number_neighbors != 0
        rmv_class_member!(class[s-1], x)
      end
    end
  end

  return nothing
end

# update neighbors due to x interacting with y, where x changes from i to j
function update_classes_neighbors!(lattice, enumeration, x, y, i, j)
  # N = lattice.N
  nbtype = topology(lattice)

  # check the coordinates around x under the given neighborhood type
  for coord in eachdir(nbtype, x)
    # if there is a particle at the coordinate...
    if istracked(lattice, coord)
      neighbor_site = get_site(lattice, coord)

      # if the neighbor site is not the one that interacted with x...
      if neighbor_site != y
        # change neighborhood classes
        k_old = get_neighbor_class(neighbor_site)
        k_new = get_new_nbclass(enumeration, k_old, i, j)

        change_neighbor_class!(neighbor_site, k_new)

        # change sample classes
        l = get_ptype(neighbor_site)
        update_classes_neighbor_change!(enumeration, neighbor_site, l, k_old, i, j)

        # 3. update counts only if neighbor_site did not reaction with x
        # N[k_old, l] -= 1
        # N[k_new, l] += 1
      end
    # otherwise, there is an open site at the coordinate
    else
      neighbor_site = spawn_new_site(lattice, coord)
      l = get_ptype(neighbor_site)

      # set neighborhood class and assign sample classes
      k = derive_nbclass(enumeration, j)
      change_neighbor_class!(neighbor_site, k)
      sample_class_additions!(enumeration, neighbor_site, l, k)

      # add the site
      add_site!(lattice, neighbor_site)

      # update the neighborhood
      for coord2 in eachdir(nbtype, neighbor_site)
        if istracked(lattice, coord2)
          z = get_site(lattice, coord2)

          is_neighbor = label(z) in neighborhood(lattice, neighbor_site)

          if !is_neighbor
            add_neighbor!(lattice, neighbor_site, z)
            add_neighbor!(lattice, z, neighbor_site)
          end
        end
      end
    end
  end

  return nothing
end

# we could make this use a binary search if neighbors
# are stored in a sorted fashion
function sample_neighbor(x, lattice, enum, l)
  # build an iterable object that returns the type of site for each neighbor
  site = lattice.site
  nb   = neighborhood(lattice, x)
  iter = (get_ptype(site[id]) for id in nb)

  # extract the number of type l neighbors
  k = get_neighbor_class(x)
  number_l = enum.composition[k][l]

  # generate a random index into nb
  ix = sample_neighbor(iter, number_l, l)

  if ix == 0
    msg = """
    $(get_ptype(x))|$(get_neighbor_class(x)) @ $(label(x)) / $(x) does not have a neighbor of type $(l).
      Adjacent types: $(map(get_ptype, lattice.site[nb]))
    """
    error(msg)
  end

  # return the chosen Site object
  get_site(lattice, nb[ix])
end

function sample_neighbor(iter, number_l, l)
  c = number_l * rand()
  s = 0 # tracks the number of type l neighbors we've checked
  j = 0 # tracks the random index we want to generate

  for val in iter
    s >= c && break
    isequal(val, l) && (s += 1)
    j += 1
  end

  return j
end

# # TESTING
#
# number_trials = 10_000
# x = rand(0:9, 20)
# l = rand(0:9)
# number_l = count(isequal(l), x)
# iter = (x[j] for j in eachindex(x))
#
# # 1. we always return an index mapping to l
#
# for k in 1:number_trials
#   i = sample_neighbor(iter, number_l, l)
#   @assert x[i] == l
# end
#
# # 2. every index is mapping to l is touched
#
# idxs = Int[]
#
# for k in 1:number_trials
#    push!(idxs, sample_neighbor(iter, number_l, l))
# end
# @assert length(unique(idxs)) == number_l
#
# # 3. the distribution should be uniform
# histogram(idxs)
#


text:
module Assembly # Bukdu

# https://github.com/MikeInnes/Charlotte.jl
# https://github.com/MikeInnes/WebAssembly.jl

end # module Bukdu.Assembly


text:
#############################################################################
# SCS.jl
# Wrapper around the SCS solver https://github.com/cvxgrp/scs
# See http://github.com/karanveerm/SCS.jl
#############################################################################
# test/mpb.jl
# Test the MathProgBase.jl interface for the SCS.jl solver wrapper
#############################################################################

using Base.Test
using MathProgBase
using SCS

solver = SCSSolver()
objtol = 1e-4
primaltol = 1e-4

# min -x
# s.t. 2x + y <= 1.5
# x,y >= 0
# solution is (0.75,0) with objval -0.75
sol = linprog([-1,0],[2 1],'<',1.5,solver)
@test sol.status == :Optimal
@test_approx_eq_eps sol.objval -0.75 objtol
@test_approx_eq_eps norm(sol.sol - [0.75,0.0]) 0 primaltol

sol = linprog([-1,0],sparse([2 1]),'<',1.5,solver)
@test sol.status == :Optimal
@test_approx_eq_eps sol.objval -0.75 objtol
@test_approx_eq_eps norm(sol.sol - [0.75,0.0]) 0 primaltol

# test infeasible problem:
# min x
# s.t. 2x+y <= -1
# x,y >= 0
sol = linprog([1.0,0.0],[2.0 1.0],'<',-1.0,solver)
@test sol.status == :Infeasible

# test unbounded problem:
# min -x-y
# s.t. -x+2y <= 0
# x,y >= 0
sol = linprog([-1.,-1.],[-1. 2.],'<',[0.0],solver)
@test sol.status == :Unbounded

# More complicated through loadproblem!
A = [1.0 1.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 1.0;
     0.0 0.0 1.0 1.0 1.0]
collb = [0.0, 0.0, 0.0, 0.0, 0.0]
colub = [Inf, Inf, Inf, Inf, Inf]
obj   = [3.0, 4.0, 4.0, 9.0, 5.0]
sense = :Max
rowlb = [-Inf, -Inf, -Inf]
rowub = [ 5.0,  3.0,  9.0]
m = MathProgBase.LinearQuadraticModel(solver)
MathProgBase.loadproblem!(m, A, collb, colub, obj, rowlb, rowub, sense)
MathProgBase.optimize!(m)
@test_approx_eq_eps MathProgBase.getobjval(m) 99.0 1e-3
x = MathProgBase.getsolution(m)
@test_approx_eq_eps x[1] 2.0 1e-4
@test_approx_eq_eps x[2] 3.0 1e-4
@test_approx_eq_eps x[3] 0.0 1e-4
@test_approx_eq_eps x[4] 9.0 1e-4
@test_approx_eq_eps x[5] 0.0 1e-4


text:
#####
##### Calculate the areas, volumes and spacings of different kind of `AbstractGrid`
##### Currently for RegularRectilinearGrid and RegularLatLonGrid
#####

#####
##### RegularRectilinearGrid
#####
@inline ΔxC(i, j, k, g::RegularRectilinearGrid) = g.Δx
@inline ΔyC(i, j, k, g::RegularRectilinearGrid) = g.Δy
@inline ΔzC(i, j, k, g::RegularRectilinearGrid) = g.Δz

@inline ΔxF(i, j, k, g::RegularRectilinearGrid) = g.Δx
@inline ΔyF(i, j, k, g::RegularRectilinearGrid) = g.Δy
@inline ΔzF(i, j, k, g::RegularRectilinearGrid) = g.Δz

@inline Ax(i, j, k, g::RegularRectilinearGrid) = g.Δy * g.Δz
@inline Ay(i, j, k, g::RegularRectilinearGrid) = g.Δx * g.Δz
@inline Az(i, j, k, g::RegularRectilinearGrid) = g.Δx * g.Δy

@inline volume(i, j, k, g::RegularRectilinearGrid) = Az(i, j, k, g) * g.Δz

#####
##### RegularLatLonGrid (degree to meter), halo points included
#####
@inline ΔxC(i, j, k, g::RegularLatLonGrid) = @inbounds g.dxC[i,j]
@inline ΔyC(i, j, k, g::RegularLatLonGrid) = @inbounds g.dyC[i,j]
@inline ΔzC(i, j, k, g::RegularLatLonGrid) = g.Δz

@inline ΔxF(i, j, k, g::RegularLatLonGrid) = @inbounds g.dxF[i,j]
@inline ΔyF(i, j, k, g::RegularLatLonGrid) = @inbounds g.dyF[i,j]
@inline ΔzF(i, j, k, g::RegularLatLonGrid) = ΔzC(i, j, k, g)

@inline Ax(i, j, k, g::RegularLatLonGrid) = @inbounds g.Ax[i,j]
@inline Ay(i, j, k, g::RegularLatLonGrid) = @inbounds g.Ay[i,j]
@inline Az(i, j, k, g::RegularLatLonGrid) = @inbounds g.Az[i,j]

@inline volume(i, j, k, g::RegularLatLonGrid) = @inbounds g.Vol[i,j]

#####
##### VerticallyStretchedLatLonGrid (degree to meter), halo points included
#####
@inline ΔxC(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.dxC[i,j]
@inline ΔyC(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.dyC[i,j]
@inline ΔzC(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.dzC[k]

@inline ΔxF(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.dxF[i,j]
@inline ΔyF(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.dyF[i,j]
@inline ΔzF(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.dzF[k]

@inline Ax(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.Ax[i,j,k]
@inline Ay(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.Ay[i,j,k]
@inline Az(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.Az[i,j]

@inline volume(i, j, k, g::VerticallyStretchedLatLonGrid) = @inbounds g.Vol[i,j,k]


text:
using Test, RecursiveArrayTools, StaticArrays

using DiffEqBase: UNITLESS_ABS2, recursive_length,  ODE_DEFAULT_NORM

n = UNITLESS_ABS2(3.0+4.0im)
@test n == 25.0
@test typeof(n)<:Real

@test ODE_DEFAULT_NORM(3.0+4.0im, 0.0) == 5.0

u1 = ones(3)
@test UNITLESS_ABS2(u1) == 3.0
@test recursive_length(u1) == 3
@test ODE_DEFAULT_NORM(u1, 0.0) == 1.0

u2 = [SA[1.0 1.0; 1.0 1.0] for i = 1:3]
@test UNITLESS_ABS2(u2) == 12.0
@test recursive_length(u2) == 12
@test ODE_DEFAULT_NORM(u2, 0.0) == 1.0

u3 = VectorOfArray([ones(5), ones(5)])
@test UNITLESS_ABS2(u3) == 10.0
@test recursive_length(u3) == 10
@test ODE_DEFAULT_NORM(u3, 0.0) == 1.0

u4 = ArrayPartition(u1, u2, u3)
@test UNITLESS_ABS2(u4) == 25.0
@test recursive_length(u4) == 25
@test ODE_DEFAULT_NORM(u4, 0.0) == 1.0

u5 = ArrayPartition(u4, u4)
@test UNITLESS_ABS2(u5) == 50.0
@test recursive_length(u5) == 50
@test ODE_DEFAULT_NORM(u5, 0.0) == 1.0


text:
using CUDAdrv
using CUDAnative

mlock!(p::Ptr, size::Integer) = ccall(:mlock, Cint, (Ptr{Cvoid}, Csize_t), p, size)
munlock!(p::Ptr, size::Integer) = ccall(:unmlock, Cint, (Ptr{Cvoid}, Csize_t), p, size)

function pin!(a)
    ad = Mem.register(Mem.Host, pointer(a), sizeof(a))
    finalizer(_ -> Mem.unregister(ad), a)
end

"""
    empiricalbandwidth(nbytes=2*1024^3; devicenumber=0, numtests=10)

Compute the emperical bandwidth in GB/s of a CUDA device using `nbytes` of
memory.  The device to test can be slected with `devicenumber` and the
bandwidth is an average of `ntests`.
"""
function empiricalbandwidth(nbytes=1024^2; devicenumber=0, ntests=10, pinned=false)
    device!(devicenumber)
    stm = CuStream(CUDAdrv.STREAM_NON_BLOCKING)

    a = rand(Char, nbytes÷sizeof(Char))
    pinned && pin!(a)

    b = Mem.alloc(Mem.Device, nbytes)

    Mem.copy!(pointer(a), b, nbytes)
    Mem.copy!(b, pointer(a), nbytes)

    t_dtoh = CUDAdrv.@elapsed stm for n = 1:ntests
        Mem.copy!(pointer(a), b, nbytes, async=true, stream=stm)
    end
    bandwidth_dtoh = nbytes*ntests/(t_dtoh*1e9)

    t_htod = CUDAdrv.@elapsed stm for n = 1:ntests
        Mem.copy!(b, pointer(a), nbytes, async=true, stream=stm)
    end
    bandwidth_htod = nbytes*ntests/(t_htod*1e9)

    Mem.free(b)

    (bandwidth_dtoh, bandwidth_htod)
end


text:
#!/usr/bin/env julia
include("src/NamedArrays.jl")
using NamedArrays

cd("test")
include("test/runtests.jl")


text:
######################IO################

vec2String(args...; keyargs...) = error("Depreciated- replace with join")
texTable(args...; keyargs...) = error("Depreciated- replace with textable")
writeTexTable(args...; keyargs...) = error("Depreciated- replace with writetextable")
array2String(args...; keyargs...) = error("Depreciated- replace with array2string")

#functions for converting Floats to Strings
function num2str(x::T, decimals::Int=DEFAULT_DECIMALS;
    scalefactor::Float64 = 1.0,
    scalefactoronhurdle::Float64 = 1.0,
    Ints::Bool = false, scalehurdle::U=nothing, texequation::Bool = false)::String where {
      T<:Union{Real,Missing,Nothing}, U<:Union{Nothing,Real}}

  local outstr::String

  x*=scalefactor

  #this controls for both missing and nothing values
  if ismissing(something(x,missing))
    outstr = ""
  elseif scalehurdle ≠ nothing && x≥scalehurdle
    outstr = "$(Int(round(x*scalefactoronhurdle)))"
  elseif decimals == 0 || (Ints && round(x) == x)
    outstr = "$(Int(round(x)))"
  else
    outstr = format("{:.$(decimals)f}",x)
  end

  outstr = texequation ? "\$$outstr\$" : outstr

  return outstr
end

#depreciated functions
#array2String
  #=takes content of an array, formats each cell into a string, appends a
  prefix and suffix, and sends back an array of strings with the same
  dimension as the input
  IN: A general array, optionally a prefix, suffix, scaling factor, and
   set number of decimal digits
  OUT: An array of strings with the same dimensions as the input=#
function array2string(m::Array{Float64}; prefix::String = "",
  suffix::String = "", scaling::Float64 = 1.0,
  decimaldigits::Int = 2)::Array{String}

  #broadcast the prefix, value in m, and suffix and concatenate respectively
  stringmat::Array{String} =
      ((x::Float64)->(prefix*"$(round(x,decimaldigits))"*suffix)).(m)

  return stringmat
end

#textable
#=    Creates a well-formed Latex table from string components
  The content matrix consists of an array of matrices, which are overlayed
  via a one line offset. For example, the first matrix might be the values,
  the second matrix the errors. Observational rows can be put in
  desccontent, which is a vector of rows. Similarly, colnames is a vector
  of rows. So that column headers can span multiple columns, the widthcolnames
  and widthdesccontent provides the option to specify dimentions
  IN: colnames, rownames, content, desccontent (rows of
      descriptive data), notes, optionally widthcolnames and widthDescConent
      (which are the width of each entry in columns)
  OUT: The tex table as a string=#

function textable(;
  colnames::Vector{Vector{String}} = Vector{Vector{String}}(),
  contentrownames::Vector{String} = Vector{String}(),
  content::Vector{Matrix{String}} = Vector{Matrix{String}}(),
  descrownames::Vector{String} = Vector{String}(),
  desccontent::Vector{Vector{String}}=Vector{Vector{String}}(),
  notes::Vector{String} = Vector{String}(),
  linespacer::String = "\\\\",
  summarymathmode::Bool = true,
  widthcolnames::Vector{Vector{Int}} = #contains the number of columns for each entry
    broadcast((i::Int)->ones(Int,length(colnames[i])),1:length(colnames)),
  alignmentcolnames::Vector{Vector{String}} = #contains the number of columns for each entry
    broadcast((i::Int)->["r" for i ∈ 1:length(colnames[i])],1:length(colnames)),
  widthdesccontent::Vector{Vector{Int}} = #contains the number of columns for each entry
    broadcast((i::Int)->ones(Int,length(desccontent[i])),1:length(desccontent)),
  colheadername::Vector{String} = ["" for i ∈ 1:length(colnames)],
  alignmentstring::String = string(" l", join(["r" for i ∈ 1:length(desccontent[1])])),
  rowlabelheader::Bool = false)

  #size parameters for content
  numcontentrows::Int = length(contentrownames)

  numcontentcols::Int = sum(widthcolnames[1]) - rowlabelheader

  numcontentsubrows::Int = length(content)

      #intiate the stream
  b::IOBuffer = IOBuffer()
  write(b, """
      \\begin{tabular}{$alignmentstring""")

    write(b,"}\n \\toprule")#filler tex

  for r::Int ∈ 1:length(colnames) #for each row of column headings
      #write(b,"\n\\\\[-1.8ex]") #use this if the below doesn't work
    write(b, colheadername[r])
    if length(colnames[r]) ≠ numcontentcols
      for c::Int ∈ 1:length(colnames[r]) #for each column heading
        (c≠1 || (!rowlabelheader)) && write(b,"\t&")
        write(b,"\t\\multicolumn{$(widthcolnames[r][c])}{$(alignmentcolnames[r][c])}{$(colnames[r][c])}")
      end
    else
      for c::Int ∈ 1:length(colnames[r]) #for each column heading
        (c≠1 || (!rowlabelheader)) && write(b,"\t&")
        write(b,"\t$(colnames[r][c])")
      end
    end
      write(b,"\n\\\\")
  end
  write(b, " \\midrule\n ")

  #now write out the table content
    if numcontentrows > 0
        write(b," \t \t ")
        for r::Int ∈ 1:numcontentrows
          write(b, "$(contentrownames[r])") #the row label
          for s::Int ∈ 1:numcontentsubrows #print the sub-rows for each row
            for c::Int ∈ 1:numcontentcols #print the columns for each sub-row
              write(b,"\t&\t$((content[s])[r, c])")
            end
            write(b,"\n $(linespacer) \t\t") #line-break stylistic formatting
          end
        end
        write(b, "\n \\midrule\n")
    end
    sumMathFlag::String = summarymathmode ? "\$" : ""

    if length(desccontent)>0
        for r::Int ∈ 1:length(desccontent) #for each description row
          write(b, "$(descrownames[r])") #the row label
          if length(widthdesccontent[r]) ≠ numcontentcols
            for c::Int ∈ 1:length(desccontent[r]) #for each descriptive row
              if length(desccontent[r][c]) ≥ 1
                write(b,"\t&\t\\multicolumn{$(widthdesccontent[r][c])}{r}{$(sumMathFlag)$(desccontent[r][c])$(sumMathFlag)}")
              else
                write(b,"\t&\t\\multicolumn{$(widthdesccontent[r][c])}{r}{}")
              end
            end
          else
            for c::Int ∈ 1:length(desccontent[r]) #for each descriptive row
              if length(desccontent[r][c]) ≥ 1
                write(b,"\t&\t$(sumMathFlag)$(desccontent[r][c])$(sumMathFlag)")
              else
                write(b,"\t&\t")
              end
            end
          end
            write(b,"\n $(linespacer) \t\t") #line-break stylistic formatting
        end
    end

      write(b,""" \\bottomrule""")
  #if we have footnotes
  if length(notes) > 0
      write(b,"""\n \\\\[-1.0ex] \\textit{Notes:} \t """)
      for r ∈ length(notes)
          write(b, "\t&\t \\multicolumn{$numcontentcols}{l}{$(notes[r])}\n \\\\")
      end
  end
  write(b, "\t \\end{tabular}\n")
  return String(take!(copy(b)))
end

#=this version formats a table from a set of linear models
While customizable, the idea here is to have sensible default arguments
where possible.
IN: A collection of models, the rows to capture (returns an empty string in
the table if a row is not present, optionally column names, their width,
names of the rows, names of the descriptive rows, descriptive content,
a boolean switch for including signficance stars, signficance levels for the stars,
a scaling factor, the number of digits in values, and customizable
notes
OUT: A latex table string
=#
function textable(models::Vector{FMLM},
    getΣ::T where T<:Union{Function, Vector{Function}},
    rows::Vector{String};
    colnames::Vector{Vector{String}} = [["($i)" for i ∈ 1:length(models)]],
    contentrownames::Vector{String} = String.(rows),
    descrownames::Vector{String}=Vector{String}(),
    desccontent::Vector{Vector{String}}=Vector{Vector{String}}(),
    notes::Vector{String} = Vector{String}(),
    linespacer::String = "\\\\",
    summarymathmode::Bool = true,
    widthcolnames::Vector{Vector{Int}} =
        broadcast((i::Int)->ones(Int,length(colnames[i])),1:length(colnames)),
    alignmentcolnames::Vector{Vector{String}} = #contains the number of columns for each entry
      broadcast((i::Int)->["r" for i ∈ 1:length(colnames[i])],1:length(colnames)),
    widthdesccontent::Vector{Vector{Int}} =
        broadcast((i::Int)->ones(Int,length(desccontent[i])),1:length(desccontent)),
    stars::Bool=true,
    starlvls::Vector{Float64} = [.9, .95, .99],
    starlegend::String = stars ? DEFAULT_STAR_LEGEND : "",
    starstrings::Vector{String} =
      ["\\ensuremath{^\\text{*}}","\\ensuremath{^\\text{**}}","\\ensuremath{^\\text{***}}"],
    scaling::Vector{Float64}=ones(length(rows)),
    decimaldigits::Int = 2,
    colheadername::Vector{String} = ["" for i::Int ∈ 1:length(colnames)],
    alignmentstring::String = string(" l | ", join(["r" for i ∈ 1:length(desccontent[1])])),
    rowlabelheader::Bool = false)

  Ncols = length(models)
  numcontentrows = length(rows)

  #Pre-allocate the vector of SE errors
  σs::Vector{Vector{Float64}} =
      [Vector{Float64}(undef, models[i].K) for i∈1:Ncols]

  #pull out the β coefficients and N
  βs::Vector{Vector{Float64}} = [models[i].β for i::Int ∈ 1:Ncols]
  Ns::Vector{Int} = [models[i].N for i ∈ 1:Ncols]
  Xnames::Vector{Vector{String}} = [models[i].Xnames for i::Int ∈ 1:Ncols]

  #get the standard errors
  for c ∈ 1:Ncols
      if typeof(getΣ) <: Function
        σs[c] .= sqrt.(diag(getΣ(models[c])))
      else
        σs[c] .= sqrt.(diag(getΣ[c](models[c])))
      end
  end

  if length(starlegend) > 0
      notes = [starlegend; notes]
  end

  #get the content matrices
  content::Vector{Matrix{String}} = getcontentmatrices!( βs=βs, σs=σs,
    Xnames=Xnames, Ns=Ns, rows=rows,
    stars=stars, starlvls=starlvls, scaling=scaling,
    decimaldigits=decimaldigits, starstrings=starstrings)

  return textable(colnames=colnames,
    contentrownames=contentrownames, content=content, descrownames=descrownames,
    desccontent=desccontent, notes=notes, linespacer=linespacer,
    summarymathmode=summarymathmode, widthcolnames=widthcolnames,
    widthdesccontent=widthdesccontent, colheadername=colheadername,
    alignmentcolnames=alignmentcolnames,
    alignmentstring=alignmentstring, rowlabelheader=rowlabelheader)
end

#=getcontentmatrices!
this is a bit of utiltiy code for making tex tables. It is not model specific
hence why it was extracted into its own function
#IN: a vector of βs, σs, names, rows (coefficeints selected). Optional parameters
include a pre-allcoation of the string matrix, a switch for the inclusion of
stars, a scaling factor, and the number of digits (rounding level)=#
#OUT: Writes to and returns the content matrix
function getcontentmatrices!(;
    βs::Vector{Vector{Float64}} = error("βs is a required argument for getcontentmatrices"),
    σs::Vector{Vector{Float64}} = error("σs is a required argument for getcontentmatrices"),
    Xnames::Vector{Vector{String}} = error("Xnames is a required argument for getcontentmatrices"),
    Ns::Vector{Int} = error("Ns is a required argument for getcontentmatrices"),
    rows::Vector{String} = error("rows is a required argument for getcontentmatrices"),
    stars::Bool=true, #whether to display signficance stars
    starlvls::Vector{Float64} = [.9, .95, .99],  #cutoffs for signficance (must be sorted)
    starstrings::Vector{String} =
      ["\\ensuremath{^\\text{*}}","\\ensuremath{^\\text{**}}","\\ensuremath{^\\text{***}}"],
    scaling::Vector{Float64}=ones(length(rows)), # an optional scaling factor
    decimaldigits::Int = 2) #number of decimal digits

    content::Vector{Matrix{String}} = #will hold the coefficients and errors
        [fill("",length(rows),length(βs)) for i ∈ 1:2];

    #iterate through the models
    for c ∈ 1:length(βs)

        #build a dictionary of the names
        XNameTbl::Dict = Dict(Xnames[c][i] => i for i::Int ∈ 1:length(βs[c]))

        for r ∈ 1:length(rows)
            if haskey(XNameTbl, rows[r]) #need to check if it exists
                ind::Int = XNameTbl[rows[r]]
                p::Float64 =
                    cdf(TDist(Ns[c]), βs[c][ind]/σs[c][ind]) #get CDF from T distribution
                p = p > .5 ? 1-(1.0 - p)*2.0 : 1.0 - p*2.0 #calc 2-tailed p value
                sigLevel::Int = sum(p.>starlvls)
                if sigLevel > 0 && stars
                    starString::String = starstrings[sum(p.>starlvls)]
                else
                    starString = ""
                end

                #scale, round and write the β coefficeint and σ into the string matrices
                content[1][r,c] =
                  #"$(round(scaling[r]*βs[c][ind],decimaldigits))^{$starString}"
                  "\$$(num2str(βs[c][ind], decimaldigits, scalefactor=scaling[r]))\$$starString"
                content[2][r,c] =
                  "(\$$(num2str(σs[c][ind], decimaldigits, scalefactor=scaling[r]))\$)"
                  #"($(round(scaling[r]*σs[c][ind],decimaldigits)))"
            end
        end
    end

  return content
end

#=Writes the fully formed string to the file
  IN: A string with the output, optionally a path to the file, and the
  output name
  OUT: Writes the string to a file=#
function writetextable(tablestring::String;
  path::String=pwd(),
  outname::String = "$path\\results$(Dates.format(now(),"yymmdd_THMS")).tex")

  oStream::IOStream = open("$path\\$outname","w+")
  write(oStream, tablestring)
  close(oStream)
end

function iotest()
  #test parameters
  nrowsContent = 8
  ncols = 4
  nsecs = ncols ÷ 2
  nrowsDesc = 2
  nsubRows = 2

  path = pwd() * "\\results"
  outname = "test.tex"
  footerName = "footer.tex"
  headerName = "header.tex"

  #test table
  colnames::Vector{Vector{String}} =
      [broadcast(i->"secs$i",1:nsecs),broadcast(i->"ncols$i",1:ncols)]
  contentrownames::Vector{String} = ["nrows$i" for i ∈ 1:nrowsContent]
  content::Vector{Matrix{String}} =
      (string).(broadcast((i::Int)->i .* ones(nrowsContent,ncols), 1:nsubRows))
  descrownames::Vector{String} = ["Desc. row $i" for i ∈ 1:nrowsDesc]
  desccontent::Vector{Vector{String}} =
      [broadcast(i->"desc-r1c$i",1:ncols),broadcast(i->"desc-r1c$i",1:ncols)]
  notes::Vector{String} = ["a note"]
  widthcolnames::Vector{Vector{Int}} =
      [broadcast(i->(ncols ÷ nsecs), 1:nsecs),broadcast(i->1, 1:ncols)]

  s=textable(colnames,  contentrownames, content, descrownames,
      desccontent, notes, widthcolnames=widthcolnames)

  println(array2string(fill(.45/π,3,3,4),decimaldigits=4))

  print(s)
end

max_stars_repo_name: clintonTE/Finometrics

max_issues_repo_name: clintonTE/Finometrics

max_forks_repo_name: clintonTE/Finometrics


text:
using JuliaParser
using FactCheck

const Parser = JuliaParser.Parser
const Lexer  = JuliaParser.Lexer

include("ast.jl")
include("util.jl")

# TODO: Julia does not really have a good way to do this so
# this only currently works from source builds.
if length(ARGS) > 0
    const BASEPATH = ARGS[1]
else
    const BASEPATH = abspath(joinpath(JULIA_HOME, "..", ".."))
end

passed = 0
failed = 0
errors = 0

ptime = 0.0
btime = 0.0

function testall(srcdir::AbstractString)
    global passed
    global failed
    global errors
    global ptime
    global btime

    dirs, files = [], []

    for fname in sort(readdir(srcdir))
        path = joinpath(srcdir, fname)
        if isdir(path)
            push!(dirs, path)
            continue
        end
        _, ext = splitext(fname)
        if ext == ".jl"
            push!(files, path)
        end
    end

    if !isempty(files)
        println(bold("test $srcdir"))

        for jlpath in files

            fname = splitdir(jlpath)[end]

            buf = IOBuffer()
            write(buf, "begin\n")
            write(buf, open(readstring, jlpath))
            write(buf, "\nend")

            src = String(buf)
            try
                tic()
                past = Parser.parse(src)
                t1 = toq()

                tic()
                bast = Base.parse(src)
                t2 = toq()

                if norm_ast(past) == norm_ast(bast)
                    println(green("OK:     $fname"))
                    passed += 1
                    ptime += t1
                    btime += t2
                else
                    println(red("FAILED: $fname"))
                    failed += 1
                end
            catch ex
                println(bold(red("ERROR:  $fname")))
                errors += 1
            end
        end
        println()
    end
    for dir in dirs
        testall(dir)
    end
end

if isdir(BASEPATH) && isdir(joinpath(BASEPATH, "base"))
    testall(joinpath(BASEPATH, "examples"))
    testall(joinpath(BASEPATH, "test"))
    testall(joinpath(BASEPATH, "base"))
else
    warn("""
Could not find julia base sources in $BASEPATH,
perhaps you are using a Julia not built from source?""")
end

println("\n"^2)
println(bold("TOTAL: $(passed + failed + errors)"))
println(green("Passed:\t$passed"))
println(red("Failed:\t$failed"))
println(red("Errors:\t$errors"))
println()

pct = ptime / btime * 100.0

if pct >= 100
    pctstr = @sprintf("%0.2f", pct-100)
    println(bold("Parser is ~$pctstr% slower than base"))
else
    pctstr = @sprintf("%0.2f", 100-pct)
    println(bold("Parser is ~$pctstr% faster than base"))
end

#exit(failed + errors)


text:
using Catlab.Graphics.Graphviz
import Catlab.Graphics.Graphviz: Graph, Edge
import Base.Iterators: flatten
using StatsBase

export Graph

function edgify(δ, transition, reverse::Bool)
  return [Edge(reverse ? ["T_$transition", "S_$k"] :
                         ["S_$k", "T_$transition"],
               Attributes(:label=>"$(δ[k])", :labelfontsize=>"6"))
           for k in collect(keys(δ)) if δ[k] != 0]
end

function Graph(p::AbstractPetriNet)
  statenodes = [Node(string("S_$(sname(p, s))"), Attributes(:shape=>"circle", :color=>"#6C9AC3")) for s in 1:ns(p)]
  transnodes = [Node(string("T_$(tname(p, k))"), Attributes(:shape=>"square", :color=>"#E28F41")) for k in 1:nt(p)]

  graph_attrs = Attributes(:rankdir=>"LR")
  node_attrs  = Attributes(:shape=>"plain", :style=>"filled", :color=>"white")
  edge_attrs  = Attributes(:splines=>"splines")

  stmts = vcat(statenodes, transnodes)

  edges = map(1:nt(p)) do k
    k_name = tname(p, k)
    vcat(edgify(countmap(map(x->sname(p,x), inputs(p, k))), k_name, false),
         edgify(countmap(map(x->sname(p, x), outputs(p, k))), k_name, true))
  end |> flatten |> collect

  stmts = vcat(stmts, edges)
  g = Graphviz.Digraph("G", stmts; graph_attrs=graph_attrs, node_attrs=node_attrs, edge_attrs=edge_attrs)
  return g
end

function Graph(op::Union{OpenPetriNet, OpenLabelledPetriNetUntyped, OpenReactionNet, OpenLabelledReactionNetUntyped})
    Graph(apex(op))
end


text:
#=
This program benchmarks the performance of a ternary conditional expression
in place of using an if-else block. This is done by iterating over 10^7 random
floats that are normally distributed between [0, 1). A retrun value is incrumented
or decrumented at each iteration, using a conditional block or a ternary statement.

This benchmark is performed in C and Julia, where the results are reported to the
console. There appears to be no major performance differences on my machine.

author: Daniel Henderson
=#
using BenchmarkTools
using Libdl

V = rand(10^7)

# ---------------------- Julia declaration
function ternary(V)
    ret = 0
    for v in V
        v < 0.5 ? (ret += 1) : (ret -= 1)
    end
    ret
end

function conditional(V)
    ret = 0
    for v in V
        if v < 0.5
            ret += 1
        else
            ret -= 1
        end
    end
    ret
end

# ---------------------- C declaration
C_code_cond = """
#include <stddef.h>
int c_conditional(size_t n, double *V) {
    int ret = 0;
    for (size_t i = 0; i < n; ++i) {
        if (V[i] < 0.5) {
            ret++;
        } else {
            ret--;
        }
    }
    return ret;
}
"""

C_code_ternary = """
#include <stddef.h>
int c_ternary(size_t n, double *V) {
    int ret = 0;
    for (size_t i = 0; i < n; ++i) {
        (V[i] < 0.5) ? ret++ : ret--;
    }
    return ret;
}
"""

# compile C code to a temp-shared directory w/ gcc and no optimization
const Clib_c = tempname()
open(`gcc -fPIC -msse3 -xc -shared -o $(Clib_c * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code_cond)
end
const Clib_t = tempname()
open(`gcc -fPIC -msse3 -xc -shared -o $(Clib_t * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code_ternary)
end

# create a Julia aliases to call the C functions
c_conditional(V::Array{Float64}) = ccall(("c_conditional", Clib_c), Int64, (Csize_t, Ptr{Float64}), length(V), V)
c_ternary(V::Array{Float64}) = ccall(("c_ternary", Clib_t), Int64, (Csize_t, Ptr{Float64}), length(V), V)

# ---------------------- Benchmark and Report
println("\nIs ternary() similar to conditional()? ", ternary(V) ≈ conditional(V))
t_bench_Julia = @benchmark ternary(V)
c_bench_Julia = @benchmark conditional(V)

println("\nRESULT OF ternary() -------------------------------------------------")
display(t_bench_Julia)

println("\n\nRESULT OF conditional() ---------------------------------------------")
display(c_bench_Julia)

println("\n\n\nIs c_ternary() similar to c_conditional()? ", c_ternary(V) ≈ c_conditional(V))
t_bench_C = @benchmark c_ternary(V)
c_bench_C = @benchmark c_conditional(V)

println("\nRESULT OF c_ternary() -------------------------------------------------")
display(t_bench_C)

println("\n\nRESULT OF c_conditional() ---------------------------------------------")
display(c_bench_C)
println()

max_stars_repo_name: danphenderson/benchmarks

max_issues_repo_name: danphenderson/benchmarks

max_forks_repo_name: danphenderson/benchmarks


text:
using Printf

function maths_examples()
    println("--------------------------------------")
    println("Matha examples")

    println("Lots of obvious maths functions")
    @printf("5 + 4 = %d\n", (5 + 4))
    @printf("5 - 4 = %d\n", (5 - 4))
    @printf("5 * 4 = %d\n", (5 * 4))
    @printf("5 / 4 = %d\n", (5 / 4))
    @printf("5 %% 4 = %d\n", (5 % 4)) # Note we had to escape the '%' symbol here because it's usually used for identifying parameters!
    @printf("5 ^ 4 = %d\n", (5 ^ 4))
    @printf("round(3.5) = %d\n", round(3.5))
    @printf("floor(3.5) = %d\n", floor(3.5))
    @printf("ceil(3.5) = %d\n", round(3.5))
    @printf("abs(-3.5) = %d\n", abs(-3.5))
    @printf("sqrt(100) = %d\n", sqrt(100))
    @printf("cbrt(100) = %d\n", cbrt(100))
    @printf("hypot(90,90) = %d\n", hypot(90,90))
    @printf("exp(2) = %d\n", exp(2))
    @printf("log(100) = %d\n", log(100))
    @printf("log2(100) = %d\n", log2(100))
    @printf("log10(100) = %d\n", log10(100))

    # Trig Functions : sin, cos, tan, cot, sec, csc
    # sinh, cosh, tanh, coth, sech, csch
    # asin, acos, atan, acot, asec, acsc
    # asinh, acosh, atanh, acoth, asech, acsch
    # sinc, cosc
    # Shortcuts available += -= *= /= ^=

    println("We can multiply explicitly ('2x' instead of '2 * x')")
    x = 5
    println(2x)

    println("The 'Dot' operators perform operations on arrays, e.g. '.* 3' multiply each element in array by 3")
    a13 = [1,2,3]
    println(a13 .* 3)

end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump