text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

using ValidatedNumerics
using FactCheck

facts("Constructing intervals") do
    setprecision(Interval, 53)
    @fact ValidatedNumerics.parameters.precision --> 53

    setprecision(Interval, Float64)
    @fact ValidatedNumerics.parameters.precision --> 53

    # There is an inexplicable error on 0.5 with the following:
    @pending precision(BigFloat) --> 53
    @pending precision(Interval) --> (Float64, 53)

    # Checks for parameters
    @fact ValidatedNumerics.parameters.precision_type --> Float64
    @fact ValidatedNumerics.parameters.precision --> 53
    @fact ValidatedNumerics.parameters.rounding --> :narrow
    @fact ValidatedNumerics.parameters.pi --> @biginterval(pi)

    # Naive constructors, with no conversion involved
    @fact Interval(1) --> Interval(1.0, 1.0)
    @fact Interval(big(1)) --> Interval(1.0, 1.0)
    @fact Interval(eu) --> Interval(1.0*eu)
    @fact Interval(1//10) --> Interval{Rational{Int}}(1//10, 1//10)
    @fact Interval(BigInt(1)//10) --> Interval{Rational{BigInt}}(1//10, 1//10)
    @fact Interval( (1.0, 2.0) ) --> Interval(1.0, 2.0)

    @fact Interval{Rational{Int}}(1) --> Interval(1//1)
    #@fact Interval{Rational{Int}}(pi) --> Interval(rationalize(1.0*pi))

    @fact Interval{BigFloat}(1) --> Interval{BigFloat}(big(1.0),big(1.0))
    @fact Interval{BigFloat}(pi) -->
        Interval{BigFloat}(big(3.1415926535897931), big(3.1415926535897936))

    # Disallowed conversions with a > b

    @fact_throws ArgumentError Interval(2, 1)
    @fact_throws ArgumentError Interval(big(2), big(1))
    @fact_throws ArgumentError Interval(BigInt(1), 1//10)
    @fact_throws ArgumentError Interval(1, 0.1)
    @fact_throws ArgumentError Interval(big(1), big(0.1))

    @fact_throws ArgumentError @interval(2, 1)
    @fact_throws ArgumentError @interval(big(2), big(1))
    @fact_throws ArgumentError @interval(big(1), 1//10)
    @fact_throws ArgumentError @interval(1, 0.1)
    @fact_throws ArgumentError @interval(big(1), big(0.1))

    # Conversions; may involve rounding
    # @fact convert(Interval, 1) --> Interval(1.0)
    # @fact convert(Interval, pi) --> @interval(pi)
    # @fact convert(Interval, eu) --> @interval(eu)
    # @fact convert(Interval, BigInt(1)) --> Interval(BigInt(1))
    # @fact convert(Interval, 1//10) --> @interval(1//10)
    # @fact convert(Interval, 0.1) --> Interval(0.09999999999999999, 0.1)
    # @fact convert(Interval, BigFloat(0.1)) --> Interval(big(0.1))

    @fact convert(Interval{Rational{Int}}, 0.1) --> Interval(1//10)
    # @fact convert(Interval{Rational{BigInt}}, pi) --> Interval{Rational{BigInt}}(pi)

    ## promotion
    @fact promote(Interval(2//1,3//1), Interval(1, 2)) -->
        (Interval(2.0,3.0), Interval(1.0,2.0))
    @fact promote(Interval(1.5), parse(BigFloat, "2.1")) -->
        (Interval(BigFloat(1.5)), Interval(BigFloat(2.1)))
    @fact promote(Interval(1.0), pi) --> (Interval(1.0), @interval(pi))

    # Constructors from the macros @interval, @floatinterval @biginterval
    setprecision(Interval, 53)

    a = @interval(0.1)
    b = @interval(pi)

    @fact nextfloat(a.lo) --> a.hi
    @fact typeof(a) --> Interval{BigFloat}
    @fact a --> @biginterval("0.1")
    @fact convert(Interval{Float64}, a) --> @floatinterval(0.1)
    @fact nextfloat(b.lo) --> b.hi

    @fact b --> @biginterval(pi)
    x = 10238971209348170283710298347019823749182374098172309487120398471029837409182374098127304987123049817032984712039487
    @fact @interval(x) --> @biginterval(x)
    @fact isthin(@interval(x)) --> false

    x = 0.1
    a = @interval(x)
    @fact nextfloat(a.lo) --> a.hi

    setprecision(Interval, Float64)
    a = @interval(0.1)
    b = @interval(pi)

    @fact a --> @floatinterval("0.1")
    @fact typeof(a) --> Interval{Float64}
    @fact nextfloat(a.lo) --> a.hi
    @fact b --> @floatinterval(pi)
    @fact nextfloat(b.lo) --> b.hi
    @fact convert(Interval{Float64}, @biginterval(0.1)) --> a
    x = typemax(Int)
    @fact @interval(x) --> @floatinterval(x)
    @fact isthin(@interval(x)) --> false
    x = rand()
    c = @interval(x)
    @fact nextfloat(c.lo) --> c.hi

    a = @interval("[0.1, 0.2]")
    b = @interval(0.1, 0.2)

    @fact a ⊆ b --> true

    @fact_throws ArgumentError @interval("[0.1]")
    @fact_throws ArgumentError @interval("[0.1, 0.2")

    for precision in (64, Float64)
        setprecision(Interval, precision)
        d = big(3)
        f = @interval(d, 2d)
        @fact @interval(3, 6) ⊆ f --> true
    end

    for rounding in (:wide, :narrow)
        a = @interval(0.1, 0.2)
        @fact a ⊆ Interval(0.09999999999999999, 0.20000000000000004) --> true

        b = @interval(0.1)
        @fact b ⊆ Interval(0.09999999999999999, 0.10000000000000002) --> true

        b = setprecision(Interval, 128) do
            @interval(0.1, 0.2)
        end
        @fact b ⊆ Interval(0.09999999999999999, 0.20000000000000004) --> true

        @fact float(b) ⊆ a --> true

        c = @interval("0.1", "0.2")
        @fact c ⊆ a --> true  # c is narrower than a
        @fact Interval(1//2) == Interval(0.5) --> true
        @fact Interval(1//10).lo == rationalize(0.1) --> true
    end

    @fact string(emptyinterval()) == "∅" --> true

    params = ValidatedNumerics.IntervalParameters()
    @fact params.precision_type == BigFloat --> true
    @fact params.precision == 256 --> true
    @fact params.rounding == :narrow --> true

    setprecision(Interval, 53)
    a = big(1)//3
    @fact @interval(a) --> Interval(big(3.3333333333333331e-01), big(3.3333333333333337e-01))

end

facts("Big intervals") do
    a = @floatinterval(3)
    @fact typeof(big(a)) --> Interval{BigFloat}

    @fact @floatinterval(123412341234123412341241234) --> Interval(1.234123412341234e26, 1.2341234123412342e26)
    @fact @interval(big"3") --> @floatinterval(3)
    @fact @floatinterval(big"1e10000") --> Interval(1.7976931348623157e308, ∞)

    a = big(10)^10000
    @fact @floatinterval(a) --> Interval(1.7976931348623157e308, ∞)
    setprecision(Interval, 53)
    @fact @biginterval(a) --> Interval(big"9.9999999999999994e+9999", big"1.0000000000000001e+10000")

end

facts("Complex intervals") do
    a = @floatinterval(3 + 4im)
    @fact a --> Interval(3) + im*Interval(4)

    b = exp(a)
    @fact real(b) --> Interval(-13.12878308146216, -13.128783081462153)
    @fact imag(b) --> Interval(-15.200784463067956, -15.20078446306795)
end

facts("± tests") do
    setprecision(Interval, Float64)

    @fact 3 ± 0.5 --> Interval(2.5, 3.5)
    @fact 3 ± 0.1 --> Interval(2.9, 3.1)
    @fact 0.5 ± 1 --> Interval(-0.5, 1.5)
end


text:
mutable struct WorkStealingScheduler <: AbstractScheduler
    @const unfairspins::Int
    @const yieldspins::Int
    @const workerindices::Vector{Int}
    @const threadids::Vector{Int}
    @const workers::Vector{Worker}
    @const queues::Vector{WorkStealingDeque{Task}}
    @const externalqueue::ConcurrentQueue{Task}
    @atomic closed::Bool
end

# default_yieldspins() = Threads.nthreads() == 1 ? -1 : 50
default_yieldspins() = -1

function WorkStealingScheduler(
    threadids::Vector{Int};
    unfairspins::Integer = -1,
    yieldspins::Integer = default_yieldspins(),
)
    workerindices = workerindices_from_threadids(threadids)
    ref = Ref{WorkStealingScheduler}()
    workers = create_workers(threadids, ref)
    queues = [WorkStealingDeque{Task}() for _ in threadids]
    scheduler = WorkStealingScheduler(
        unfairspins,
        yieldspins,
        workerindices,
        threadids,
        workers,
        queues,
        ConcurrentQueue{Task}(),
        false,
    )
    ref[] = scheduler
    return scheduler
end

function Base.push!(scheduler::WorkStealingScheduler, task::Task)
    threadid = Threads.threadid()
    workerid = scheduler.workerindices[threadid]
    workerid == 0 && errorwith(; scheduler, threadid) do io, (; scheduler, threadid)
        println(io, "push!(scheduler, task) called outside of worker threads")
        println(io, "called at: threadid = ", threadid)
        print(io, "scheduler ", objid(scheduler), " uses worker thread(s): ")
        join(io, scheduler.threadids, ", ", ", and")
    end
    push!(scheduler.queues[workerid], task)
    return scheduler
end

function pushat!(scheduler::WorkStealingScheduler, workerid::Integer, task::Task)
    push!(scheduler.queues[workerid], task)
    return scheduler
end

function get_next_task!(scheduler::WorkStealingScheduler, isretry::Bool)
    let y = maybepopfirst!(scheduler.externalqueue)
        if y isa Some
            return something(y)::Task
        end
    end
    threadid = Threads.threadid()
    workerid = scheduler.workerindices[threadid]
    if !isretry
        task = maybepop!(scheduler.queues[workerid])
        task === nothing || return something(task)
        length(scheduler.threadids) == 1 && return nothing
        for _ in 1:length(scheduler.threadids)
            w = rand(1:length(scheduler.threadids)-1)
            w += scheduler.threadids[w] == threadid
            task = maybepopfirst!(scheduler.queues[w])
            task === nothing || return something(task)
        end
    end
    for (tid, queue) in zip(scheduler.threadids, scheduler.queues)
        tid == threadid && continue
        task = maybepopfirst!(queue)
        task === nothing || return something(task)
    end
    return nothing
end

ext: jl


text:
using MakieCore

@time begin
    using MakieCore, GeometryBasics
    @time begin
        s = MakieCore.CairoScreen(500, 500)
        x = MakieCore.Scatter(rand(Point2f0, 10))
        MakieCore.draw_atomic(s, x)
    end
end


text:
# some trivial helpers

"""
    frequency_to_wavenumber(x)

Converts a given frequency in Hz to wavenumber in ``cm^{-1}``
"""
frequency_to_wavenumber(x::Real) = x / (1e2*c_c_SI)

"""
    wavelength_to_wavenumber(x)

Converts a given wavelength in m to wavenumber in ``cm^{-1}``
"""
wavelength_to_wavenumber(x::Real) = 1e-2 / x

# float tolerant leq method
≲(a, b) = (a <= b) || a ≈ b


text:
"""
    Bond

Represents a bond between atom indices `i` and `j`, with integer displacement of
`n[1] ... n[D]` unit cells along each dimension.
"""
struct Bond
    i :: Int
    j :: Int
    n :: SVector{3, Int}
end

"Represents a bond expressed as two fractional coordinates"
struct BondRaw
    ri::SVector{3, Float64}
    rj::SVector{3, Float64}
end

function Bond(cryst::Crystal, b::BondRaw)
    i = position_to_index(cryst, b.ri)
    j = position_to_index(cryst, b.rj)
    ri = cryst.positions[i]
    rj = cryst.positions[j]
    n = round.(Int, (b.rj-b.ri) - (rj-ri))
    return Bond(i, j, n)
end

function BondRaw(cryst::Crystal, b::Bond)
    return BondRaw(cryst.positions[b.i], cryst.positions[b.j]+b.n)
end

function Base.show(io::IO, mime::MIME"text/plain", bond::Bond)
    print(io, "Bond($(bond.i), $(bond.j), $(bond.n))")
end

function position_to_index(cryst::Crystal, r::Vec3)
    return findfirst(r′ -> is_same_position(r, r′; symprec=cryst.symprec), cryst.positions)
end

"""    displacement(cryst::Crystal, b::Bond)

The displacement vector ``𝐫_j - 𝐫_i`` in global coordinates between atoms
`b.i` and `b.j`, accounting for the integer offsets `b.n` between unit cells.
"""
function displacement(cryst::Crystal, b::BondRaw)
    return cryst.lat_vecs * (b.rj - b.ri)
end

function displacement(cryst::Crystal, b::Bond)
    return displacement(cryst, BondRaw(cryst, b))
end

"""    distance(cryst::Crystal, b::Bond)

The global distance between atoms in bond `b`. Equivalent to
`norm(displacement(cryst, b))`.
"""
function distance(cryst::Crystal, b::BondRaw)
    return norm(displacement(cryst, b))
end

function distance(cryst::Crystal, b::Bond)
    return norm(displacement(cryst, b))
end

function transform(s::SymOp, b::BondRaw)
    return BondRaw(transform(s, b.ri), transform(s, b.rj))
end

function transform(cryst::Crystal, s::SymOp, b::Bond)
    return Bond(cryst, transform(s, BondRaw(cryst, b)))
end

function reverse(b::BondRaw)
    return BondRaw(b.rj, b.ri)
end


text:
module FLICameras
using Printf
using Images
using Dates
using Base: @propagate_inbounds
using Base.Threads
using Distributed
import Base:
    VersionNumber,
    axes,
    copy,
    deepcopy,
    eachindex,
    eltype,
    empty!,
    fill!,
    firstindex,
    getindex,
    getproperty,
    isvalid,
    isreadable,
    iswritable,
    isequal,
    islocked,
    iterate,
    IndexStyle,
    last,
    lastindex,
    length,
    lock,
    ndims,
    propertynames,
    parent,
    parse,
    reset,
    reshape,
    size,
    show,
    showerror,
    similar,
    stride,
    setproperty!,
    setindex!,
    timedwait,
    trylock,
    unlock,
    wait

# TAO bindings
using Statistics
using ArrayTools
using ResizableArrays
import Base.Libc: TimeVal
using Base: @propagate_inbounds

# include dependents
begin deps = normpath(joinpath(@__DIR__, "../deps/deps.jl"))
    isfile(deps) || error(
        "File \"$deps\" does not exits, see \"README.md\" for installation.")
    include(deps)
end

# prepare files
function __init__()
    img_fname = "img_config.txt"
    path = "/tmp/FLICameras/"
    shmid_fname = "shmids.txt"
    try
        mkdir(path)
        touch(joinpath(path,img_fname))
        touch(joinpath(path,shmid_fname))
    catch InitError
        @warn "files already exist"
    finally
        @info "image configuration and shmids files are created in $path"
    end
    nothing
end

include("macros.jl")
include("types.jl")
include("errors.jl")
include("methods.jl")
include("cameras.jl")

end


text:
# SELECT clause.

struct SelectTop
    limit::Int
    with_ties::Bool

    SelectTop(; limit, with_ties = false) =
        new(limit, with_ties)
end

Base.convert(::Type{SelectTop}, limit::Integer) =
    SelectTop(limit = limit)

Base.convert(::Type{SelectTop}, t::NamedTuple) =
    SelectTop(; t...)

function PrettyPrinting.quoteof(t::SelectTop)
    if !t.with_ties
        t.limit
    else
        Expr(:tuple, Expr(:(=), :limit, t.limit),
                     Expr(:(=), :with_ties, t.with_ties))
    end
end

mutable struct SelectClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    top::Union{SelectTop, Nothing}
    distinct::Bool
    list::Vector{SQLClause}

    SelectClause(;
                 over = nothing,
                 top = nothing,
                 distinct = false,
                 list) =
        new(over, top, distinct, list)
end

SelectClause(list...; over = nothing, top = nothing, distinct = false) =
    SelectClause(over = over, top = top, distinct = distinct, list = SQLClause[list...])

"""
    SELECT(; over = nothing, top = nothing, distinct = false, list)
    SELECT(list...; over = nothing, top = nothing, distinct = false)

A `SELECT` clause.  Unlike raw SQL, `SELECT()` should be placed at the end of a
clause chain.

Set `distinct` to `true` to add a `DISTINCT` modifier.

# Examples

```jldoctest
julia> c = SELECT(true, false);

julia> print(render(c))
SELECT TRUE, FALSE
```

```jldoctest
julia> c = FROM(:location) |>
           SELECT(distinct = true, :zip);

julia> print(render(c))
SELECT DISTINCT "zip"
FROM "location"
```
"""
SELECT(args...; kws...) =
    SelectClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(SELECT), pats::Vector{Any}) =
    dissect(scr, SelectClause, pats)

function PrettyPrinting.quoteof(c::SelectClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(SELECT))
    if c.top !== nothing
        push!(ex.args, Expr(:kw, :top, quoteof(c.top)))
    end
    if c.distinct !== false
        push!(ex.args, Expr(:kw, :distinct, c.distinct))
    end
    if isempty(c.list)
        push!(ex.args, Expr(:kw, :list, Expr(:vect)))
    else
        append!(ex.args, quoteof(c.list, ctx))
    end
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::SelectClause, c′) =
    SelectClause(over = rebase(c.over, c′), top = c.top, distinct = c.distinct, list = c.list)


text:
export QuadDirect

struct QuadDirect end

function __solve(prob::OptimizationProblem, opt::QuadDirect; splits = nothing, maxiters = nothing, kwargs...)

    local x, _loss

    if !(isnothing(maxiters)) && maxiters <= 0.0
        error("The number of maxiters has to be a non-negative and non-zero number.")
    elseif !(isnothing(maxiters))
        maxiters = convert(Int, maxiters)
    end

    if splits === nothing
        error("You must provide the initial locations at which to evaluate the function in `splits` (a list of 3-vectors with values in strictly increasing order and within the specified bounds).")
    end

    _loss = function(θ)
        x = prob.f(θ, prob.p)
        return first(x)
    end

    t0 = time()

    root, x0 = !(isnothing(maxiters)) ? QuadDIRECT.analyze(_loss, splits, prob.lb, prob.ub; maxevals = maxiters, kwargs...) : QuadDIRECT.analyze(_loss, splits, prob.lb, prob.ub; kwargs...)
    box = minimum(root)
    t1 = time()

    SciMLBase.build_solution(prob, opt, QuadDIRECT.position(box, x0), QuadDIRECT.value(box); original=root)
end


text:
#=
 Train a simple network to predict formation energy per atom (downloaded from Materials Project). =#
using CSV, DataFrames
using Random, Statistics
using Flux
using Flux: @epochs, glorot_uniform, Zygote, gpu
using ChemistryFeaturization
using AtomicGraphNets
using Formatting

include("models.jl")

function train_formation_energy(;
    num_pts = 100,
    num_epochs = 25,
    data_dir = joinpath(@__DIR__, "data"),
    verbose = true,
    sample_size = 10,
    complexity_constant = 1e-6,
)
    println("Setting things up...")

    # data-related options
    train_frac = 0.8 # what fraction for training?
    num_train = Int32(round(train_frac * num_pts))
    num_test = num_pts - num_train
    prop = "formation_energy_per_atom"
    id = "task_id" # field by which to label each input material

    # set up the featurization
    featurization = GraphNodeFeaturization([
        "Group",
        "Row",
        "Block",
        "Atomic mass",
        "Atomic radius",
        "X",
    ])
    num_features =
        sum(ChemistryFeaturization.FeatureDescriptor.output_shape.(featurization.features)) # TODO: update this with cleaner syntax once new version of CF is tagged that has it

    # model hyperparameters – keeping it pretty simple for now
    num_conv = 3 # how many convolutional layers?
    crys_fea_len = 32 # length of crystal feature vector after pooling (keep node dimension constant for now)
    num_hidden_layers = 1 # how many fully-connected layers after convolution and pooling?

    # dataset...first, read in outputs
    info = CSV.read(string(data_dir, "/", prop, ".csv"), DataFrame)
    y = Array(Float32.(info[!, Symbol(prop)]))

    # shuffle data and pick out subset
    indices = shuffle(1:size(info, 1))[1:num_pts]
    info = info[indices, :]
    output = y[indices]

    # next, make and featurize graphs
    if verbose
        println("Building graphs and feature vectors from structures...")
    end
    inputs = FeaturizedAtoms[]

    # for r in eachrow(info)
    cifpaths = [
        joinpath(data_dir, format("{}_cifs", prop), string(r[Symbol(id)], ".cif")) for
        r in eachrow(info)
    ]
    graphs = skipmissing(AtomGraph.(cifpaths))
    inputs = [featurize(gr, featurization) for gr in graphs]

    # pick out train/test sets
    if verbose
        println("Dividing into train/test sets...")
    end
    train_output = output[1:num_train]
    test_output = output[num_train+1:end]
    train_input = inputs[1:num_train]
    test_input = inputs[num_train+1:end]
    train_data = zip(train_input, train_output)
    test_data = zip(test_input, test_output)

    # build the model
    if verbose
        println("Building the network...")
    end
    model = BayesianCGCNN(
        num_features,
        num_conv = num_conv,
        atom_conv_feature_length = crys_fea_len,
        pooled_feature_length = (Int(crys_fea_len / 2)),
        num_hidden_layers = 1,
    )
    model = gpu(model)
    ps = Flux.params(model)
    opt = ADAM(0.01) # optimizer

    # define loss function and a callback to monitor progress
    loss(x, y) = Flux.Losses.mse(model(x), y)

    if verbose
        println("Training!")
    end

    function test(test_data)
        total_loss = 0
        count = 0
        for (x, y) in test_data
            x, y = x |> gpu, y |> gpu
            total_loss += loss(x, y)
            count += 1
        end
        return (total_loss / count)
    end

    function kl_loss_calc(model)
        layers = Zygote.@ignore Flux.modules(model)
        return sum(normal_kl_divergence.(layers))
    end
    for epoch = 1:num_epochs
        for (x, y) in train_data
            x, y = x |> gpu, y |> gpu
            gs = Flux.gradient(ps) do
                for sample in sample_size
                    total_loss = loss(x, y)
                    kl_loss = kl_loss_calc(model)
                    total_loss += complexity_constant * kl_loss
                end
                total_loss /= sample_size
                return total_loss
            end
            Flux.Optimise.update!(opt, ps, gs)
        end
        test_loss = test(test_data)
        println(format("Test Loss: {}", test_loss))
    end

    return model
end

train_formation_energy()


text:
#!/usr/bin/env julia

# Unary vectorised math operations

fun_arr_arr = [(+), (.*), (./), (.^), hypot, min, mod]
fun_arr_arr_name = ["add", "mul", "div", "pow", "hypot", "min2", "fmod"]

# Binary vectorised math operations

function sigmoid(z)
    return 1.0 ./ (1.0 .+ exp.(-z))
end
function sqrt_a(z)
    return sqrt.(z)
end
function exp_a(z)
    return exp.(z)
end

fun_arr = [copy, abs, exp_a, log, sqrt_a, cbrt, sin, tan,
  asin, sinh, asinh, round, sort, sigmoid]
fun_arr_name = ["copy", "abs", "exp", "log", "sqrt", "cbrt", "sin", "tan",
  "asin", "sinh", "asinh", "round", "sort", "sigmoid"]

# Fold and scan operations

fun_axis_arr = [maximum, sum, prod, cumprod, cummax]
fun_axis_arr_name = ["max", "sum", "prod", "cumprod", "cummax"]

fun_axes_arr = [sum]
fun_axes_arr_name = ["sum_reduce"]

# Repeat operations

function rep(x, s)
  return repeat(x, inner=s)
end
function tile(x, s)
  return repeat(x, outer=s)
end
fun_repeat = [rep, tile]
fun_repeat_name = ["repeat", "tile"]

# Linear Algebra operation

function matm(x) return x * x end
fun_linalg = [matm, inv, eigvals, svd, lu, qr]
fun_linalg_name = ["matmul", "inv", "eigvals", "svd", "lu", "qr"]

# Timing functions

function time_fun(fn)
    return (@elapsed fn()) * 1000
end

function remove_outlier(arr)
    fp = quantile(arr, 0.25)
    tp = quantile(arr, 0.75)
    return filter(x -> (x >= fp) && (x <= tp), arr)
end

function timing(fn, msg)
    c = 30
    ts = Float32[]
    for i = 1:c
        t = fn()
        push!(ts, t)
    end
    ts = remove_outlier(ts)
    m = mean(ts)
    s = std(ts)
    @printf("| %s :\t mean = %.5f \t std = %.5f\n", msg, m, s)
    return m, s
end

function evalop_arr_arr(fn, name, sz)
    function f()
        inp1 = rand(Float32, sz)
        inp2 = rand(Float32, sz)
        discard = fn(inp1, inp2)
        function g()
            return fn(inp1, inp2)
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (%d)", name, sz))
end

function evalop_arr(fn, name, sz)
    function f()
        inp = rand(Float32, sz)
        discard = fn(inp)
        function g()
            return fn(inp)
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (%d)", name, sz))
end

function evalop_axis_arr(axis, fn, name, sz)
    function f()
        inp = rand(Float32, sz)
        discard = fn(inp, axis)
        function g()
            return fn(inp, axis)
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (axis=%d, %s)", name, axis, string(sz)))
end

function evalop_axes_arr(axis, fn, name, sz)
    function f()
        inp = rand(Float32, sz)
        discard = fn(inp, axis)
        function g()
            return fn(inp, axis)
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (axis=%s, %s)", name, string(axis), string(sz)))
end

function evalop_repeat(axes, fn, name, sz)
    function f()
        inp = rand(Float32, sz)
        discard = fn(inp, axes)
        function g()
            return fn(inp, axes)
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (axis=%s, %s)", name, string(axes), string(sz)))
end

# TODO : more general unpacking
function evalop_slice(idx, idx_str, sz)
    function f()
        inp = rand(Float32, sz)
        function g()
            i1, i2, i3 = idx
            return inp[i1, i2, i3]
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (%s)", "get_slice", idx_str))
end

function evalop_linalg(fn, name, sz)
    function f()
        inp = rand(Float32, sz)
        discard = fn(inp)
        function g()
            return fn(inp)
        end
        return time_fun(g)
    end
    return timing(f, @sprintf("%s (%d*%d)", name, sz[1], sz[2]))
end

# Evaluate simple arr and arr_arr operations

function evaluate_simple()
    sz = [10, 100, 1000, 10000, 100000, 200000, 400000, 600000, 800000, 1000000]
    sz_str = "10,,100,,1000,,1e4,,1e5,,2e5,,4e5,,6e5,,8e5,,1e6"

    result_str = "," * sz_str * "\n"
    for i = 1:length(fun_arr)
        result_str *= fun_arr_name[i] * ","
        for j = 1:length(sz)
            mu, std = evalop_arr(fun_arr[i], fun_arr_name[i], sz[j])
            result_str *= @sprintf("%.4f, %.4f,", mu, std)
        end
        result_str *= "\n"
    end

    for i = 1:length(fun_arr_arr)
        result_str *= fun_arr_arr_name[i] * ","
        for j = 1:length(sz)
            mu, std = evalop_arr_arr(fun_arr_arr[i], fun_arr_arr_name[i], sz[j])
            result_str *= @sprintf("%.4f, %.4f,", mu, std)
        end
        result_str *= "\n"
    end
    return result_str
end

# Evaluate axis operations

function evaluate_axis()
    sz = [(10, 10, 10, 10), (20, 20, 20, 20), (30, 30, 30, 30),
        (40, 40, 40, 40), (50, 50, 50, 50), (60, 60, 60, 60)]
    sz_str = "10,,20,,30,,40,,50,,60"
    axis = [1,4]
    result_str = "," * sz_str * "\n"
    for i = 1:length(fun_axis_arr)
        for k = 1:length(axis)
            result_str *= @sprintf("%s(axis=%d),", fun_axis_arr_name[i], (axis[k] - 1))
            for j = 1:length(sz)
                mu, std = evalop_axis_arr(axis[k], fun_axis_arr[i],
                    fun_axis_arr_name[i], sz[j])
                result_str *= @sprintf("%.4f, %.4f,", mu, std)
            end
            result_str *= "\n"
        end
    end
    return result_str
end

# Evaluate axes operations

function evaluate_axes()
    sz = [(10, 10, 10, 10), (20, 20, 20, 20), (30, 30, 30, 30),
        (40, 40, 40, 40), (50, 50, 50, 50), (60, 60, 60, 60), (70, 70, 70, 70)]
    sz_str = "10,,20,,30,,40,,50,,60,,70"
    axis = [(1,4), (1,3)]
    axis_for_str = [(0,3), (0,2)]
    result_str = "," * sz_str * "\n"
    for i = 1:length(fun_axes_arr)
        for k = 1:length(axis)
            axes_str = join(axis_for_str[k], '*')
            result_str *= @sprintf("%s(axes=%s),", fun_axes_arr_name[i], axes_str)
            for j = 1:length(sz)
                mu, std = evalop_axes_arr(axis[k], fun_axes_arr[i],
                    fun_axes_arr_name[i], sz[j])
                result_str *= @sprintf("%.4f, %.4f,", mu, std)
            end
            result_str *= "\n"
        end
    end
    return result_str
end

# Evaluate repeat operations

function evaluate_repeat()
    sz = [(10, 10, 10, 10), (15, 15, 15, 15), (20, 20, 20, 20),
        (25, 25, 25, 25), (30, 30, 30, 30), (35, 35, 35, 35)]
    sz_str = "10,,15,,20,,25,,30,,35"
    axes = [(1,1,1,5), (1,4,4,1), (3,3,3,1), (2,2,2,2)]

    result_str = "," * sz_str * "\n"
    for i = 1:length(fun_repeat)
        for k = 1:length(axes)
            axes_str = join(axes[k], '*')
            result_str *=  @sprintf("%s(axes=%s),", fun_repeat_name[i], axes_str)
            for j = 1:length(sz)
                mu, std = evalop_repeat(axes[k], fun_repeat[i], fun_repeat_name[i], sz[j])
                result_str *= @sprintf("%.4f, %.4f,", mu, std)
            end
            result_str *= "\n"
        end
    end
    return result_str
end

# Evaluate repeat operations

# TODO: use general unpack operations; the index is not flexible
function evaluate_slice()
    sz = [(10, 300, 3000), (3000, 300, 10)]
    sz_str = "10*300*3000,,3000*300*10"

    index1 = [
        [1:10, :, :],
        [10:-1:1, 1:2, :],
        [10:-1:1, 300:-1:1, 1],
        [10, 300:-1:1, :],
        [:, 300:-1:1, :],
        [:, 1:300, 3000:-1:1],
        [:, 300:-1:1, 1:2],
        [:, 1:300, 3000:-2:1]
    ]
    index2 = [
        [1:3000, :, :],
        [3000:-1:1, 1:2, :],
        [3000:-1:1, 300:-1:1, 1],
        [3000, 300:-1:1, :],
        [:, 300:-1:1, :],
        [:, 1:300, 10:-1:1],
        [:, 300:-1:1, 1:2],
        [:, 1:300, 10:-2:1]
    ]
    index = [index1, index2]

    index_str = [
        "[[0;-1]; []; []]", "[[-1;0]; [0;1]; []]",
        "[[-1;0]; [-1;0]; [0]]", "[[-1]; [-1;0];[]]",
        "[[]; [-1;0]; []]", "[[]; [0;-1]; [-1;0]]",
        "[[]; [-1;0]; [0;1]]", "[[]; [0;-1]; [-1;0;-2]]"]

    result_str = "," * sz_str * "\n"
    for i = 1:length(index[1])
        result_str *= @sprintf("%s(index=%s),", "get_slice", index_str[i])
        for j = 1:length(sz)
            mu, std = evalop_slice(index[j][i], index_str[i], sz[j])
            result_str *= @sprintf("%.4f, %.4f,", mu, std)
        end
        result_str *= "\n"
    end
    return result_str
end

# Evaluate linear algebra operations

function evaluate_linalg()
    sz = [(10, 10), (50, 50), (100, 100),
        (150, 150), (200, 200), (300, 300), (400, 400),
        (600, 600), (800, 800), (1000, 1000)]
    sz_str = "10,,50,,100,,150,,200,,300,,400,,600,,800,,1000"

    result_str = "," * sz_str * "\n"
    for i = 1:length(fun_linalg)
        result_str *= @sprintf("%s,", fun_linalg_name[i])
        for j = 1:length(sz)
            mu, std = evalop_linalg(fun_linalg[i],
                fun_linalg_name[i], sz[j])
            result_str *= @sprintf("%.4f, %.4f,", mu, std)
        end
        result_str *= "\n"
    end
    return result_str
end

function write_file(fname, output_str)
    open(fname, "w") do f
        write(f, output_str)
    end
end

write_file("simple_julia.csv", evaluate_simple())
write_file("axis_julia.csv",   evaluate_axis())
write_file("axes_julia.csv",   evaluate_axes())
# write_file("repeat_julia.csv", evaluate_repeat())
write_file("slice_julia.csv",  evaluate_slice())
write_file("linalg_julia.csv", evaluate_linalg())

max_stars_repo_name: owlbarn/benchmark

max_issues_repo_name: owlbarn/benchmark

max_forks_repo_name: owlbarn/benchmark


text:
# Julia wrapper for header: laszip_api.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const laszip = replace(liblaszip, "_api" => "")
@info "Using laszip library at $liblaszip"

# function renameat(arg1::Cint, arg2, arg3::Cint, arg4)
#     ccall((:renameat, laszip), Cint, (Cint, Cstring, Cint, Cstring), arg1, arg2, arg3, arg4)
# end

# function renamex_np(arg1, arg2, arg3::UInt32)
#     ccall((:renamex_np, laszip), Cint, (Cstring, Cstring, UInt32), arg1, arg2, arg3)
# end

# function renameatx_np(arg1::Cint, arg2, arg3::Cint, arg4, arg5::UInt32)
#     ccall((:renameatx_np, laszip), Cint, (Cint, Cstring, Cint, Cstring, UInt32), arg1, arg2, arg3, arg4, arg5)
# end

# function clearerr(arg1)
#     ccall((:clearerr, laszip), Nothing, (Ptr{Cvoid},), arg1)
# end

# function fclose(arg1)
#     ccall((:fclose, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function feof(arg1)
#     ccall((:feof, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function ferror(arg1)
#     ccall((:ferror, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function fflush(arg1)
#     ccall((:fflush, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function fgetc(arg1)
#     ccall((:fgetc, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function fgetpos(arg1, arg2)
#     ccall((:fgetpos, laszip), Cint, (Ptr{Cvoid}, Ptr{Int64}), arg1, arg2)
# end

# function fgets(arg1, arg2::Cint, arg3)
#     ccall((:fgets, laszip), Cstring, (Cstring, Cint, Ptr{Cvoid}), arg1, arg2, arg3)
# end

# function fopen(__filename, __mode)
#     ccall((:fopen, laszip), Ptr{Cvoid}, (Cstring, Cstring), __filename, __mode)
# end

# function fputc(arg1::Cint, arg2)
#     ccall((:fputc, laszip), Cint, (Cint, Ptr{Cvoid}), arg1, arg2)
# end

# function fputs(arg1, arg2)
#     ccall((:fputs, laszip), Cint, (Cstring, Ptr{Cvoid}), arg1, arg2)
# end

# function fread(__ptr, __size::Csize_t, __nitems::Csize_t, __stream)
#     ccall((:fread, laszip), Csize_t, (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}), __ptr, __size, __nitems, __stream)
# end

# function freopen(arg1, arg2, arg3)
#     ccall((:freopen, laszip), Ptr{Cvoid}, (Cstring, Cstring, Ptr{Cvoid}), arg1, arg2, arg3)
# end

# function fseek(arg1, arg2::Clong, arg3::Cint)
#     ccall((:fseek, laszip), Cint, (Ptr{Cvoid}, Clong, Cint), arg1, arg2, arg3)
# end

# function fsetpos(arg1, arg2)
#     ccall((:fsetpos, laszip), Cint, (Ptr{Cvoid}, Ptr{Int64}), arg1, arg2)
# end

# function ftell(arg1)
#     ccall((:ftell, laszip), Clong, (Ptr{Cvoid},), arg1)
# end

# function fwrite(__ptr, __size::Csize_t, __nitems::Csize_t, __stream)
#     ccall((:fwrite, laszip), Csize_t, (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}), __ptr, __size, __nitems, __stream)
# end

# function getc(arg1)
#     ccall((:getc, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function getchar()
#     ccall((:getchar, laszip), Cint, ())
# end

# function gets(arg1)
#     ccall((:gets, laszip), Cstring, (Cstring,), arg1)
# end

# function perror(arg1)
#     ccall((:perror, laszip), Nothing, (Cstring,), arg1)
# end

# function putc(arg1::Cint, arg2)
#     ccall((:putc, laszip), Cint, (Cint, Ptr{Cvoid}), arg1, arg2)
# end

# function putchar(arg1::Cint)
#     ccall((:putchar, laszip), Cint, (Cint,), arg1)
# end

# function puts(arg1)
#     ccall((:puts, laszip), Cint, (Cstring,), arg1)
# end

# function remove(arg1)
#     ccall((:remove, laszip), Cint, (Cstring,), arg1)
# end

# function rename(__old, __new)
#     ccall((:rename, laszip), Cint, (Cstring, Cstring), __old, __new)
# end

# function rewind(arg1)
#     ccall((:rewind, laszip), Nothing, (Ptr{Cvoid},), arg1)
# end

# function setbuf(arg1, arg2)
#     ccall((:setbuf, laszip), Nothing, (Ptr{Cvoid}, Cstring), arg1, arg2)
# end

# function setvbuf(arg1, arg2, arg3::Cint, arg4::Csize_t)
#     ccall((:setvbuf, laszip), Cint, (Ptr{Cvoid}, Cstring, Cint, Csize_t), arg1, arg2, arg3, arg4)
# end

# function tmpfile()
#     ccall((:tmpfile, laszip), Ptr{Cvoid}, ())
# end

# function tmpnam(arg1)
#     ccall((:tmpnam, laszip), Cstring, (Cstring,), arg1)
# end

# function ungetc(arg1::Cint, arg2)
#     ccall((:ungetc, laszip), Cint, (Cint, Ptr{Cvoid}), arg1, arg2)
# end

# function vfprintf(arg1, arg2, arg3::Cstring)
#     ccall((:vfprintf, laszip), Cint, (Ptr{Cvoid}, Cstring, Cstring), arg1, arg2, arg3)
# end

# function vprintf(arg1, arg2)
#     ccall((:vprintf, laszip), Cint, (Cstring, Ptr{__va_list_tag}), arg1, arg2)
# end

# function vsprintf(arg1, arg2, arg3)
#     ccall((:vsprintf, laszip), Cint, (Cstring, Cstring, Ptr{__va_list_tag}), arg1, arg2, arg3)
# end

# function ctermid(arg1)
#     ccall((:ctermid, laszip), Cstring, (Cstring,), arg1)
# end

# function fdopen(arg1::Cint, arg2)
#     ccall((:fdopen, laszip), Ptr{Cvoid}, (Cint, Cstring), arg1, arg2)
# end

# function fileno(arg1)
#     ccall((:fileno, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function pclose(arg1)
#     ccall((:pclose, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function popen(arg1, arg2)
#     ccall((:popen, laszip), Ptr{Cvoid}, (Cstring, Cstring), arg1, arg2)
# end

# function flockfile(arg1)
#     ccall((:flockfile, laszip), Nothing, (Ptr{Cvoid},), arg1)
# end

# function ftrylockfile(arg1)
#     ccall((:ftrylockfile, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function funlockfile(arg1)
#     ccall((:funlockfile, laszip), Nothing, (Ptr{Cvoid},), arg1)
# end

# function getc_unlocked(arg1)
#     ccall((:getc_unlocked, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function getchar_unlocked()
#     ccall((:getchar_unlocked, laszip), Cint, ())
# end

# function putc_unlocked(arg1::Cint, arg2)
#     ccall((:putc_unlocked, laszip), Cint, (Cint, Ptr{Cvoid}), arg1, arg2)
# end

# function putchar_unlocked(arg1::Cint)
#     ccall((:putchar_unlocked, laszip), Cint, (Cint,), arg1)
# end

# function getw(arg1)
#     ccall((:getw, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function putw(arg1::Cint, arg2)
#     ccall((:putw, laszip), Cint, (Cint, Ptr{Cvoid}), arg1, arg2)
# end

# function tempnam(__dir, __prefix)
#     ccall((:tempnam, laszip), Cstring, (Cstring, Cstring), __dir, __prefix)
# end

# function fseeko(__stream, __offset::Int64, __whence::Cint)
#     ccall((:fseeko, laszip), Cint, (Ptr{Cvoid}, Int64, Cint), __stream, __offset, __whence)
# end

# function ftello(__stream)
#     ccall((:ftello, laszip), Int64, (Ptr{Cvoid},), __stream)
# end

# function vfscanf(__stream, __format, arg1)
#     ccall((:vfscanf, laszip), Cint, (Ptr{Cvoid}, Cstring, Ptr{__va_list_tag}), __stream, __format, arg1)
# end

# function vscanf(__format, arg1)
#     ccall((:vscanf, laszip), Cint, (Cstring, Ptr{__va_list_tag}), __format, arg1)
# end

# function vsnprintf(__str, __size::Culong, __format, arg1)
#     ccall((:vsnprintf, laszip), Cint, (Cstring, Culong, Cstring, Ptr{__va_list_tag}), __str, __size, __format, arg1)
# end

# function vsscanf(__str, __format, arg1)
#     ccall((:vsscanf, laszip), Cint, (Cstring, Cstring, Ptr{__va_list_tag}), __str, __format, arg1)
# end

# function vdprintf(arg1::Cint, arg2, arg3::Cstring)
#     ccall((:vdprintf, laszip), Cint, (Cint, Cstring, Cstring), arg1, arg2, arg3)
# end

# function getdelim(__linep, __linecapp, __delimiter::Cint, __stream)
#     ccall((:getdelim, laszip), Int32, (Ptr{Cstring}, Ptr{Csize_t}, Cint, Ptr{Cvoid}), __linep, __linecapp, __delimiter, __stream)
# end

# function getline(__linep, __linecapp, __stream)
#     ccall((:getline, laszip), Int32, (Ptr{Cstring}, Ptr{Csize_t}, Ptr{Cvoid}), __linep, __linecapp, __stream)
# end

# function ctermid_r(arg1)
#     ccall((:ctermid_r, laszip), Cstring, (Cstring,), arg1)
# end

# function fgetln(arg1, arg2)
#     ccall((:fgetln, laszip), Cstring, (Ptr{Cvoid}, Ptr{Csize_t}), arg1, arg2)
# end

# function fmtcheck(arg1, arg2)
#     ccall((:fmtcheck, laszip), Cstring, (Cstring, Cstring), arg1, arg2)
# end

# function fpurge(arg1)
#     ccall((:fpurge, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function setbuffer(arg1, arg2, arg3::Cint)
#     ccall((:setbuffer, laszip), Nothing, (Ptr{Cvoid}, Cstring, Cint), arg1, arg2, arg3)
# end

# function setlinebuf(arg1)
#     ccall((:setlinebuf, laszip), Cint, (Ptr{Cvoid},), arg1)
# end

# function vasprintf(arg1, arg2, arg3::Cstring)
#     ccall((:vasprintf, laszip), Cint, (Ptr{Cstring}, Cstring, Cstring), arg1, arg2, arg3)
# end

# function zopen(arg1, arg2, arg3::Cint)
#     ccall((:zopen, laszip), Ptr{Cvoid}, (Cstring, Cstring, Cint), arg1, arg2, arg3)
# end

# function funopen(arg1, arg2, arg3, arg4, arg5)
#     ccall((:funopen, laszip), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5)
# end

function laszip_get_version(version_major, version_minor, version_revision, version_build)
    ccall((:laszip_get_version, laszip), Int32, (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt16}, Ptr{UInt32}), version_major, version_minor, version_revision, version_build)
end

function laszip_create(pointer::Ref{Ptr{Cvoid}})
    ccall((:laszip_create, laszip), Int32, (Ptr{Ptr{Cvoid}},), pointer)
end

function laszip_get_error(pointer::Ptr{Cvoid}, error)
    ccall((:laszip_get_error, laszip), Int32, (Ptr{Cvoid}, Ptr{Cstring}), pointer, error)
end

function laszip_get_warning(pointer::Ptr{Cvoid}, warning)
    ccall((:laszip_get_warning, laszip), Int32, (Ptr{Cvoid}, Ptr{Cstring}), pointer, warning)
end

function laszip_clean(pointer::Ptr{Cvoid})
    ccall((:laszip_clean, laszip), Int32, (Ptr{Cvoid},), pointer)
end

function laszip_destroy(pointer::Ptr{Cvoid})
    ccall((:laszip_destroy, laszip), Int32, (Ptr{Cvoid},), pointer)
end

function laszip_get_header_pointer(pointer::Ptr{Cvoid}, header_pointer)
    ccall((:laszip_get_header_pointer, laszip), Int32, (Ptr{Cvoid}, Ref{Ptr{LazHeader}}), pointer, header_pointer)
end

function laszip_get_point_pointer(pointer::Ptr{Cvoid}, point_pointer)
    ccall((:laszip_get_point_pointer, laszip), Int32, (Ptr{Cvoid}, Ref{Ptr{LazPoint}}), pointer, point_pointer)
end

function laszip_get_point_count(pointer::Ptr{Cvoid}, count)
    ccall((:laszip_get_point_count, laszip), Int32, (Ptr{Cvoid}, Ptr{Int64}), pointer, count)
end

function laszip_set_header(pointer::Ptr{Cvoid}, header)
    ccall((:laszip_set_header, laszip), Int32, (Ptr{Cvoid}, Ptr{LazHeader}), pointer, header)
end

# function laszip_set_point_type_and_size(pointer::Ptr{Cvoid}, point_type::UInt8, point_size::UInt16)
#     ccall((:laszip_set_point_type_and_size, laszip), Int32, (Ptr{Cvoid}, UInt8, UInt16), pointer, point_type, point_size)
# end

# function laszip_check_for_integer_overflow(pointer::Ptr{Cvoid})
#     ccall((:laszip_check_for_integer_overflow, laszip), Int32, (Ptr{Cvoid},), pointer)
# end

# function laszip_auto_offset(pointer::Ptr{Cvoid})
#     ccall((:laszip_auto_offset, laszip), Int32, (Ptr{Cvoid},), pointer)
# end

function laszip_set_point(pointer::Ptr{Cvoid}, point)
    ccall((:laszip_set_point, laszip), Int32, (Ptr{Cvoid}, Ptr{LazPoint}), pointer, point)
end

# function laszip_set_coordinates(pointer::Ptr{Cvoid}, coordinates)
#     ccall((:laszip_set_coordinates, laszip), Int32, (Ptr{Cvoid}, Ptr{Float64}), pointer, coordinates)
# end

# function laszip_get_coordinates(pointer::Ptr{Cvoid}, coordinates)
#     ccall((:laszip_get_coordinates, laszip), Int32, (Ptr{Cvoid}, Ptr{Float64}), pointer, coordinates)
# end

# function laszip_set_geokeys(pointer::Ptr{Cvoid}, number::UInt32, key_entries)
#     ccall((:laszip_set_geokeys, laszip), Int32, (Ptr{Cvoid}, UInt32, Ptr{Cvoid}), pointer, number, key_entries)
# end

# function laszip_set_geodouble_params(pointer::Ptr{Cvoid}, number::UInt32, geodouble_params)
#     ccall((:laszip_set_geodouble_params, laszip), Int32, (Ptr{Cvoid}, UInt32, Ptr{Float64}), pointer, number, geodouble_params)
# end

# function laszip_set_geoascii_params(pointer::Ptr{Cvoid}, number::UInt32, geoascii_params)
#     ccall((:laszip_set_geoascii_params, laszip), Int32, (Ptr{Cvoid}, UInt32, Ptr{UInt8}), pointer, number, geoascii_params)
# end

# function laszip_add_attribute(pointer::Ptr{Cvoid}, _type::UInt32, name, description, scale::Float64, offset::Float64)
#     ccall((:laszip_add_attribute, laszip), Int32, (Ptr{Cvoid}, UInt32, Ptr{UInt8}, Ptr{UInt8}, Float64, Float64), pointer, _type, name, description, scale, offset)
# end

# function laszip_add_vlr(pointer::Ptr{Cvoid}, user_id, record_id::UInt16, record_length_after_header::UInt16, description, data)
#     ccall((:laszip_add_vlr, laszip), Int32, (Ptr{Cvoid}, Ptr{UInt8}, UInt16, UInt16, Ptr{UInt8}, Ptr{UInt8}), pointer, user_id, record_id, record_length_after_header, description, data)
# end

# function laszip_remove_vlr(pointer::Ptr{Cvoid}, user_id, record_id::UInt16)
#     ccall((:laszip_remove_vlr, laszip), Int32, (Ptr{Cvoid}, Ptr{UInt8}, UInt16), pointer, user_id, record_id)
# end

# function laszip_create_spatial_index(pointer::Ptr{Cvoid}, create::Cint, append::Cint)
#     ccall((:laszip_create_spatial_index, laszip), Int32, (Ptr{Cvoid}, Cint, Cint), pointer, create, append)
# end

# function laszip_preserve_generating_software(pointer::Ptr{Cvoid}, preserve::Cint)
#     ccall((:laszip_preserve_generating_software, laszip), Int32, (Ptr{Cvoid}, Cint), pointer, preserve)
# end

# function laszip_request_native_extension(pointer::Ptr{Cvoid}, request::Cint)
#     ccall((:laszip_request_native_extension, laszip), Int32, (Ptr{Cvoid}, Cint), pointer, request)
# end

# function laszip_request_compatibility_mode(pointer::Ptr{Cvoid}, request::Cint)
#     ccall((:laszip_request_compatibility_mode, laszip), Int32, (Ptr{Cvoid}, Cint), pointer, request)
# end

# function laszip_set_chunk_size(pointer::Ptr{Cvoid}, chunk_size::UInt32)
#     ccall((:laszip_set_chunk_size, laszip), Int32, (Ptr{Cvoid}, UInt32), pointer, chunk_size)
# end

function laszip_open_writer(pointer::Ptr{Cvoid}, file_name, compress::Cint)
    ccall((:laszip_open_writer, laszip), Int32, (Ptr{Cvoid}, Ptr{UInt8}, Cint), pointer, file_name, compress)
end

function laszip_write_point(pointer::Ptr{Cvoid})
    ccall((:laszip_write_point, laszip), Int32, (Ptr{Cvoid},), pointer)
end

# function laszip_write_indexed_point(pointer::Ptr{Cvoid})
#     ccall((:laszip_write_indexed_point, laszip), Int32, (Ptr{Cvoid},), pointer)
# end

function laszip_update_inventory(pointer::Ptr{Cvoid})
    ccall((:laszip_update_inventory, laszip), Int32, (Ptr{Cvoid},), pointer)
end

function laszip_close_writer(pointer::Ptr{Cvoid})
    ccall((:laszip_close_writer, laszip), Int32, (Ptr{Cvoid},), pointer)
end

# function laszip_exploit_spatial_index(pointer::Ptr{Cvoid}, exploit::Cint)
#     ccall((:laszip_exploit_spatial_index, laszip), Int32, (Ptr{Cvoid}, Cint), pointer, exploit)
# end

# function laszip_decompress_selective(pointer::Ptr{Cvoid}, decompress_selective::UInt32)
#     ccall((:laszip_decompress_selective, laszip), Int32, (Ptr{Cvoid}, UInt32), pointer, decompress_selective)
# end

function laszip_open_reader(pointer::Ptr{Cvoid}, file_name, is_compressed)
    ccall((:laszip_open_reader, laszip), Int32, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Cint}), pointer, file_name, is_compressed)
end

# function laszip_has_spatial_index(pointer::Ptr{Cvoid}, is_indexed, is_appended)
#     ccall((:laszip_has_spatial_index, laszip), Int32, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}), pointer, is_indexed, is_appended)
# end

# function laszip_inside_rectangle(pointer::Ptr{Cvoid}, min_x::Float64, min_y::Float64, max_x::Float64, max_y::Float64, is_empty)
#     ccall((:laszip_inside_rectangle, laszip), Int32, (Ptr{Cvoid}, Float64, Float64, Float64, Float64, Ptr{Cint}), pointer, min_x, min_y, max_x, max_y, is_empty)
# end

function laszip_seek_point(pointer::Ptr{Cvoid}, index::Int64)
    ccall((:laszip_seek_point, laszip), Int32, (Ptr{Cvoid}, Int64), pointer, index)
end
laszip_seek_point(pointer::Ptr{Cvoid}, index::Int32) = laszip_seek_point(pointer, convert(Int64, index))

function laszip_read_point(pointer::Ptr{Cvoid})
    ccall((:laszip_read_point, laszip), Int32, (Ptr{Cvoid},), pointer)
end

# function laszip_read_inside_point(pointer::Ptr{Cvoid}, is_done)
#     ccall((:laszip_read_inside_point, laszip), Int32, (Ptr{Cvoid}, Ptr{Cint}), pointer, is_done)
# end

function laszip_close_reader(pointer::Ptr{Cvoid})
    ccall((:laszip_close_reader, laszip), Int32, (Ptr{Cvoid},), pointer)
end


text:
using DifferentialEquations
using Plots
using CPUTime

function main()
     a = 3.0  # パラメータ

     """
     微分方程式
     """
     dy(y, t, a) = a * y
     y0 = 1.0  # 初期値を設定
     tspan = (0.0, 1.0)  # 解く範囲を設定

     prob = ODEProblem(dy , y0 , tspan)
     sol = solve(prob)  # ソルバで解く

     # plot
     plot(
          sol,
          linewidth=5,
          xaxis="x",
          yaxis="y(x)",
          label="solution"
     )

     # savefig("solution_of_example_2.png")  #  保存
end

"""
オイラー法

・参考にしました -> [url](https://twitter.com/genkuroki/status/1301832571131633665/photo/1)
"""
function jisaku_solve_euler(dx, x₀, t_span, Δt)
     t = range(t_span..., step = Δt)  # 時間軸
     x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

     x[1] = x₀  # 初期値
     for i in 1:length(x)-1
          x[i+1] = x[i] + dx(t[i], x[i])Δt
     end

     t, x
end

"""ルンゲクッタ法（4次）"""
function jisaku_solve_RungeKutta(dx, x₀, t_span, Δt)
     t = range(t_span..., step = Δt)  # 時間軸
     x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

     x[1] = x₀  # 初期値
     for i in 1:length(x)-1
          k₁ = dx(t[i], x[i])
          k₂ = dx(t[i]+Δt/2, x[i]+k₁*Δt/2)
          k₃ = dx(t[i]+Δt/2, x[i]+k₂*Δt/2)
          k₄ = dx(t[i]+Δt, x[i]+k₃*Δt)
          x[i+1] = x[i] + (k₁ + 2k₂ + 2k₃ +k₄)Δt/6
     end

     t, x
end

function main()
     a = 3.0  # パラメータ
     dx(t, x) = a * x  # 微分方程式

     x₀ = 1.0
     t_span = (0.0, 5.0)
     Δt = 0.01
     t, x = jisaku_solve_RungeKutta(dx, x₀, t_span, Δt)
     plot(t, x, label="RungeKutta")

     t, x = jisaku_solve_euler(dx, x₀, t_span, Δt)
     plot!(t, x, label="euler")

     x = exp.(a*t)
     plot!(t, x, label="exact sol exp(at)")
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Images

import Base.Order: Ordering, ForwardOrdering, ReverseOrdering
import Base: ==, .==, +, -, *, /, .+, .-, .*, ./, .^, .<, .>
import Base: atan2, clamp, convert, copy, copy!, ctranspose, delete!, done, eltype,
             fft, float32, float64, get, getindex, haskey, hypot, ifft, length, map, map!,
             maximum, mimewritable, minimum, next, ndims, one, parent, permutedims, reinterpret,
             setindex!, show, showcompact, similar, size, slice, sqrt, squeeze,
             start, strides, sub, sum, write, writemime, zero
# "deprecated imports" are below

using Color, FixedPointNumbers, Compat
import Color: Fractional
if VERSION < v"0.4.0-dev+3275"
    using Base.Graphics
    import Base.Graphics: width, height, Point
else
    using Graphics
    import Graphics: width, height, Point
end
import FixedPointNumbers: ufixed8, ufixed10, ufixed12, ufixed14, ufixed16

using Base.Cartesian
include("compatibility/forcartesian.jl")

# if isdefined(module_parent(Images), :Grid)
#     import ..Grid.restrict
# end

const is_little_endian = ENDIAN_BOM == 0x04030201
immutable TypeConst{N} end  # for passing compile-time constants to functions

include("colortypes.jl")
using .ColorTypes

include("core.jl")
include("map.jl")
include("overlays.jl")
const have_imagemagick = include("ioformats/libmagickwand.jl")
@osx_only include("ioformats/OSXnative.jl")
include("io.jl")
include("labeledarrays.jl")
include("algorithms.jl")
include("connected.jl")
include("edge.jl")

__init__() = LibMagick.init()

function precompile()
    for T in (Uint8, Uint16, Int, Float32, Float64)
        Tdiv = typeof(one(T)/2)
        for N = 2:3
            precompile(restrict!, (Array{Tdiv, N}, Array{T,N}, Int))
            precompile(imfilter, (Array{T,N}, Array{Float64,N}))
            precompile(imfilter, (Array{T,N}, Array{Float64,2}))
            precompile(imfilter, (Array{T,N}, Array{Float32,N}))
            precompile(imfilter, (Array{T,N}, Array{Float32,2}))
        end
    end
    for T in (Float32, Float64)
        for N = 2:3
            precompile(_imfilter_gaussian!, (Array{T,N}, Vector{Float64}))
            precompile(_imfilter_gaussian!, (Array{T,N}, Vector{Int}))
            precompile(imfilter_gaussian_no_nans!, (Array{T,N}, Vector{Float64}))
            precompile(imfilter_gaussian_no_nans!, (Array{T,N}, Vector{Int}))
            precompile(fft, (Array{T,N},))
        end
    end
    for T in (Complex{Float32}, Complex{Float64})
        for N = 2:3
            precompile(ifft, (Array{T,N},))
        end
    end
end

export # types
    Gray,
    GrayAlpha,
    ARGB,
    BGRA,
    RGBA,
    RGB1,
    RGB4,
    YIQ,

    AbstractImage,
    AbstractImageDirect,
    AbstractImageIndexed,
    Image,
    ImageCmap,
    BitShift,
    ClampMin,
    ClampMax,
    ClampMinMax,
    Clamp,
    LabeledArray,
    MapInfo,
    MapNone,
    Overlay,
    OverlayImage,
    ScaleAutoMinMax,
    ScaleMinMax,
    ScaleSigned,
    SliceData,

    # constants
    palette_fire,
    palette_gray32,
    palette_gray64,
    palette_rainbow,

    # core functions
    assert2d,
    assert_scalar_color,
    assert_timedim_last,
    assert_xfirst,
    assert_yfirst,
    colordim,
    colorspace,
    coords,
    coords_spatial,
    copyproperties,
    data,
    dimindex,
    dimindexes,
    getindexim,
    grayim,
    colorim,
    height,
    isdirect,
    isxfirst,
    isyfirst,
    limits,
    maxabsfinite,
    maxfinite,
    minfinite,
    nchannels,
    ncolorelem,
    nimages,
    pixelspacing,
    properties,
    rerange!,
    reslice!,
    restrict,
    sdims,
    size_spatial,
    shareproperties,
    sliceim,
    spacedirections,
    spatialorder,
    spatialpermutation,
    spatialproperties,
    storageorder,
    subim,
    timedim,
    width,
    widthheight,
    raw,

    # io functions
    add_image_file_format,
    imread,
    imwrite,
    loadformat,

    # iterator functions
    first_index,
    iterate_spatial,
    parent,

    # color-related functions
    alphaval,
    blueval,
    greenval,
    redval,
    red,
    green,
    blue,
    hsi2rgb,
    imadjustintensity,
    indexedcolor,
    lut,
    rgb2hsi,
    rgb2ycbcr,
    separate,
    uint32color,
    uint32color!,
    ycbcr2rgb,

    # Scaling of intensity
    sc,
    scale,
    mapinfo,
    uint8sc,
    uint16sc,
    uint32sc,

    # algorithms
    ando3,
    ando4,
    ando5,
    backdiffx,
    backdiffy,
    dilate,
    erode,
    opening,
    closing,
    forwarddiffx,
    forwarddiffy,
    gaussian2d,
    imaverage,
    imcomplement,
    imdog,
    imedge,
    imfilter,
    imfilter!,
    imfilter_fft,
    imfilter_gaussian,
    imfilter_gaussian!,
    imfilter_LoG,
    imgaussiannoise,
    imgradients,
    imlaplacian,
    imlineardiffusion,
    imlog,
    imROF,
    imstretch,
#     imthresh,
    label_components,
    label_components!,
    magnitude,
    magnitude_phase,
    meanfinite,
    ncc,
    orientation,
    padarray,
    phase,
    prewitt,
    sad,
    sadn,
    sobel,
    ssd,
    ssdn,
    thin_edges,
    thin_edges_subpix,
    thin_edges_nonmaxsup,
    thin_edges_nonmaxsup_subpix,

    # phantoms
    shepp_logan

export # Deprecated exports
    ClipMin,
    ClipMax,
    ClipMinMax,
    ScaleInfo,
    climdefault,
    float32sc,
    float64sc,
    ntsc2rgb,
    rgb2gray,
    rgb2ntsc,
    scaleinfo,
    scaleminmax,
    scalesigned

import Base: scale, scale!  # delete when deprecations are removed
@deprecate scaleminmax  ScaleMinMax
@deprecate scaleminmax(img::AbstractArray, min::Real, max::Real)  ScaleMinMax(RGB24, img, min, max)
@deprecate float32sc    float32
@deprecate float64sc    float64
@deprecate uint8sc      ufixed8sc
@deprecate uint16sc(img)  ufixedsc(Ufixed16, img)
@deprecate ClipMin      ClampMin
@deprecate ClipMax      ClampMax
@deprecate ClipMinMax   ClampMinMax
@deprecate climdefault(img) zero(eltype(img)), one(eltype(img))
@deprecate ScaleMinMax{T<:Real}(img::AbstractArray{T}, mn, mx) ScaleMinMax(Ufixed8, img, mn, mx)
@deprecate ScaleMinMax{T<:ColorValue}(img::AbstractArray{T}, mn, mx) ScaleMinMax(RGB{Ufixed8}, img, mn, mx)
@deprecate ntsc2rgb(A::AbstractArray)  convert(Array{RGB}, A)
@deprecate ntsc2rgb(A::AbstractImage)  convert(Image{RGB}, A)
@deprecate rgb2ntsc(A::AbstractArray)  convert(Array{YIQ}, A)
@deprecate rgb2ntsc(A::AbstractImage)  convert(Image{YIQ}, A)
@deprecate scaleinfo    mapinfo
@deprecate scale(mapi::MapInfo, A) map(mapi, A)                # delete imports above when eliminated
@deprecate scale!(dest, mapi::MapInfo, A) map!(mapi, dest, A)  #   "
@deprecate rgb2gray(img::AbstractArray)  convert(Array{Gray}, img)
@deprecate copy(A::AbstractArray, B::AbstractArray) copyproperties(A, B)
@deprecate share(A::AbstractArray, B::AbstractArray) shareproperties(A, B)

const ScaleInfo = MapInfo  # can't deprecate types?

end


text:
include("totalreward.jl")
include("featuredgraph/fg.jl")
include("geometricflux/geometricflux.jl")


text:
"""
    GeomBMModel(startdate, startprice, interestrate, carryrate, volatility)

A model for `SingleStock`, following a geometric Brownian motion.
* `startdate`
* `startprice`: initial price at `startdate`
* `interestrate`: risk free rate of return.
* `carryrate`: the carry rate, i.e. the net return for holding the asset:
  - for stocks this is typically positive (i.e. dividends)
  - for commodities this is typically negative (i.e. cost-of-carry)
* `volatility`:

The `interestrate`, `carryrate` and `volatility` are all specified on a continously compounded, Actual/365 basis.

The price is assumed to follow the PDE:

```math
d S_t = \\kappa S_t dt + \\sigma S_t dW_t
```

or

```math
d \\log S_t = (\\kappa - \\sigma^2/2) dt + \\sigma dW_t
```

where ``W_t`` is a Wiener process, and `κ = interestrate - carryrate`.
"""
struct GeomBMModel{C,V} <: AbstractModel
    core::C
    volatility::V
end

GeomBMModel(startdate, startprice, yieldrate, carryrate, volatility) =
    GeomBMModel(CoreModel(startdate, startprice, yieldrate, carryrate),
                volatility)

numeraire(m::GeomBMModel) = numeraire(m.core)
startdate(m::GeomBMModel) = startdate(m.core)

yearfractionto(m::GeomBMModel, dt::Date) =
    yearfractionto(m.core,  dt)
forwardprice(m::GeomBMModel, s::SingleStock, dt::Date) =
    forwardprice(m.core, s, dt)

value(m::GeomBMModel, c::SingleStock) = value(m.core, c)
value(m::GeomBMModel, c::WhenAt{SingleStock}) = value(m.core, c)
value(m::GeomBMModel, c::WhenAt{Receive{T}}) where {T} = value(m.core,c)

function value(m::GeomBMModel, c::EuropeanCall{SingleStock,T}) where T
    dt = maturitydate(c)
    σ = m.volatility
    black(value(m, WhenAt(dt, SingleStock())),
          value(m, WhenAt(dt, Receive(strikeprice(c)))),
          σ*sqrt(yearfractionto(m, dt)))
end
function value(m::GeomBMModel, c::EuropeanPut{SingleStock,T}) where T
    dt = maturitydate(c)
    σ = m.volatility
    black(value(m, WhenAt(dt, Receive(strikeprice(c)))),
          value(m, WhenAt(dt, SingleStock())),
          σ*sqrt(yearfractionto(m, dt)))
end

# partial derivatives
function black_p1(S′,K′,σ′)
    x  = log(S′/K′)
    d₊ = x/σ′ + σ′/2
    normcdf(d₊)
end
function black_p2(S′,K′,σ′)
    x  = log(S′/K′)
    d₋ = x/σ′ - σ′/2
    -normcdf(d₋)
end
function black_p3(S′,K′,σ′)
    x  = log(S′/K′)
    d₊ = x/σ′ + σ′/2
    S′*normpdf(d₊)
end

function delta(m::GeomBMModel, c::EuropeanCall{SingleStock,T}) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    S = m.core.startprice
    K = strikeprice(c)
    σ = m.volatility
    iQ = discount(m.core.carrycurve, dt)
    iR = discount(m.core.yieldcurve, dt)
    iQ * black_p1(S*iQ, K*iR, σ*sqrt(t))
end
function delta(m::GeomBMModel, c::EuropeanPut{SingleStock,T}) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    S = m.core.startprice
    K = strikeprice(c)
    σ = m.volatility
    iQ = discount(m.core.carrycurve, dt)
    iR = discount(m.core.yieldcurve, dt)
    iQ * black_p2(K*iR, S*iQ, σ*sqrt(t))
end

function vega(m::GeomBMModel, c::EuropeanCall{SingleStock,T}) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    S = m.core.startprice
    K = strikeprice(c)
    σ = m.volatility
    sqrt(t) * black_p3(S*discount(m.core.carrycurve, dt), K*discount(m.core.yieldcurve, dt), σ*sqrt(t))
end
function vega(m::GeomBMModel, c::EuropeanPut{SingleStock,T}) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    S = m.core.startprice
    K = strikeprice(c)
    σ = m.volatility
    sqrt(t) * black_p3(K*discount(m.core.yieldcurve, dt), S*discount(m.core.carrycurve, dt), σ*sqrt(t))
end

ivol(m::GeomBMModel, c::Contract) = m.volatility

"""
    ivol(m::CoreModel, c::Contract, price)

Compute the Black-Scholes implied volatility of contract `c` at `price`, under the assumptions of model `m` (ignoring the volatility value of `m`).

See also: `fit`
"""
function ivol(m::CoreModel, c::EuropeanCall{SingleStock,T}, price) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    S = m.startprice
    K = strikeprice(c)
    black_ivol(S*discount(m.carrycurve, dt), K*discount(m.yieldcurve, dt), price)/sqrt(t)
end
function ivol(m::CoreModel, c::EuropeanPut{SingleStock,T}, price) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    S = m.startprice
    K = strikeprice(c)
    black_ivol(K*discount(m.yieldcurve, dt), S*discount(m.carrycurve, dt), price)/sqrt(t)
end

function ivol(m::CoreForwardModel, c::EuropeanCall{SingleStock,T}, price) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    F = forwardprice(m, SingleStock(), dt)
    K = strikeprice(c)
    D = discount(m.yieldcurve, dt)
    black_ivol(F*D, K*D, price)/sqrt(t)
end
function ivol(m::CoreForwardModel, c::EuropeanPut{SingleStock,T}, price) where T
    dt = maturitydate(c)
    t = yearfractionto(m, dt)
    F = forwardprice(m, SingleStock(), dt)
    K = strikeprice(c)
    D = discount(m.yieldcurve, dt)
    black_ivol(K*D, F*D, price)/sqrt(t)
end

"""
    fit(GeomBMModel, m::Union{CoreModel,CoreForwardModel}, c::Contract, price)

Fit a `GeomBMModel` using the implied volatility of `c` at `price`, using the parameters
of `m`.

See also: `ivol`
"""
fit(::Type{GeomBMModel}, mcore::AbstractCoreModel, c::Union{EuropeanCall{SingleStock,T},EuropeanPut{SingleStock,T}}, price) where {T} =
    GeomBMModel(mcore, ivol(mcore, c, price))


text:
println("Welcome to rate-plot.jl, where rates for CIRIS are plotted")
using DataFrames
using Gadfly

path = "/scratch/cns7ae/losalamos/"
outfile = "rates-plot.pdf"
xmax = 1e15
xmin = 1e13
speciesName = "O3"
speciesNum  = 7

plottheme = Gadfly.Theme(
                        highlight_width=0pt,
                        default_point_size=2pt,
                        line_width=2pt,
                        minor_label_font_size=13pt,
                        minor_label_color=colorant"black",
                        major_label_font_size=15pt,
                        major_label_color=colorant"black",
                        key_label_color=colorant"black",
                        key_title_color=colorant"black",
                        key_title_font_size=15pt,
                        key_label_font_size=13pt,
                        grid_color=colorant"black")

reactions_outpule_filename = "ozone_reactions.wsv"
reactionOutput = path*reactions_outpule_filename
println("The reactions output file is $reactionOutput")

speciesInput_filename     = "species.dat"
speciesInput = path*speciesInput_filename
println("The species input file is $speciesInput")

reactionsInput_filename   = "reactions.dat"
reactionsInput = path*reactionsInput_filename
println("The reactions input file is $reactionsInput")

#=
This function finds the index of something in a vector array in which it is
unique
=#
function findindex(p,itr)
  if p == "0"
    return 0
  end
  for n = 1:size(itr,1)
    if p == itr[n]
      return n
    end
  end
  return 666
end

#=
This functions Counts how many of something is in a vector array
=#
function howmany(thing,list)
  count = 0
  for n = 1:size(list,1)
    println("Thing is $thing and list is $(list[n])")
    if typeof(thing) == Int
      if thing == list[n]
        count =+ 1
      end
    else
      if stripAndAscii(thing) == stripAndAscii(list[n])
        count =+ 1
      end
    end
  end
  return count
end

function in_check(in_r,in_p,out_r,out_p)
  if size(union(in_r,out_r),1) == size(unique(in_r),1)
    test1 = true
  else
    test1 = false
  end

  if size(union(in_p,out_p),1) == size(unique(in_p),1)
    test2 = true
  else
    test2 = false
  end

  return test1*test2
end

#=
This function determines the class of a reactions
=#
function reactionclass(reactants,products,species)
  testReact = howmany(species,reactants)
#  println("testReact=$testReact")
  testProd = howmany(species,products)
  println("testProd=$testProd")
  if testReact > testProd
#    println("Destruction")
    return "Destruction"
  elseif testProd > testReact
#    println("Production")
    return "Production"
  else
#    println("Neutral")
    return "Neutral"
  end
end

#=
This function takes the species names and makes proper superscripts and subscripts
=#
function formatMolecule(specieslist)
  i = 1
  for species in specieslist
    species = replace(species,"+","+")
    species = replace(species,"-","-")
    species = replace(species,"*","*")
    for letter in species
      test = tryparse(Int64,"$letter")
      if ( isnull(test) == false )
        species = replace(species,"$letter","$letter")
      end
    end
    specieslist[i] = species
    i += 1
  end
  return specieslist
end

function stripAndAscii(string)
  return ascii(strip(string))
end

#=
Generate the dataframes containing the input and output dataframes
=#
println("Generating initial dataframes")

ozoneReactions = readtable(
    reactionOutput,
    header = true,
    names = [:R1,:R2,:P1,:P2,:P3,:Fluence],
    eltypes = [Int64,Int64,Int64,Int64,Int64,Float64]
                            )

println("Generated dataframe of reactions output")

println("Now generating species dataframe")

species = readtable(
    speciesInput,
    names = [:Name, :E_D],
    eltypes = [String,Float64],
    separator = ',',
    header = false,
    allowcomments = true,
    commentmark = '!'
)

println("Generated species input dataframe")

inreactions = readtable(
    reactionsInput,
    names = [:R1,:R2,:P1,:P2,:P3,:α, :β, :γ, :rtype, :source],
    separator = ',',
    eltypes = [String,String,String,String,String,Float64,Float64,Float64,Int64,String],
    header = false,
    allowcomments = true,
    ignorepadding = true,
    commentmark = '!'
)

println("Generated inreactions dataframe")

reactions = DataFrame()
reactions[:R1] = inreactions[:R1]
reactions[:R2] = inreactions[:R2]
reactions[:P1] = inreactions[:P1]
reactions[:P2] = inreactions[:P2]
reactions[:P3] = inreactions[:P3]

println("Generated reactions input dataframe")

# Format strings in each dataframe
reactions[:R1] = formatMolecule(reactions[:R1])
reactions[:R2] = formatMolecule(reactions[:R2])
reactions[:P1] = formatMolecule(reactions[:P1])
reactions[:P2] = formatMolecule(reactions[:P2])
reactions[:P3] = formatMolecule(reactions[:P3])
species[:Name] = formatMolecule(species[:Name])

#=
Add reactions Type column to reactions dataframe and determine which one
each Type is
=#
println("Adding new column to reactions dataframe")
reactions[:Type] = "Neutral"
for i in 1:size(reactions[:Type],1)
  reactants = [reactions[:R1][i],reactions[:R2][i]]
  products  = [reactions[:P1][i],reactions[:P2][i],reactions[:P3][i]]
  reactions[:Type][i] = reactionclass(reactants,products,speciesName)
end

println("Adding Count column to reactions dataframe")
reactions[:Count] = 0
fluenceLimit = 1.0E5
fluenceArr = zeros(0)
rateAnalytics = Array(Int64,(0,size(reactions[:R1],1)))
#=
Go through the reactions output, add any new reactions, and Count up rates
=#
for i in 1:size(ozoneReactions[:R1],1)
  if i%1000 == 0
    println("$i of $(size(ozoneReactions[:R1],1))")
  end
  outReacts = [ozoneReactions[:R1][i];ozoneReactions[:R2][i]]
  outProds  = [ozoneReactions[:P1][i];ozoneReactions[:P2][i];ozoneReactions[:P3][i]]
  in_network = false
  # Go through reactions for each line in the output reaction file
  for j in 1:size(reactions[:R1],1)
    in_reacts = [findindex(reactions[:R1][j],species[:Name]);
                 findindex(reactions[:R2][j],species[:Name])]
    inProds  = [findindex(reactions[:P1][j],species[:Name]);
                 findindex(reactions[:P2][j],species[:Name]);
                 findindex(reactions[:P3][j],species[:Name])]
    in_network = in_check(in_reacts,inProds,outReacts,outProds)
    if in_network == true
      reactions[:Count][j] += 1
      break
    end
  end
  # If the reaction isn't in the network, add it
  if in_network == false
    println("$outReacts, $outProds not in network")
    r1 = species[:Name][outReacts[1]]
    r2 = species[:Name][outReacts[2]]
    if outProds[1] == 0
      p1 = "0"
    else
      p1 = species[:Name][outProds[1]]
    end
    if outProds[2] == 0
      p2 = "0"
    else
      p2 = species[:Name][outProds[2]]
    end
    if outProds[3] == 0
      p3 = "0"
    else
      p3 = species[:Name][outProds[3]]
    end
    # Update the reactions dataframe
    push!(reactions,[r1;r2;p1;p2;p3;reactionclass(outReacts,outProds,speciesNum);0])
    # Expand the rateAnalytics array
    rateAnalytics = hcat(rateAnalytics,zeros(Int64,size(rateAnalytics,1)))
  end
  # Test to see if the fluence condition has been met, if so, save outputs and
  # re-initialize Counters
  if ozoneReactions[:Fluence][i] > fluenceLimit
    rateAnalytics = vcat(rateAnalytics,transpose(reactions[:Count]))
    reactions[:Count] = 0
    fluenceLimit += 0.35*fluenceLimit
    fluenceArr = vcat(fluenceArr,fluenceLimit)
#    println("New fluence limit = $fluenceLimit")
  end
end

# Get the indices of all production and destruction reactions
prodReacts = Array(Int64,0)
destReacts = Array(Int64,0)
neutReacts = Array(Int64,0)
for i in 1:size(reactions[:Type],1)
  if reactions[:Type][i] == "Production"
    push!(prodReacts,i)
  elseif reactions[:Type][i] == "Destruction"
    push!(destReacts,i)
  elseif reactions[:Type][i] == "Neutral"
    push!(neutReacts,i)
  end
end
println("prodReacts=$prodReacts")
println("destReacts=$destReacts")
println("neutReacts=$neutReacts")

println("Now initializing the percentRates array")
totProdArr = zeros(0)
totDestArr = zeros(0)
totArr     = zeros(0)
percentRates = zeros(Float64,(size(rateAnalytics,1),size(rateAnalytics,2)))
println("percentRates array now initialized")
for i in 1:size(rateAnalytics,1)
  totalProd = 0
  totalDest = 0
  totalNeut = 0
  # Go through row and add up reaction types
  for j in 1:size(rateAnalytics,2)
    if in(j,prodReacts)
      totalProd += rateAnalytics[i,j]
    elseif in(j,destReacts)
      totalDest += rateAnalytics[i,j]
    elseif in(j,neutReacts)
      totalNeut += rateAnalytics[i,j]
    end
  end
  # Go back through row and get percentage rates
  for j in 1:size(rateAnalytics,2)
    if in(j,prodReacts)
      if totalProd != 0
        percentRates[i,j] = Float64(rateAnalytics[i,j])/Float64(totalProd)
      end
    elseif in(j,destReacts)
      if totalDest != 0
        percentRates[i,j] = Float64(rateAnalytics[i,j])/Float64(totalDest)
      end
    elseif in(j,neutReacts)
      if totalNeut != 0
        percentRates[i,j] = Float64(rateAnalytics[i,j])/Float64(totalNeut)
      end
    end
  end
  totProdArr = vcat(totProdArr,totalProd)
  totDestArr = vcat(totDestArr,totalDest)
  totArr     = vcat(totArr,totalProd+totalDest)
end

proddf = DataFrame()
destdf = DataFrame()
neutdf = DataFrame()
println("Plot dataframes initialized")
for i in 1:size(percentRates,2)
  r1 = reactions[:R1][i]
  r2 = reactions[:R2][i]
  prods = [reactions[:P1][i]; reactions[:P2][i];reactions[:P3][i]]
  rname = "$r1 + $r2 -> "
  reactants = "$r1 + $r2"
  pcount = 0
  for q in 1:3
    if pcount == 0
      if prods[q] != "0"
        rname = rname*"$(prods[q]) "
        pcount += 1
      end
    else
      if prods[q] != "0"
        rname = rname*"+ $(prods[q])"
      end
    end
  end
  if in(i,prodReacts)
    proddf = vcat(proddf,DataFrame(
                                   Fluence=fluenceArr,
                                   PercentTotal=percentRates[:,i],
                                  #  label="$rname"
                                   label="$reactants"
                                   ))
  elseif in(i,destReacts)
    destdf = vcat(destdf,DataFrame(
                                   Fluence=fluenceArr,
                                   PercentTotal=percentRates[:,i],
                                  #  label="$rname"
                                   label="$reactants"
                                   ))
  elseif in(i,neutReacts)
    neutdf = vcat(neutdf,DataFrame(
                                   Fluence=fluenceArr,
                                   PercentTotal=percentRates[:,i],
                                  #  label="$rname"
                                   label="$reactants"
                                   ))
  end
end

println("Now generating plots")
p1 = plot(
          proddf,
#          proddf[xmax.>proddf[:Fluence].>xmin,:],
#          proddf,
          x="Fluence",
          y="PercentTotal",
          color="label",
#          Geom.smooth,
          Geom.line,
          Guide.colorkey("Reactants"),
          # Guide.title("Production Reactions"),
          Guide.xlabel("Fluence"),
          Guide.ylabel("Fraction of total"),
#          Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
#          Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#          Scale.y_log10
         )

push!(p1,plottheme)
println("Plotted p1")

p2 = plot(
#          destdf[xmax.>destdf[:Fluence].>xmin,:],
          destdf,
          x="Fluence",
          y="PercentTotal",
          color="label",
#          Geom.smooth,
          Geom.line,
          # Guide.title("Destruction Reactions"),
          Guide.xlabel("Fluence"),
          Guide.colorkey("Reactants"),
          Guide.ylabel("Fraction of total"),
#          Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
#          Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#          Scale.y_log10
         )

push!(p2,plottheme)
println("Plotted p2")

totProdPercent = zeros(0)
totDestPercent = zeros(0)

for i in 1:size(totProdArr,1)
  if totArr[i] != 0
    totProdPercent = vcat(totProdPercent,totProdArr[i]/totArr[i])
    totDestPercent = vcat(totDestPercent,totDestArr[i]/totArr[i])
  else
    totProdPercent = vcat(totProdPercent,0.0)
    totDestPercent = vcat(totDestPercent,0.0)
  end
end
totdf = vcat(
             DataFrame(
                       Fluence=fluenceArr,
                       PercentTotal=totProdPercent,
                       label="Total production"
             ),
             DataFrame(
                       Fluence=fluenceArr,
                       PercentTotal=totDestPercent,
                       label="Total destruction"
             )
)

p3 = plot(
          totdf,
#          totdf[xmax.>totdf[:Fluence].>xmin,:],
          x="Fluence",
          y="PercentTotal",
          color="label",
          Guide.colorkey("Type"),
#          Geom.smooth,
          Geom.line,
          # Guide.title("Relative Production and Destruction"),
          Guide.xlabel("Fluence"),
          Guide.ylabel("Fraction of total"),
          #Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
#          Scale.x_log10(minvalue=xmin,maxvalue=xmax)
         )
println("Plotted p3")

push!(p3,plottheme)
println("Generating final output pdf")
draw(PDF(path*"f4a.pdf",6inch,4.5inch),p1)
draw(PDF(path*"f4b.pdf",6inch,4.5inch),p2)
draw(PDF(path*"f4c.pdf",6inch,4.5inch),p3)
draw(PDF(path*outfile,6inch,9inch),vstack(p1,p2,p3))

println("Ending rate-plot.jl!")
quit()


text:
# Device type and auxiliary functions

export
    devcount,
    CuDevice, name, totalmem, attribute, capability,
    list_devices

function devcount()
    count_box = ptrbox(Cint)
    @cucall(:cuDeviceGetCount, (Ptr{Cint},), count_box)
    ptrunbox(count_box)
end

immutable CuDevice
    ordinal::Cint
    handle::Cint

    function CuDevice(i::Integer)
        ordinal = convert(Cint, i)
        handle_box = ptrbox(Cint)
        @cucall(:cuDeviceGet, (Ptr{Cint}, Cint), handle_box, ordinal)
        new(ordinal, ptrunbox(handle_box))
    end
end

function name(dev::CuDevice)
    const buflen = 256
    buf = Array(Cchar, buflen)
    @cucall(:cuDeviceGetName, (Ptr{Cchar}, Cint, Cint),
                              buf, buflen, dev.handle)
    return bytestring(pointer(buf))
end

function totalmem(dev::CuDevice)
    mem_box = ptrbox(Csize_t)
    @cucall(:cuDeviceTotalMem, (Ptr{Csize_t}, Cint), mem_box, dev.handle)
    return ptrunbox(mem_box)
end

function attribute(dev::CuDevice, attrcode::Integer)
    value_box = ptrbox(Csize_t)
    @cucall(:cuDeviceGetAttribute, (Ptr{Cint}, Cint, Cint),
                                   value_box, attrcode, dev.handle)
    return ptrunbox(value_box)
end

function capability(dev::CuDevice)
    major_box = ptrbox(Cint)
    minor_box = ptrbox(Cint)
    @cucall(:cuDeviceComputeCapability, (Ptr{Cint}, Ptr{Cint}, Cint),
                                        major_box, minor_box, dev.handle)
    return VersionNumber(ptrunbox(major_box), ptrunbox(minor_box))
end

function list_devices()
    cnt = devcount()
    if cnt == 0
        println("No CUDA-capable device found.")
        return
    end

    for i = 0:cnt-1
        dev = CuDevice(i)
        nam = name(dev)
        tmem = iround(totalmem(dev) / (1024^2))
        cap = capability(dev)

        println("device[$i]: $(nam), capability $(cap.major).$(cap.minor), total mem = $tmem MB")
    end
end


text:
using Dates

using Pumas, GlobalSensitivity, CairoMakie, PumasPlots

peg_inf_model = @model begin

  @param begin
    logKa   ∈  RealDomain()
    logKe   ∈  RealDomain()
    logVd   ∈  RealDomain()
    logn    ∈  RealDomain()
    logd    ∈  RealDomain()
    logc    ∈  RealDomain()
    logEC50 ∈  RealDomain()
    ω²Ka   ∈ RealDomain(lower=0.0)
    ω²Ke   ∈ RealDomain(lower=0.0)
    ω²Vd   ∈ RealDomain(lower=0.0)
    ω²n    ∈ RealDomain(lower=0.0)
    ω²δ    ∈ RealDomain(lower=0.0)
    ω²c    ∈ RealDomain(lower=0.0)
    ω²EC50 ∈ RealDomain(lower=0.0)
    σ²PK ∈ RealDomain(lower=0.0)
    σ²PD ∈ RealDomain(lower=0.0)
  end

  @random begin
    ηKa   ~ Normal(0.0, sqrt(ω²Ka))
    ηKe   ~ Normal(0.0, sqrt(ω²Ke))
    ηVd   ~ Normal(0.0, sqrt(ω²Vd))
    ηn    ~ Normal(0.0, sqrt(ω²n))
    ηδ    ~ Normal(0.0, sqrt(ω²δ))
    ηc    ~ Normal(0.0, sqrt(ω²c))
    ηEC50 ~ Normal(0.0, sqrt(ω²EC50))
  end

  @pre begin
    # constants
    p = 100.0
    d = 0.001
    e = 1e-7
    s = 20000.0
    log_Ka   = logKa   + ηKa
    log_Ke   = logKe   + ηKe
    log_Vd   = logVd   + ηVd
    log_n    = logn    + ηn
    log_d    = logd    + ηδ
    log_c    = logc    + ηc
    log_EC50 = logEC50 + ηEC50
  end

  @init begin
    T = exp(log_c + log_d)/(p*e)
    I = (s*e*p - d*exp(log_c + log_d))/(p*exp(log_d)*e)
    W = (s*e*p - d*exp(log_c + log_d))/(exp(log_c + log_d)*e)
  end

  @dynamics begin
    X' = -exp(log_Ka)*X
    A' = exp(log_Ka)*X - exp(log_Ke)*A
    T' = s - T*(e*W + d)
    I' = e*W*T - exp(log_d)*I
    W' = p/((A/exp(log_Vd)/exp(log_EC50) + 1e-100)^exp(log_n) + 1)*I - exp(log_c)*W
  end

  @derived begin
    conc   = @. A/exp(log_Vd)
    log10W = @. log10(max.(W, 1e-5))
    yPK ~ @. Normal(A/exp(log_Vd), sqrt(σ²PK))
    yPD ~ @. Normal(log10W, sqrt(σ²PD))
    nca1 := @nca conc
    cmax_pk = NCA.cmax(nca1)
    nca2 := @nca log10W
    cmax_pd = NCA.cmax(nca2)
  end

end

peg_inf_dr = DosageRegimen(180.0)
t = collect(0.0:1.0:28.0)
_pop = map(i -> Subject(id=i, observations=(yPK=[], yPD=[]), events=peg_inf_dr, time=t), 1:3)

param_PKPD = (
  logKa   = log(0.80),
  logKe   = log(0.15),
  logVd   = log(100.0),
  logn    = log(2.0),
  logd    = log(0.20),
  logc    = log(7.0),
  logEC50 = log(0.12),
  # random effects variance parameters, set to 0 to neutralize effect of random effects
  ω²Ka   = 0.0,
  ω²Ke   = 0.0,
  ω²Vd   = 0.0,
  ω²n    = 0.0,
  ω²δ    = 0.0,
  ω²c    = 0.0,
  ω²EC50 = 0.0,
  # variance parameter in proportional error model
  σ²PK = 0.04,
  σ²PD = 0.04)

simdata = simobs(peg_inf_model, _pop, param_PKPD, ensemblealg=EnsembleThreads())
sim_plot(peg_inf_model ,simdata, observations=[:yPD,:yPK], all = true)

p_range_low =  (logKa   = log(0.80)/2,
                logKe   = log(0.15)/2,
                logVd   = log(100.0)/2,
                logn    = log(2.0)/2,
                logd    = log(0.20)/2,
                logc    = log(7.0)/2,
                logEC50 = log(0.12)/2, )

p_range_high = (logKa   = log(0.80)*2,
                logKe   = log(0.15)*2,
                logVd   = log(100.0)*2,
                logn    = log(2.0)*2,
                logd    = log(0.20)*2,
                logc    = log(7.0)*2,
                logEC50 = log(0.12)*2, )

morris_ = Pumas.gsa(
  peg_inf_model,
  _pop,param_PKPD,
  GlobalSensitivity.Morris(),
  [:cmax_pk,:cmax_pd],
  p_range_low,
  p_range_high,
  ensemblealg=EnsembleThreads())

keys_ = [keys(p_range_low)...]
cmax_pk_meansstar = [morris_.means_star[1,:][key] for key in keys_]
cmax_pd_meansstar = [morris_.means_star[2,:][key] for key in keys_]

fig = Figure(resolution = (1200, 400))
plot_pk = barplot(fig[1,1], string.(keys_), cmax_pk_meansstar, axis = (xlabel ="Parameters", title="Cmax PK Morris means"))
plot_pd = barplot(fig[1,2], string.(keys_), cmax_pd_meansstar, axis = (xlabel ="Parameters", title="Cmax PD Morris means"))
display(fig)

sobol_ = Pumas.gsa(
  peg_inf_model,
  _pop,
  param_PKPD,
  GlobalSensitivity.Sobol(),
  [:cmax_pk,:cmax_pd],
  p_range_low,
  p_range_high,
  N=4000,
  ensemblealg=EnsembleThreads())

cmax_pk_s1 = [sobol_.first_order[1,:][key] for key in keys_]
cmax_pd_s1 = [sobol_.first_order[2,:][key] for key in keys_]
cmax_s1 = hcat(cmax_pk_s1, cmax_pd_s1)

fig = Figure(resolution = (1200, 400))
plot_s1 = heatmap(fig[1,1], 1:2, 1:7, cmax_s1, axis = (xticks = (1:2, ["Cmax_PK","Cmax_PD"]), yticks = (1:7, string.(keys_)), title="First Order Indices"), colormap = "darkrainbow")

cmax_pk_st = [sobol_.total_order[1,:][key] for key in keys_]
cmax_pd_st = [sobol_.total_order[2,:][key] for key in keys_]
cmax_st = hcat(cmax_pk_st, cmax_pd_st)

plot_st = heatmap(fig[1,2], 1:2, 1:7, cmax_st,  axis = (xticks = (1:2, ["Cmax_PK","Cmax_PD"]), yticks = (1:7, string.(keys_)), title="First Order Indices"), colormap = "darkrainbow")
display(fig)


text:
# ---
# title: 1406. Stone Game III
# id: problem1406
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/stone-game-iii/description/>
# hidden: true
# ---
#
# Alice and Bob continue their games with piles of stones. There are several
# stones  **arranged in a row** , and each stone has an associated value which
# is an integer given in the array `stoneValue`.
#
# Alice and Bob take turns, with **Alice** starting first. On each player's
# turn, that player can take **1, 2 or 3 stones**  from the **first** remaining
# stones in the row.
#
# The score of each player is the sum of values of the stones taken. The score
# of each player is **0**  initially.
#
# The objective of the game is to end with the highest score, and the winner is
# the player with the highest score and there could be a tie. The game continues
# until all the stones have been taken.
#
# Assume Alice and Bob  **play optimally**.
#
# Return _" Alice"_ if Alice will win, _" Bob"_ if Bob will win or _" Tie"_ if
# they end the game with the same score.
#
#
#
# **Example 1:**
#
#
#
#     Input: values = [1,2,3,7]
#     Output: "Bob"
#     Explanation: Alice will always lose. Her best move will be to take three piles and the score become 6. Now the score of Bob is 7 and Bob wins.
#
#
# **Example 2:**
#
#
#
#     Input: values = [1,2,3,-9]
#     Output: "Alice"
#     Explanation: Alice must choose all the three piles at the first move to win and leave Bob with negative score.
#     If Alice chooses one pile her score will be 1 and the next move Bob's score becomes 5. The next move Alice will take the pile with value = -9 and lose.
#     If Alice chooses two piles her score will be 3 and the next move Bob's score becomes 3. The next move Alice will take the pile with value = -9 and also lose.
#     Remember that both play optimally so here Alice will choose the scenario that makes her win.
#
#
# **Example 3:**
#
#
#
#     Input: values = [1,2,3,6]
#     Output: "Tie"
#     Explanation: Alice cannot win this game. She can end the game in a draw if she decided to choose all the first three piles, otherwise she will lose.
#
#
# **Example 4:**
#
#
#
#     Input: values = [1,2,3,-1,-2,-3,7]
#     Output: "Alice"
#
#
# **Example 5:**
#
#
#
#     Input: values = [-1,-2,-3]
#     Output: "Tie"
#
#
#
#
# **Constraints:**
#
#   * `1 <= values.length <= 50000`
#   * `-1000 <= values[i] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using CompScienceMeshes, BEAST

function blkdiagm(blks...)
    m = sum(size(b,1) for b in blks)
    n = sum(size(b,2) for b in blks)
    A = zeros(eltype(first(blks)), m, n)
    o1 = 1
    o2 = 1
    for b in blks
        m1 = size(b,1)
        n1 = size(b,2)
        A[o1:o1+m1-1,o2:o2+n1-1] .= b
        o1 += m1
        o2 += n1
    end
    A
end

width, height, h = 1.0, 0.5, 0.05
G1 = meshrectangle(width, height, h)
G2 = CompScienceMeshes.rotate(G1, 0.5π * x̂)
G3 = CompScienceMeshes.rotate(G1, 1.0π * x̂)

G12 = weld(G1,-G2)
G23 = weld(G2,-G3)
G31 = weld(G3,-G1)

X12 = lagrangec0d1(G12)
X23 = lagrangec0d1(G23)
X31 = lagrangec0d1(G31)

Y12 = duallagrangecxd0(G12)
Y23 = duallagrangecxd0(G23)
Y31 = duallagrangecxd0(G31)

κ = 1.0
SL = Helmholtz3D.singlelayer(wavenumber=κ)
HS = Helmholtz3D.hypersingular(gamma=κ*im)
N = Identity()

X = X12 × X23 × X31
Y = Y12 × Y23 × Y31

E = Helmholtz3D.planewave(direction=ẑ, wavenumber=κ)
e = strace(E, G12)
ex = assemble(e, X)

Sxx = assemble(SL, X, X)
N1 = assemble(Identity(), X12, Y12)
N2 = assemble(Identity(), X23, Y23)
N3 = assemble(Identity(), X31, Y31)
Nxy = blkdiagm(N1,N2,N3)
Syy = assemble(HS, Y, Y)

Dyx = inv(Nxy)
Q = transpose(Dyx) * Syy * Dyx * Sxx
R = transpose(Dyx) * Syy * Dyx * ex

u1, ch1 = solve(BEAST.GMRESSolver(Sxx),ex)
u2, ch2 = solve(BEAST.GMRESSolver(Q),R)

@show ch1.iters
@show ch2.iters

using LinearAlgebra
cond(Sxx)
cond(Q)

w1 = eigvals(Sxx)
w2 = eigvals(Q)

using Plots
plot()
scatter!(w1)
scatter!(w2)


text:
using CUDA
using CUDA.CUSPARSE
using Printf
using FiniteDiff
using ForwardDiff
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using ExaPF
import ExaPF: PowerSystem, AutoDiff

const PS = PowerSystem

@testset "Polar formulation" begin
    datafile = joinpath(dirname(@__FILE__), "..", "data", "case9.m")
    tolerance = 1e-8
    pf = PS.PowerNetwork(datafile, 1)

    if has_cuda_gpu()
        ITERATORS = zip([CPU(), CUDADevice()], [Array, CuArray])
    else
        ITERATORS = zip([CPU()], [Array])
    end

    @testset "Initiate polar formulation on device $device" for (device, M) in ITERATORS
        polar = PolarForm(pf, device)
        # Test printing
        println(devnull, polar)

        b = bounds(polar, State())
        b = bounds(polar, Control())

        # Test getters
        nᵤ = ExaPF.get(polar, NumberOfControl())
        nₓ = ExaPF.get(polar, NumberOfState())

        # Get initial position
        x0 = ExaPF.initial(polar, State())
        u0 = ExaPF.initial(polar, Control())
        p = ExaPF.initial(polar, Parameters())

        @test length(u0) == nᵤ
        @test length(x0) == nₓ
        for v in [x0, u0, p]
            @test isa(v, M)
        end

        cache = ExaPF.get(polar, ExaPF.PhysicalState())
        @testset "NetworkState cache" begin
            @test isa(cache.vmag, M)
            ExaPF.transfer!(polar, cache, x0, u0, p)
        end
        @testset "Polar model API" begin
            xₖ = copy(x0)
            # Init AD factory
            jx, ju, ∂obj = ExaPF.init_autodiff_factory(polar, cache)

            # Test powerflow with cache signature
            conv = powerflow(polar, jx, cache, tol=tolerance)
            ExaPF.get!(polar, State(), xₖ, cache)
            # Refresh power of generators in cache
            for Power in [PS.ActivePower, PS.ReactivePower]
                ExaPF.refresh!(polar, PS.Generator(), Power(), cache)
            end

            # Bounds on state and control
            u_min, u_max = ExaPF.bounds(polar, Control())
            x_min, x_max = ExaPF.bounds(polar, State())

            @test isequal(u_min, [0.9, 0.1, 0.1, 0.9, 0.9])
            @test isequal(u_max, [1.1, 3.0, 2.7, 1.1, 1.1])
            @test isequal(x_min, [-Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9])
            @test isequal(x_max, [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1])

            # Test callbacks
            ## Power Balance
            ExaPF.power_balance!(polar, cache)
            g = cache.balance
            @test isa(g, M)
            # As we run powerflow before, the balance should be below tolerance
            @test norm(g, Inf) < tolerance
            ## Cost Production
            c2 = ExaPF.cost_production(polar, cache.pg)
            @test isa(c2, Real)
            ## Inequality constraint
            for cons in [ExaPF.state_constraint, ExaPF.power_constraints]
                m = ExaPF.size_constraint(polar, cons)
                @test isa(m, Int)
                g = M{Float64, 1}(undef, m) # TODO: this signature is not great
                fill!(g, 0)
                cons(polar, g, cache)

                g_min, g_max = ExaPF.bounds(polar, cons)
                @test length(g_min) == m
                @test length(g_max) == m
                # Are we on the correct device?
                @test isa(g_min, M)
                @test isa(g_max, M)
                # Test constraints are consistent
                @test isless(g_min, g_max)
            end
        end
    end
end


text:
using Revise: convert
using Revise

includet("../src/IPhreeqc.jl")

iph = IPhreeqc.CreateIPhreeqc()

IPhreeqc.LoadDatabase(iph, "../database/phreeqc.dat")
# iph = convert(Int, iph)

IPhreeqc.AccumulateLine(iph, "SOLUTION 1")
IPhreeqc.AccumulateLine(iph, "Na  0.1")
IPhreeqc.AccumulateLine(iph, "Cl  0.1")
IPhreeqc.AccumulateLine(iph, "pH  9.0  charge")
IPhreeqc.AccumulateLine(iph, "END")

IPhreeqc.SetOutputStringOn(iph, 1)
IPhreeqc.RunAccumulated(iph)

res = IPhreeqc.GetOutputString(iph)

print(res)


text:
module process

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh