text:
# Copyright 2017 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

export MultinomialKernel

mutable struct MultinomialKernel <: NaiveBayesKernel
    acc_given_y::Matrix{Float64}
    smoothing::Float64
end

function MultinomialKernel(X::AbstractVector{ItemType}, y::AbstractVector{Int}, nlabels::Int; smoothing::Float64=0.0) where ItemType
    C = zeros(Float64, length(X[1]), nlabels)

    @inbounds for i in 1:length(X)
        label = y[i]
        for (j, x) in enumerate(X[i])
            C[j, label] += x
        end
    end

    MultinomialKernel(C, smoothing)
end

function kernel_prob(nbc::NaiveBayesClassifier, kernel::MultinomialKernel, x::AbstractVector{Float64})::Vector{Float64}
    n = length(nbc.le.labels)
    scores = zeros(Float64, n)
    for i in 1:n
        pxy = 1.0
        py = nbc.probs[i]
        for j in 1:length(x)
            den = kernel.acc_given_y[j, i] + kernel.smoothing * n
            pxy *= (x[j] + kernel.smoothing) / den
        end
        scores[i] = py * pxy
    end

    scores
end


text:
module MarkovChains

include("diffusionProcesses.jl")
include("CTMC.jl")

using LinearAlgebra
using SparseArrays
using Distributions

import Base:
	show

# Export types
export
	ContinuousTimeMarkovChain,
	ItoDiffusionProcess

# Export methods
export
	stationary_distribution,
	random_sample

"""
A Julia package for simulating Markov Chains.

"""

end

ext: jl


text:
"""
    ESS

Elliptical slice sampling algorithm.

# Examples
```jldoctest; setup = :(Random.seed!(1))
julia> @model gdemo(x) = begin
           m ~ Normal()
           x ~ Normal(m, 0.5)
       end
gdemo (generic function with 2 methods)

julia> sample(gdemo(1.0), ESS(), 1_000) |> mean
Mean

│ Row │ parameters │ mean     │
│     │ Symbol     │ Float64  │
├─────┼────────────┼──────────┤
│ 1   │ m          │ 0.824853 │
```
"""
struct ESS{space} <: InferenceAlgorithm end

ESS() = ESS{()}()
ESS(space::Symbol) = ESS{(space,)}()

# always accept in the first step
function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:ESS},
    vi::AbstractVarInfo;
    kwargs...
)
    # Sanity check
    vns = _getvns(vi, spl)
    length(vns) == 1 ||
        error("[ESS] does only support one variable ($(length(vns)) variables specified)")
    for vn in vns[1]
        dist = getdist(vi, vn)
        EllipticalSliceSampling.isgaussian(typeof(dist)) ||
            error("[ESS] only supports Gaussian prior distributions")
    end

    return Transition(vi), vi
end

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:ESS},
    vi::AbstractVarInfo;
    kwargs...
)
    # obtain previous sample
    f = vi[spl]

    # define previous sampler state
    # (do not use cache to avoid in-place sampling from prior)
    oldstate = EllipticalSliceSampling.ESSState(f, getlogp(vi), nothing)

    # compute next state
    sample, state = AbstractMCMC.step(
        rng,
        EllipticalSliceSampling.ESSModel(
            ESSPrior(model, spl, vi), ESSLogLikelihood(model, spl, vi),
        ),
        EllipticalSliceSampling.ESS(),
        oldstate,
    )

    # update sample and log-likelihood
    vi[spl] = sample
    setlogp!(vi, state.loglikelihood)

    return Transition(vi), vi
end

# Prior distribution of considered random variable
struct ESSPrior{M<:Model,S<:Sampler{<:ESS},V<:AbstractVarInfo,T}
    model::M
    sampler::S
    varinfo::V
    μ::T

    function ESSPrior{M,S,V}(model::M, sampler::S, varinfo::V) where {
        M<:Model,S<:Sampler{<:ESS},V<:AbstractVarInfo
    }
        vns = _getvns(varinfo, sampler)
        μ = mapreduce(vcat, vns[1]) do vn
            dist = getdist(varinfo, vn)
            EllipticalSliceSampling.isgaussian(typeof(dist)) ||
                error("[ESS] only supports Gaussian prior distributions")
            vectorize(dist, mean(dist))
        end
        return new{M,S,V,typeof(μ)}(model, sampler, varinfo, μ)
    end
end

function ESSPrior(model::Model, sampler::Sampler{<:ESS}, varinfo::AbstractVarInfo)
    return ESSPrior{typeof(model),typeof(sampler),typeof(varinfo)}(
        model, sampler, varinfo,
    )
end

# Ensure that the prior is a Gaussian distribution (checked in the constructor)
EllipticalSliceSampling.isgaussian(::Type{<:ESSPrior}) = true

# Only define out-of-place sampling
function Base.rand(rng::Random.AbstractRNG, p::ESSPrior)
    sampler = p.sampler
    varinfo = p.varinfo
    vns = _getvns(varinfo, sampler)
    set_flag!(varinfo, vns[1][1], "del")
    p.model(rng, varinfo, sampler)
    return varinfo[sampler]
end

# Mean of prior distribution
Distributions.mean(p::ESSPrior) = p.μ

# Evaluate log-likelihood of proposals
struct ESSLogLikelihood{M<:Model,S<:Sampler{<:ESS},V<:AbstractVarInfo}
    model::M
    sampler::S
    varinfo::V
end

function (ℓ::ESSLogLikelihood)(f)
    sampler = ℓ.sampler
    varinfo = ℓ.varinfo
    varinfo[sampler] = f
    ℓ.model(varinfo, sampler)
    return getlogp(varinfo)
end

function DynamicPPL.tilde(rng, ctx::DefaultContext, sampler::Sampler{<:ESS}, right, vn::VarName, inds, vi)
    if inspace(vn, sampler)
        return DynamicPPL.tilde(rng, LikelihoodContext(), SampleFromPrior(), right, vn, inds, vi)
    else
        return DynamicPPL.tilde(rng, ctx, SampleFromPrior(), right, vn, inds, vi)
    end
end

function DynamicPPL.tilde(ctx::DefaultContext, sampler::Sampler{<:ESS}, right, left, vi)
    return DynamicPPL.tilde(ctx, SampleFromPrior(), right, left, vi)
end

function DynamicPPL.dot_tilde(rng, ctx::DefaultContext, sampler::Sampler{<:ESS}, right, left, vn::VarName, inds, vi)
    if inspace(vn, sampler)
        return DynamicPPL.dot_tilde(rng, LikelihoodContext(), SampleFromPrior(), right, left, vn, inds, vi)
    else
        return DynamicPPL.dot_tilde(rng, ctx, SampleFromPrior(), right, left, vn, inds, vi)
    end
end

function DynamicPPL.dot_tilde(rng, ctx::DefaultContext, sampler::Sampler{<:ESS}, right, left, vi)
    return DynamicPPL.dot_tilde(rng, ctx, SampleFromPrior(), right, left, vi)
end


text:
using Reexport

@reexport module Simplices
	import Distributions: Uniform

	# Abstract type and functions that should work on all composite types
	include("AbstractSimplex.jl")

	# Simplices whose vertices are represented by vectors
	include("composite_types/Simplex.jl")
	include("composite_types/MutableSimplex.jl")

	# Simplices whose vertices are represented by SVectors or MVectors
	include("composite_types/SSimplex.jl")
	include("composite_types/MutableSSimplex.jl")

	# Various functions
	include("point_contained_in_simplex.jl")

	# Plot recipes
	include("plotting/plot_recipes.jl")

end # module


text:
# Based on this guide: https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf

using Atomistic
using AtomsBase
using InteratomicPotentials
using NBodySimulator
using Plots
using Unitful
using UnitfulAtomic

N = 864
element = :Ar
box_size = 3.47786u"nm"
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = DynamicSystem(initial_bodies, box_size)
eq_simulator = NBSimulator(
    Δt = Δt,
    steps = 2000,
    t₀ = 0.0u"s",
    thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
)
potential = LennardJones(
    austrip(1.657e-21u"J"),
    austrip(0.34u"nm")
)
eq_result = @time simulate(initial_system, eq_simulator, potential)

temp = @time plot_temperature(eq_result, 10)
energy = @time plot_energy(eq_result, 10)

prod_simulator = NBSimulator(
    Δt = Δt,
    steps = 5000,
    t₀ = eq_simulator.steps * Δt
)
prod_result = @time simulate(get_system(eq_result), prod_simulator, potential)

prod_stride = 10

display(@time plot_temperature!(temp, prod_result, 10))
display(@time plot_energy!(energy, prod_result, 10))

# rdf = @time plot_rdf(prod_result, potential.σ, 0.05)
# display(rdf)
# savefig(rdf, "artifacts/argon_lj_rdf.svg")

;

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
using PrinciplesCTImaging
using Test

include("lineintegralsandprojections.jl")
include("parallelreconstruction.jl")


text:
using HDF5

#file = ".data/PTBLM/ptb.test.txt"
#x = read(file)
#h5write("ptblm.h5", "test", x)
file = ".data/PTBLM/ptb.train.txt"
#x = read(file)
#h5write("ptblm.h5", "train", x)
strs = String[]
for line in open(readlines,file)
    elems = collect(split(chomp(line)))
    append!(strs, elems)
end
h5write("x.h5", "s", strs)


text:
using DrWatson

include(srcdir("datasets.jl"))
include(srcdir("models.jl"))
include(srcdir("utilities_train.jl"))
include(srcdir("utilities_evaluation.jl"))

DrWatson.datadir(args...) = joinpath("/disk/macha/data_deeptoppush", args...)

# -------------------------------------------------------------------------------
# Auxiliary functions
# -------------------------------------------------------------------------------
getdim(A::AbstractArray, d::Integer, i) = getindex(A, Base._setindex(i, d, axes(A)...)...)

function compute_scores(model, x; chunksize = 10000, device = identity)
    x_obs = ndims(x)
    n = size(x, x_obs)
    scores = zeros(eltype(x), n)

    for inds in partition(1:n, chunksize)
        scores[inds] .= model(device(getdim(x, x_obs, inds)))[:]
    end
    return scores
end

function dict_list_simple(d::Dict)
    ls = map(values(d)) do val
        typeof(val) <: AbstractVector ? length(val) : 1
    end
    if length(unique(ls)) > 2
        @error "not supported"
    else
        return map(1:maximum(ls)) do k
            Dict(key => typeof(val) <: AbstractVector ? val[k] : val for (key, val) in d)
        end
    end
end

logrange(x1, x2; kwargs...) = exp10.(range(log10(x1), log10(x2); kwargs...))

# -------------------------------------------------------------------------------
# Names....
# -------------------------------------------------------------------------------
allowedtypes(args...) = (Real, String, Symbol, DataType, Function, args...)

function dataset_savename(dataset_settings; digits = 4)
    return savename(dataset_settings; allowedtypes = allowedtypes(), digits = digits)
end

function train_savename(train_settings_in; digits = 4)
    train_settings = deepcopy(train_settings_in)
    delete!(train_settings, :epochlength)
    delete!(train_settings, :iters)
    return savename(train_settings; allowedtypes = allowedtypes(), digits = digits)
end

function model_savename(model_settings_in; digits = 4)
    model_settings = deepcopy(model_settings_in)
    delete!(model_settings, :seed)
    return savename(model_settings; allowedtypes = allowedtypes(), digits = digits)
end

function modeldir(
    dataset_settings,
    train_settings,
    model_settings;
    digits = 4,
    agg = datadir,
)

    datdir = dataset_savename(dataset_settings; digits = digits)
    traindir = train_savename(train_settings; digits = digits)
    moddir = model_savename(model_settings; digits = digits)

    return agg("models", datdir, traindir, moddir)
end

simulation_name(epoch) = string("model_epoch=", epoch, ".bson")
extract_model_type(d::Dict) = d[:model_settings][:type]
model_name(d::Dict) = model_name(d[:model_settings][:type], d[:model_settings])

function model_name(T::Type{<:Model}, d::Dict)
    vals = [string(d[key]) for key in model_args(T)]
    return string(T, "(", join(vals, ", "), ")")
end

model_args(::Type{<:Model}) = [:arg]
model_args(::Type{BaseLine}) = []
model_args(::Type{DeepTopPush}) = [:arg, :buffer]


text:
module TestExamplesDGMulti3D

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "dgmulti_3d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "DGMulti 3D" begin
  # 3d tet/hex tests
  @trixi_testset "elixir_euler_weakform.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
      l2 = [0.0010029534292051608, 0.0011682205957721673, 0.001072975385793516, 0.000997247778892257, 0.0039364354651358294] ./ sqrt(8),
      linf = [0.003660737033303718, 0.005625620600749226, 0.0030566354814669516, 0.0041580358824311325, 0.019326660236036464]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
      l2 = [0.014932088450136542, 0.017080219613061526, 0.016589517840793006, 0.015905000907070196, 0.03903416208587798] ./ sqrt(8),
      linf = [0.06856547797256729, 0.08225664880340489, 0.06925055630951782, 0.06913016119820181, 0.19161418499621874]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Hexahedral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      element_type = Hex(),
      # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
      l2 = [0.00030580190715769566, 0.00040146357607439464, 0.00040146357607564597, 0.000401463576075708, 0.0015749412434154315] ./ sqrt(8),
      linf = [0.00036910287847780054, 0.00042659774184228283, 0.0004265977427213574, 0.00042659774250686233, 0.00143803344597071]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
      l2 = [0.0010317074322517949, 0.0012277090547035293, 0.0011273991123913515, 0.0010418496196130177, 0.004058878478404962] ./ sqrt(8),
      linf = [0.003227752881827861, 0.005620317864620361, 0.0030514833972379307, 0.003987027618439498, 0.019282224709831652]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl (Hexahedral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      element_type = Hex(),
      # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
      l2 = [0.00034230612468547436, 0.00044397204714598747, 0.0004439720471461567, 0.0004439720471464591, 0.0016639410646990126] ./ sqrt(8),
      linf = [0.0003674374460325147, 0.0004253921341716982, 0.0004253921340786615, 0.0004253921340831024, 0.0014333414071048267]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl (Hexahedral elements, SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      element_type = Hex(),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
      l2 = [0.001712443468716032, 0.002491315550718859, 0.0024913155507195303, 0.002491315550720031, 0.008585818982343299] ./ sqrt(8),
      linf = [0.003810078279323559, 0.004998778644230928, 0.004998778643986235, 0.0049987786444081195, 0.016455044373650196]
    )
  end

  @trixi_testset "elixir_euler_taylor_green_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_taylor_green_vortex.jl"),
      polydeg = 3, tspan = (0.0, 1.0), cells_per_dimension = (2, 2, 2),
      l2   = [0.0003612827827560599, 0.06219350883951729, 0.062193508839503864, 0.08121963221634831, 0.07082703570808184],
      linf = [0.0007893509649821162, 0.1481953939988877, 0.14819539399791176, 0.14847291108358926, 0.21313533492212855]
      )
  end

  @trixi_testset "elixir_euler_taylor_green_vortex.jl (GaussSBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_taylor_green_vortex.jl"),
      polydeg = 3, approximation_type = GaussSBP(), tspan = (0.0, 1.0), cells_per_dimension = (2, 2, 2),
      l2 = [0.00036128278275524326, 0.062193508839511434, 0.06219350883949677, 0.08121963221635205, 0.07082703570765223],
      linf = [0.000789350964946367, 0.14819539399525805, 0.14819539399590542, 0.14847291107658706, 0.21313533492059378]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl (FD SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      element_type = Hex(),
      cells_per_dimension = (2, 2, 2),
      approximation_type = derivative_operator(
        SummationByPartsOperators.MattssonNordström2004(),
        derivative_order=1, accuracy_order=2,
        xmin=0.0, xmax=1.0, N=8),
      l2 = [0.0024092707138829925, 0.003318758964118284, 0.0033187589641182386, 0.003318758964118252, 0.012689348410504253],
      linf = [0.006118565824207778, 0.008486456080185167, 0.008486456080180282, 0.008486456080185611, 0.035113544599208346]
    )
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl (FD SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      element_type = Hex(),
      cells_per_dimension = (2, 2, 2),
      approximation_type = derivative_operator(
        SummationByPartsOperators.MattssonNordström2004(),
        derivative_order=1, accuracy_order=2,
        xmin=0.0, xmax=1.0, N=8),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      l2 = [0.0034543609010604407, 0.004944363692625066, 0.0049443636926250435, 0.004944363692625037, 0.01788695279620914],
      linf = [0.013861851418853988, 0.02126572106620328, 0.021265721066209053, 0.021265721066210386, 0.0771455289446683]
    )
  end

  @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
      l2 = [7.561896970347222e-5, 6.884047859374506e-5, 6.884047859377477e-5, 6.884047859375886e-5, 0.0002011072746170835],
      linf = [0.00013375200202747628, 0.00011571467800441937, 0.0001157146780030871, 0.00011571467800219892, 0.00034460823085824543]
    )
  end

end

# Clean up afterwards: delete Trixi output directory
@test_nowarn isdir(outdir) && rm(outdir, recursive=true)

end # module


text:
# Autogenerated wrapper script for CGAL_jll for any
export GGAL_config_cmake

using boost_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("CGAL")
JLLWrappers.@declare_file_product(GGAL_config_cmake)
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_file_product(
        GGAL_config_cmake,
        "lib/cmake/CGAL/CGALConfig.cmake",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# License is MIT: LICENSE.md

using ModuleInterfaceTools, Format

@api test StrLiterals, StrFormat

@testset "C Formatting" begin
    @testset "int" begin
        @test f"\%d(typemax(Int64))" == "9223372036854775807"
        @static if VERSION < v"1.5"
            @test f"\%a(typemax(Int64))" == "0x7.fffffffffffffffp+60"
            @test f"\%A(typemax(Int64))" == "0X7.FFFFFFFFFFFFFFFP+60"
        else
            @test f"\%a(typemax(Int64))" == "0x1p+63"
            @test f"\%A(typemax(Int64))" == "0X1P+63"
        end
    end
    @testset "printing an int value" begin
        for num in (UInt16(42), UInt32(42), UInt64(42), UInt128(42),
                    Int16(42), Int32(42), Int64(42), Int128(42))
            @test f"\%i(num)"       == "42"
            @test f"\%u(num)"       == "42"
            @test f"Test: \%i(num)" == "Test: 42"
            @test f"\%#x(num)"      == "0x2a"
            @test f"\%#o(num)"      == "052"
            @test f"\%x(num)"       == "2a"
            @test f"\%X(num)"       == "2A"
            @test f"\% i(num)"      == " 42"
            @test f"\%+i(num)"      == "+42"
            @test f"\%4i(num)"      == "  42"
            @test f"\%-4i(num)"     == "42  "
            @static if VERSION < v"1.5"
                @test f"\%a(num)"       == "0x2.ap+4"
                @test f"\%A(num)"       == "0X2.AP+4"
                @test f"\%20a(num)"     == "            0x2.ap+4"
                @test f"\%-20a(num)"    == "0x2.ap+4            "
            else
                @test f"\%a(num)"       == "0x1.5p+5"
                @test f"\%A(num)"       == "0X1.5P+5"
                @test f"\%20a(num)"     == "            0x1.5p+5"
                @test f"\%-20a(num)"    == "0x1.5p+5            "
            end
            @test f"\%f(num)"       == "42.000000"
            @test f"\%g(num)"       == "42"
        end
    end
    @testset "pointers" begin
        @static if Sys.WORD_SIZE == 64
            @test f"\%20p(0)"  == "  0x0000000000000000"
            @test f"\%-20p(0)" == "0x0000000000000000  "
        elseif Sys.WORD_SIZE == 32
            @test f"\%20p(0)"  == "          0x00000000"
            @test f"\%-20p(0)" == "0x00000000          "
        else
            @test false
        end
    end
    @testset "float / BigFloat" begin
        for num in (1.2345, big"1.2345")
            @test f"\%7.2f(num)"  == "   1.23"
            @test f"\%-7.2f(num)" == "1.23   "
            @test f"\%07.2f(num)" == "0001.23"
            @test f"\%.0f(num)"   == "1"
            @test f"\%#.0f(num)"  == "1."
            @test f"\%.4e(num)"   == "1.2345e+00"
            @test f"\%.4E(num)"   == "1.2345E+00"
            @test f"\%.2a(num)"   == "0x1.3cp+0"
            @test f"\%.2A(num)"   == "0X1.3CP+0"
        end
    end
    @testset "Inf / NaN handling" begin
        bf = big"Inf"
        bn = big"NaN"
        @test f"\%f(Inf)" == "Inf"
        @test f"\%f(NaN)" == "NaN"
        @test f"\%f(bf)"  == "Inf"
        @test f"\%f(bn)"  == "NaN"
    end
    @testset "scientific notation" begin
        fv = 3e142
        bv = big"3e142"
        @test f"\%.0e(fv)"  == "3e+142"
        @test f"\%#.0e(fv)" == "3.e+142"
        @test f"\%.0e(bv)"  == "3e+142"
        @test f"\%#.0e(bv)" == "3.e+142"
        bv = big"3e1042"
        @test f"\%.0e(bv)"   == "3e+1042"
        @test f"\%e(3e42)"   == "3.000000e+42"
        @test f"\%E(3e42)"   == "3.000000E+42"
        @test f"\%e(3e-42)"  == "3.000000e-42"
        @test f"\%E(3e-42)"  == "3.000000E-42"
        @test f"\%a(3e4)"    == "0x1.d4cp+14"
        @test f"\%A(3e4)"    == "0X1.D4CP+14"
        @test f"\%.4a(3e-4)" == "0x1.3a93p-12"
        @test f"\%.4A(3e-4)" == "0X1.3A93P-12"
    end
    @testset "%g" begin
        for (val, res) in (
            (12345678., "1.23457e+07"),
            (1234567.8, "1.23457e+06"),
            (123456.78, "123457"),
            (12345.678, "12345.7"),
            (12340000.0, "1.234e+07"),
            (big"12345678.", "1.23457e+07"),
            (big"1234567.8", "1.23457e+06"),
            (big"123456.78", "123457"),
            (big"12345.678", "12345.7"))
            @test f"\%.6g(val)" == res
        end
        for num in (123.4, big"123.4")
            @test f"\%10.5g(num)"   == "     123.4"
            @test f"\%+10.5g(num)"  == "    +123.4"
            @test f"\% 10.5g(num)"  == "     123.4"
            @test f"\%#10.5g(num)"  == "    123.40"
            @test f"\%-10.5g(num)"  == "123.4     "
            @test f"\%-+10.5g(num)" == "+123.4    "
            @test f"\%010.5g(num)"  == "00000123.4"
        end
        bv = big"-123.4"
        b2 = big"12340000.0"
        @test f"\%10.5g(-123.4)"    == "    -123.4"
        @test f"\%010.5g(-123.4)"   == "-0000123.4"
        @test f"\%.6g(12340000.0)"  == "1.234e+07"
        @test f"\%#.6g(12340000.0)" == "1.23400e+07"
        @test f"\%10.5g(bv)"        == "    -123.4"
        @test f"\%010.5g(bv)"       == "-0000123.4"
        @test f"\%.6g(b2)"          == "1.234e+07"
        @test f"\%#.6g(b2)"         == "1.23400e+07"
        @test f"\%.5g(42)"          == "42"
        @test f"\%#.2g(42)"         == "42."
        @test f"\%#.5g(42)"         == "42.000"
    end

    @testset "hex float" begin
        bv = big"1.5"
        @test f"\%a(1.5)"    == "0x1.8p+0"
        @test f"\%a(1.5f0)"  == "0x1.8p+0"
        @test f"\%a(bv)"     == "0x1.8p+0"
        @test f"\%#.0a(1.5)" == "0x2.p+0"
        @test f"\%+30a(1/3)" == "         +0x1.5555555555555p-2"
    end

    @testset "chars" begin
        @test f"\%c(65)"    == "A"
        @test f"\%c('A')"   == "A"
        @test f"\%3c('A')"  == "  A"
        @test f"\%-3c('A')" == "A  "
        @test f"\%c(248)"   == "ø"
        @test f"\%c('ø')"   == "ø"
    end

    @testset "type width specifier parsing (ignored)" begin
        @test f"\%llf(1.2)" == "1.200000"
        @test f"\%Lf(1.2)"  == "1.200000"
        @test f"\%hhu(1)"   == "1"
        @test f"\%hu(1)"    == "1"
        @test f"\%lu(1)"    == "1"
        @test f"\%llu(1)"   == "1"
        @test f"\%Lu(1)"    == "1"
        @test f"\%zu(1)"    == "1"
        @test f"\%ju(1)"    == "1"
        @test f"\%tu(1)"    == "1"
    end

    @testset "strings" begin
        s1 = "test"
        s2 = "tést"

        @test f"\%s(s1)"   == s1
        @test f"\%s(s2)"   == s2

        @test f"\%8s(s1)"  == "    test"
        @test f"\%-8s(s1)" == "test    "

        # This was broken in v0.5
        @static if VERSION >= v"0.6-"
            @test f"\%8.3s(s1)"    == "     tes"
            @test f"\%#8.3s(s1)"   == "     \"te"
            @test f"\%-8.3s(s1)"   == "tes     "
            @test f"\%#-8.3s(s1)"  == "\"te     "
            @test f"\%.3s(s1)"     == "tes"
            @test f"\%#.3s(s1)"    == "\"te"
            @test f"\%-.3s(s1)"    == "tes"
            @test f"\%#-.3s(s1)"   == "\"te"
        end
    end

    @testset "type based formatting" begin
        fmt_default!()        # resets all defaults

        # some basic functionality testing
        x = 1234.56789

        @test f"\%(x)"             == "1234.567890"
        @test f"\%(x;prec=2)"      == "1234.57"
        @test f"\%(x,10,3)"        == "  1234.568"
        @test f"\%(x,10,3,:left)"  == "1234.568  "
        @test f"\%(x,10,3,:ljust)" == "1234.568  "
        @test f"\%(x,:commas)"     == "1,234.567890"

        i = 1234567

        @test f"\%(i)" == "1234567"
        @test f"\%(i,:commas)" == "1,234,567"

        fmt_default!(Int, :commas, width = 12)
        @test f"\%(i)" == "   1,234,567"
        @test f"\%(x)" == "1234.567890"  # default hasn't changed

        fmt_default!(:commas)
        @test f"\%(i)" == "   1,234,567"
        @test f"\%(x)" == "1,234.567890"  # width hasn't changed, but added commas

        fmt_default!(Int) # resets Integer defaults
        @test f"\%(i)" == "1234567"
        @test f"\%(i,:commas)" == "1,234,567"
    end
end

@testset "Python Formatting" begin
    @testset "Format string" begin
        s = "abc"
        @test f"\{}(s)"     == "abc"
        @test f"\{s}(s)"    == "abc"
        @test f"\{2s}(s)"   == "abc"
        @test f"\{5s}(s)"   == "abc  "
        @test f"\{>5s}(s)"  == "  abc"
        @test f"\{*>5s}(s)" == "**abc"
        @test f"\{*<5s}(s)" == "abc**"
    end

    @testset "Format Char" begin
        @test f"\{}('c')"     == "c"
        @test f"\{c}('c')"    == "c"
        @test f"\{3c}('c')"   == "c  "
        @test f"\{>3c}('c')"  == "  c"
        @test f"\{*>3c}('c')" == "**c"
        @test f"\{*<3c}('c')" == "c**"
    end

    @testset "Format integer" begin
        @test f"\{}(1234)" == "1234"
        @test f"\{d}(1234)" == "1234"
        @test f"\{n}(1234)" == "1234"
        @test f"\{x}(0x2ab)" == "2ab"
        @test f"\{X}(0x2ab)" == "2AB"
        @test f"\{o}(0o123)" == "123"
        @test f"\{b}(0b1101)" == "1101"

        @test f"\{d}(0)" == "0"
        @test f"\{d}(9)" == "9"
        @test f"\{d}(10)" == "10"
        @test f"\{d}(99)" == "99"
        @test f"\{d}(100)" == "100"
        @test f"\{d}(1000)" == "1000"

        @test f"\{06d}(123)" == "000123"
        @test f"\{+6d}(123)" == "  +123"
        @test f"\{+06d}(123)" == "+00123"
        @test f"\{ d}(123)" == " 123"
        @test f"\{ 6d}(123)" == "   123"
        @test f"\{<6d}(123)" == "123   "
        @test f"\{>6d}(123)" == "   123"
        @test f"\{*<6d}(123)" == "123***"
        @test f"\{*>6d}(123)" == "***123"
        @test f"\{< 6d}(123)" == " 123  "
        @test f"\{<+6d}(123)" == "+123  "
        @test f"\{> 6d}(123)" == "   123"
        @test f"\{>+6d}(123)" == "  +123"

        @test f"\{+d}(-123)" == "-123"
        @test f"\{-d}(-123)" == "-123"
        @test f"\{ d}(-123)" == "-123"
        @test f"\{06d}(-123)" == "-00123"
        @test f"\{<6d}(-123)" == "-123  "
        @test f"\{>6d}(-123)" == "  -123"
    end

    @testset "Format floating point (f)" begin

        @test f"\{}(0.125)" == "0.125"
        @test f"\{f}(0.0)" == "0.000000"
        @test f"\{f}(0.001)" == "0.001000"
        @test f"\{f}(0.125)" == "0.125000"
        @test f"\{f}(1.0/3)" == "0.333333"
        @test f"\{f}(1.0/6)" == "0.166667"
        @test f"\{f}(-0.125)" == "-0.125000"
        @test f"\{f}(-1.0/3)" == "-0.333333"
        @test f"\{f}(-1.0/6)" == "-0.166667"
        @test f"\{f}(1234.5678)" == "1234.567800"
        @test f"\{8f}(1234.5678)" == "1234.567800"

        @test f"\{8.2f}(8.376)" == "    8.38"
        @test f"\{<8.2f}(8.376)" == "8.38    "
        @test f"\{>8.2f}(8.376)" == "    8.38"
        @test f"\{8.2f}(-8.376)" == "   -8.38"
        @test f"\{<8.2f}(-8.376)" == "-8.38   "
        @test f"\{>8.2f}(-8.376)" == "   -8.38"

        @test f"\{<08.2f}(8.376)" == "00008.38"
        @test f"\{>08.2f}(8.376)" == "00008.38"
        @test f"\{<08.2f}(-8.376)" == "-0008.38"
        @test f"\{>08.2f}(-8.376)" == "-0008.38"
        @test f"\{*<8.2f}(8.376)" == "8.38****"
        @test f"\{*>8.2f}(8.376)" == "****8.38"
        @test f"\{*<8.2f}(-8.376)" == "-8.38***"
        @test f"\{*>8.2f}(-8.376)" == "***-8.38"

        @test f"\{.2f}(0.999)" == "1.00"
        @test f"\{.2f}(0.996)" == "1.00"
        # Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
        @test (f"\{.2f}(0.995)" == "1.00" || f"\{.2f}(0.995)" == "0.99")
        @test f"\{.2f}(0.994)" == "0.99"
    end

    @testset "Format floating point (e)" begin

        @test f"\{E}(0.0)" == "0.000000E+00"
        @test f"\{e}(0.0)" == "0.000000e+00"
        @test f"\{e}(0.001)" == "1.000000e-03"
        @test f"\{e}(0.125)" == "1.250000e-01"
        @test f"\{e}(100/3)" == "3.333333e+01"
        @test f"\{e}(-0.125)" == "-1.250000e-01"
        @test f"\{e}(-100/6)" == "-1.666667e+01"
        @test f"\{e}(1234.5678)" == "1.234568e+03"
        @test f"\{8e}(1234.5678)" == "1.234568e+03"

        @test f"\{<12.2e}(13.89)" == "1.39e+01    "
        @test f"\{>12.2e}(13.89)" == "    1.39e+01"
        @test f"\{*<12.2e}(13.89)" == "1.39e+01****"
        @test f"\{*>12.2e}(13.89)" == "****1.39e+01"
        @test f"\{012.2e}(13.89)" == "00001.39e+01"
        @test f"\{012.2e}(-13.89)" == "-0001.39e+01"
        @test f"\{+012.2e}(13.89)" == "+0001.39e+01"

        @test f"\{.1e}(0.999)" == "1.0e+00"
        @test f"\{.1e}(0.996)" == "1.0e+00"
        # Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
        @test (f"\{.1e}(0.995)" == "1.0e+00" || f"\{.1e}(0.995)" == "9.9e-01")
        @test f"\{.1e}(0.994)" == "9.9e-01"
        @test f"\{.1e}(0.6)" == "6.0e-01"
        @test f"\{.1e}(0.9)" == "9.0e-01"
    end

    @testset "Format special floating point value" begin

        @test f"\{f}(NaN)" == "NaN"
        @test f"\{e}(NaN)" == "NaN"
        @test f"\{f}(NaN32)" == "NaN"
        @test f"\{e}(NaN32)" == "NaN"

        @test f"\{f}(Inf)" == "Inf"
        @test f"\{e}(Inf)" == "Inf"
        @test f"\{f}(Inf32)" == "Inf"
        @test f"\{e}(Inf32)" == "Inf"

        @test f"\{f}(-Inf)" == "-Inf"
        @test f"\{e}(-Inf)" == "-Inf"
        @test f"\{f}(-Inf32)" == "-Inf"
        @test f"\{e}(-Inf32)" == "-Inf"

        @test f"\{<5f}(Inf)" == "Inf  "
        @test f"\{>5f}(Inf)" == "  Inf"
        @test f"\{*<5f}(Inf)" == "Inf**"
        @test f"\{*>5f}(Inf)" == "**Inf"
    end

    @testset "test commas..." begin
        @test f"\%'d(1000)" == "1,000"
        @test f"\%'d(-1000)" == "-1,000"
        @test f"\%'d(100)" == "100"
        @test f"\%'d(-100)" == "-100"
        @test f"\%'f(Inf)" == "Inf"
        @test f"\%'f(-Inf)" == "-Inf"
        @test f"\%'s(1000.0)" == "1,000.0"
        @test f"\%'s(1234567.0)" == "1.234567e6"
    end # testset test commas
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 5dca09ac-523b-11eb-114b-af44490cde43
md"""
## Test1

This is an example notebook that includes markdowns and codes
- Equations: $ T = \alpha^2 $
- [julia](https://julialang.org)
"""

# ╔═╡ 700e8cc8-523b-11eb-2a95-8deb15bc31aa
function f(x,y)
    return x+y
end

# ╔═╡ 24642772-523a-11eb-1b5c-81b8a551b287
f(5,4)

# ╔═╡ 61cf7cd0-523a-11eb-1b81-b5430d6f9a59
# Blank line

# ╔═╡ d0d169f6-523b-11eb-3343-71773743c057
begin
	f(5,4)
	f(3,2)
end

# ╔═╡ d9dcff7e-523b-11eb-3d36-c17f4cda88bd
begin
	a = 5
	b = 6
	f(a,b)
end

# ╔═╡ Cell order:
# ╟─5dca09ac-523b-11eb-114b-af44490cde43
# ╠═700e8cc8-523b-11eb-2a95-8deb15bc31aa
# ╠═24642772-523a-11eb-1b5c-81b8a551b287
# ╠═61cf7cd0-523a-11eb-1b81-b5430d6f9a59
# ╠═d0d169f6-523b-11eb-3343-71773743c057
# ╠═d9dcff7e-523b-11eb-3d36-c17f4cda88bd


text:
using Compat

if VERSION >= v"0.6.0-dev.693"
    using Base.Broadcast: check_broadcast_indices, broadcast_indices
else
    using Base.Broadcast: check_broadcast_shape, broadcast_shape
    const check_broadcast_indices = check_broadcast_shape
    const broadcast_indices = broadcast_shape
end

if VERSION < v"0.6.0-dev" # Old approach needed for inference to work
    using Base: _default_eltype

    ftype(f, A) = typeof(f)
    ftype(f, A...) = typeof(a -> f(a...))
    ftype(T::DataType, A) = Type{T}
    ftype(T::DataType, A...) = Type{T}

    if isdefined(Base, :Iterators)
        using Base.Iterators: Zip2
    else
        using Base: Zip2
    end
    ziptype(A) = Tuple{eltype(eltype(A))}
    ziptype(A, B) = Zip2{Tuple{eltype(eltype(A))}, Tuple{eltype(eltype(B))}}
    @inline ziptype(A, B, C, D...) = Zip{Tuple{eltype(A)}, ziptype(B, C, D...)}

    nullable_broadcast_eltype(f, As...) =
        _default_eltype(Base.Generator{ziptype(As...), ftype(f, As...)})
else
    Base.@pure nullable_eltypestuple(a) = Tuple{eltype(eltype(a))}
    Base.@pure nullable_eltypestuple(T::Type) = Tuple{Type{eltype(T)}}
    Base.@pure nullable_eltypestuple(a, b...) =
        Tuple{nullable_eltypestuple(a).types..., nullable_eltypestuple(b...).types...}

    Base.@pure function nullable_broadcast_eltype(f, As...)
        T = Core.Inference.return_type(f, nullable_eltypestuple(As...))
        T === Union{} ? Any : T
    end
end

invoke_broadcast!{F, N}(f::F, dest, As::Vararg{NullableArray, N}) =
    invoke(broadcast!, Tuple{F, AbstractArray, Vararg{AbstractArray, N}}, f, dest, As...)

"""
    broadcast(f, As::NullableArray...)

Call `broadcast` with nullable lifting semantics and return a `NullableArray`.
Lifting means calling function `f` on the the values wrapped inside `Nullable` entries
of the input arrays, and returning null if any entry is missing.

Note that this method's signature specifies the source `As` arrays as all
`NullableArray`s. Thus, calling `broadcast` on arguments consisting
of both `Array`s and `NullableArray`s will fall back to the standard implementation
of `broadcast` (i.e. without lifting).
"""
function Base.broadcast{F}(f::F, As::NullableArray...)
    # These definitions are needed to avoid allocation due to splatting
    @inline f2(x1) = lift(f, (x1,))
    @inline f2(x1, x2) = lift(f, (x1, x2))
    @inline f2(x1, x2, x3) = lift(f, (x1, x2, x3))
    @inline f2(x1, x2, x3, x4) = lift(f, (x1, x2, x3, x4))
    @inline f2(x1, x2, x3, x4, x5) = lift(f, (x1, x2, x3, x4, x5))
    @inline f2(x1, x2, x3, x4, x5, x6) = lift(f, (x1, x2, x3, x4, x5, x6))
    @inline f2(x1, x2, x3, x4, x5, x6, x7) = lift(f, (x1, x2, x3, x4, x5, x6, x7))
    @inline f2(x...) = lift(f, x)

    T = nullable_broadcast_eltype(f, As...)
    dest = similar(NullableArray{T}, broadcast_indices(As...))
    invoke_broadcast!(f2, dest, As...)
end

"""
    broadcast!(f, dest::NullableArray, As::NullableArray...)

Call `broadcast!` with nullable lifting semantics.
Lifting means calling function `f` on the the values wrapped inside `Nullable` entries
of the input arrays, and returning null if any entry is missing.

Note that this method's signature specifies the destination `dest` array as well as the
source `As` arrays as all `NullableArray`s. Thus, calling `broadcast!` on a arguments
consisting of both `Array`s and `NullableArray`s will fall back to the standard implementation
of `broadcast!` (i.e. without lifting).
"""
function Base.broadcast!{F}(f::F, dest::NullableArray, As::NullableArray...)
    # These definitions are needed to avoid allocation due to splatting
    @inline f2(x1) = lift(f, (x1,))
    @inline f2(x1, x2) = lift(f, (x1, x2))
    @inline f2(x1, x2, x3) = lift(f, (x1, x2, x3))
    @inline f2(x1, x2, x3, x4) = lift(f, (x1, x2, x3, x4))
    @inline f2(x1, x2, x3, x4, x5) = lift(f, (x1, x2, x3, x4, x5))
    @inline f2(x1, x2, x3, x4, x5, x6) = lift(f, (x1, x2, x3, x4, x5, x6))
    @inline f2(x1, x2, x3, x4, x5, x6, x7) = lift(f, (x1, x2, x3, x4, x5, x6, x7))
    @inline f2(x...) = lift(f, x)

    invoke_broadcast!(f2, dest, As...)
end

# To fix ambiguity
function Base.broadcast!{F}(f::F, dest::NullableArray)
    f2() = lift(f)
    invoke_broadcast!(f2, dest)
end

# broadcasted ops
if VERSION < v"0.6.0-dev.1632"
    for (op, scalar_op) in (
        (:(@compat Base.:(.==)), :(==)),
        (:(@compat Base.:.!=), :!=),
        (:(@compat Base.:.<), :<),
        (:(@compat Base.:.>), :>),
        (:(@compat Base.:.<=), :<=),
        (:(@compat Base.:.>=), :>=)
    )
        @eval begin
            ($op)(X::NullableArray, Y::NullableArray) = broadcast($scalar_op, X, Y)
        end
    end
end


text:
using JuliaSet
using Base.Test

# write your own tests here
@test 1 == 1

x = collect((-10:10)/100);
y = collect((-10:10)/100);
n_iter = 10;
escape_tol = 5;

function R(z)
	return z
end

A=GenerateJuliaSet(R,x,y,n_iter);

@test all(A.==escape_tol+1);

max_stars_repo_name: restrin/JuliaSet

max_issues_repo_name: restrin/JuliaSet

max_forks_repo_name: restrin/JuliaSet


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

function testappend(nm, t, writekw, readkw, extratests)
    println("testing append: $nm")
    io = Arrow.tobuffer(t; writekw...)
    bytes = read(io)
    mktemp() do path, io
        write(io, bytes)
        close(io)

        t1 = Arrow.Table(read(path); readkw...)
        f1 = first(Tables.columns(t1))
        Arrow.append(path, t1; writekw..., readkw...)
        nparts = 0
        for t2 in Arrow.Stream(path)
            @test isequal(f1, first(Tables.columns(t2)))
            nparts += 1
        end
        @test nparts == 2
    end
end

function testappend_compression(compression_option)
    mktempdir() do path
        testdata = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
        file1 = joinpath(path, "table1.arrow")
        file2 = joinpath(path, "table2.arrow")

        open(file1, "w") do io
            Arrow.write(io, testdata; file=false, compress=compression_option)
        end
        isstream, schema, compression = open(Arrow.stream_properties, file1)
        @test isstream
        @test compression == compression_option

        open(file2, "w") do io
            Arrow.write(io, testdata; file=false)
        end

        arrow_table2 = Arrow.Table(file2)
        arrow_table2 |> Arrow.append(file1)
        arrow_table1 = Arrow.Table(file1)

        isstream, schema, compression = open(Arrow.stream_properties, file1)
        @test isstream
        @test compression == compression_option

        @test length(Tables.columns(arrow_table1)[1]) == 20
        @test length(Tables.columns(arrow_table2)[1]) == 10
    end
end

function testappend_partitions()
    mktempdir() do path
        testdata = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
        file1 = joinpath(path, "table1.arrow")
        file2 = joinpath(path, "table2.arrow")
        open(file1, "w") do io
            Arrow.write(io, testdata; file=false)
        end
        arrow_table1 = Arrow.Table(file1)
        isstream, schema, compression = open(Arrow.stream_properties, file1)
        @test isstream
        @test compression === nothing
        @test schema.names == (:col1,)
        @test schema.types == (Int64,)

        # can only append to arrow stream
        open(file2, "w") do io
            Arrow.write(io, testdata; file=true)
        end
        @test_throws ArgumentError Arrow.append(file2, arrow_table1)

        # schema must match
        testdata2 = (col2=Int64[1,2,3,4,5,6,7,8,9,10],)
        open(file2, "w") do io
            Arrow.write(io, testdata2; file=false)
        end
        @test_throws ArgumentError Arrow.append(file2, arrow_table1)

        # recreate file2 in arrow format with correct schema
        open(file2, "w") do io
            Arrow.write(io, testdata; file=false)
        end

        # start
        # arrow_table1: 1 partition, 10 rows
        # arrow_table2: 1 partition, 10 rows
        arrow_table2 = Arrow.Table(file2)
        @test length(Tables.columns(arrow_table1)[1]) == 10
        @test length(Tables.columns(arrow_table2)[1]) == 10

        @test_throws ArgumentError Arrow.append(file1, arrow_table2; ntasks = -1)
        arrow_table2 |> Arrow.append(file1)
        arrow_table1 = Arrow.Table(file1)
        # now
        # arrow_table1: 2 partitions, 20 rows
        # arrow_table2: 1 partition, 10 rows

        @test Tables.schema(arrow_table1) == Tables.schema(arrow_table2)
        @test length(Tables.columns(arrow_table1)[1]) == 20
        @test length(Tables.columns(arrow_table2)[1]) == 10
        @test length(collect(Tables.partitions(Arrow.Stream(file1)))) == 2 * length(collect(Tables.partitions(Arrow.Stream(file2))))

        Arrow.append(file2, arrow_table1; ntasks=1) # append with single task
        arrow_table2 = Arrow.Table(file2)
        # now
        # arrow_table1: 2 partitions, 20 rows
        # arrow_table2: 2 partitions, 30 rows (both partitions of table1 are appended as single partition)

        @test Tables.schema(arrow_table1) == Tables.schema(arrow_table2)
        @test length(Tables.columns(arrow_table1)[1]) == 20
        @test length(Tables.columns(arrow_table2)[1]) == 30
        @test length(collect(Tables.partitions(Arrow.Stream(file1)))) == length(collect(Tables.partitions(Arrow.Stream(file2))))

        Arrow.append(file1, Arrow.Stream(file2))
        arrow_table1 = Arrow.Table(file1)
        # now
        # arrow_table1: 4 partitions, 50 rows (partitions of file2 stream are appended without being merged)
        # arrow_table2: 2 partitions, 30 rows

        @test Tables.schema(arrow_table1) == Tables.schema(arrow_table2)
        @test length(Tables.columns(arrow_table1)[1]) == 50
        @test length(Tables.columns(arrow_table2)[1]) == 30
        @test length(collect(Tables.partitions(Arrow.Stream(file1)))) == 2 * length(collect(Tables.partitions(Arrow.Stream(file2))))
    end
end

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
using Test
using Distributed
if nworkers() == 1
	addprocs(2);
elseif nworkers() < 2
	addprocs(2 - nworkers());
end

using SparseArrays
using LinearAlgebra
@everywhere using jInv.Mesh
@everywhere using jInv.LinearSolvers
@everywhere using jInv.ForwardShare
@everywhere using jInv.InverseSolve
@everywhere using jInv.Utils
@everywhere using Printf
@everywhere using Helmholtz
@everywhere using jInvSeismic.FWI
@everywhere using jInvSeismic.Utils
@everywhere using jInvSeismic.EikonalInv
@everywhere using jInvSeismic.BasicFWI
@everywhere using Multigrid
using Statistics
using MAT
using DelimitedFiles

@testset "jInvSeismic" begin
	#include("BasicFWI/runtests.jl")
	include("EikonalInv/runtests.jl")
	include("FWI/runtests.jl")
end


text:
using Waterflow

function test_run(s::Float64, e::Float64)::Array{Float64, 1}

    a = Float64[0.0, 0.0]
    @ccall ihacres.test_ptr(a::Ptr{Cdouble}, s::Cdouble, e::Cdouble)::Cvoid

    return a
end

@info test_run(100.0, 200.0)

test_node = StreamNode{Float64}(
    100.0,  # area
    100.0,  # d
    0.1,  # d2
    0.1,  # e
    0.1,  # f
    54.35,  # a
    0.012,  # b
    2.5,  # storage_coef
    0.240195,  # alpha
    [100.0],  # storage
    [100.0],  # quickflow
    [100.0],  # slowflow
    [],  # outflow
    [],  # effective_rainfall
    [],  # et
    []   # inflow
)

@info "Outflow:" run_node(test_node, 6.0, 3.0, 50.0, 10.0)

@info test_node.slowflow

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


text:
using Base: @kwdef

"""
    Module

A data structure to store information to generate hardware for a circuit.
This structure can be manually modified if needed but typically [`@circuit`](@ref)
is used to auto-populate it.

Hardware generation traverses `dfg` and uses `handlers` to generate Verilog strings.

# Fields:
- `name::Symbol`: the name of the module
- `parameters::Dict{Symbol, Number}`: a map from the name of each parameter to its default value
- `submodules::Dict{Symbol, Symbol}`: a map from the name of each submodule to its type
- `dfg::MetaDiGraph{Int, Float64}`: a data flow graph representing the circuit to be generated
- `handlers::Dict{Operation, AbstractHandler}`: a map from operation type to a hardware generation handler.

See also: [HW.generate](@ref)
"""
@kwdef mutable struct Module
    name::Symbol
    parameters::Dict{Symbol, Number} = Dict{Symbol, Number}()
    submodules::Dict{Symbol, Symbol} = Dict{Symbol, Symbol}()
    dfg::MetaDiGraph{Int, Float64} = MetaDiGraph()
    handlers::Dict{Operation, AbstractHandler} = Dict{Operation, AbstractHandler}()
end

Base.show(io::IO, m::Module) =
    print(io, "Module $(m.name) with $(length(m.parameters)) parameters and $(length(m.submodules)) submodules.")
Base.show(io::IO, ::MIME"text/plain", m::Module) = print("""
    Module $(m.name):
        Parameters:
            $(m.parameters)
        Submodules:
            $(m.submodules)
        Number of operations: $(nv(m.dfg))
        Number of inputs: $(length(map(x -> getinputs(m.dfg, x), getroots(m.dfg))))
        Number of outputs: $(length(map(x -> getoutputs(m.dfg, x), getbuds(m.dfg))))
    """)

findnode(g::MetaDiGraph, inputs, outputs, operator) =
    collect(filter_vertices(g, (g, v) -> all(getname.(get_prop(g, v, :inputs)) .== inputs) &&
                                         all(getname.(get_prop(g, v, :outputs)) .== outputs) &&
                                         get_prop(g, v, :operator) == operator))
getroots(g::MetaDiGraph) = collect(filter_vertices(g, (g, v) -> isempty(inneighbors(g, v))))
getbuds(g::MetaDiGraph) = collect(filter_vertices(g, (g, v) -> isempty(outneighbors(g, v))))
getparents(g::MetaDiGraph, x::Variable) =
    filter_vertices(g, (g, v) -> x ∈ get_prop(g, v, :outputs))
getchildren(g::MetaDiGraph, x::Variable) =
    filter_vertices(g, (g, v) -> x ∈ get_prop(g, v, :inputs))
function traverse(g::MetaDiGraph, vs::Vector{T}, visited = T[]) where T
    levelup = unique(reduce(vcat, map(v -> outneighbors(g, v), vs)))
    parents = filter(v -> all(x -> x ∈ visited, inneighbors(g, v)), levelup)

    return parents, union(parents, visited)
end

function addnode!(m::Module, inputs::Vector{Variable}, outputs::Vector{Variable}, op::Symbol)
    add_vertex!(m.dfg)
    node = nv(m.dfg)
    set_prop!(m.dfg, node, :inputs, inputs)
    set_prop!(m.dfg, node, :outputs, outputs)
    set_prop!(m.dfg, node, :operator, op)

    for input in inputs
        for parent in getparents(m.dfg, input)
            add_edge!(m.dfg, parent, node)
        end
    end

    for output in outputs
        for child in getchildren(m.dfg, output)
            add_edge!(m.dfg, node, child)
        end
    end

    return m
end

function updatetype!(m::Module, inputs::Vector{Variable}, outputs::Vector{Variable}, op::Symbol)
    innames = getname.(inputs)
    outnames = getname.(outputs)
    v = findnode(m.dfg, innames, outnames, op)
    isempty(v) && error("""
        Could not update node because it was not found in DFG.
        inputs: $innames
        outputs: $outnames
        operator: $op
        """)
    set_prop!(m.dfg, v[1], :inputs, inputs)
    set_prop!(m.dfg, v[1], :outputs, outputs)
    set_prop!(m.dfg, v[1], :operator, op)

    return m
end

getinputs(g::MetaDiGraph, v) = get_prop(g, v, :inputs)
getoutputs(g::MetaDiGraph, v) = get_prop(g, v, :outputs)
getoperator(g::MetaDiGraph, v) = get_prop(g, v, :operator)

include("optimizations/constantreduction.jl")
include("optimizations/constantreplacement.jl")

function printdfg(m::Module)
    nodes = getroots(m.dfg)
    visited = nodes
    padding = ""

    while !isempty(nodes)
        for node in nodes
            inputs = getinputs(m.dfg, node)
            outputs = getoutputs(m.dfg, node)
            op = getoperator(m.dfg, node)
            println("$(padding)inputs: $inputs")
            println("$(padding)outputs: $outputs")
            println("$(padding)op: $op")
        end

        nodes, visited = traverse(m.dfg, nodes, visited)
        padding *= "   "
    end
end

"""
    HW.generate(m::Module, netlist::Netlist)
    HW.generate(m::Module, f)
    HW.generate(c::Tuple{Module, Function}, dut, args...)

Generate the Verilog implementation of the module.
Users will most likely call the last method form above.

# Fields:
- `m::Module`: the module to generate
- `netlist::Netlist`: the netlist for the circuit being generated
- `f`: a closure with one argument (a netlist) that calls a runtime information extraction function
- `c::Tuple{Module, Function}`: the tuple returned by [`@circuit`](@ref)
- `dut`: an instance of the circuit struct
- `args`: example arguments to circuit
"""
function generate(m::Module, netlist::Netlist)
    outstr = ""
    m = deepcopy(m)

    printdfg(m)
    println()

    # run constant reduction
    constantreduction!(m, netlist)
    printdfg(m)
    println()
    constantreplacement!(m, netlist)
    printdfg(m)
    println()

    # start at inputs
    nodes = getroots(m.dfg)
    visited = nodes

    while !isempty(nodes)
        # for each node, invoke the appropriate handler
        for node in nodes
            inputs = getinputs(m.dfg, node)
            outputs = getoutputs(m.dfg, node)
            op = Operation(gettype.(inputs), gettype.(outputs), getoperator(m.dfg, node))

            if haskey(m.handlers, op)
                handler = m.handlers[op]
            else
                handler = gethandler(op)()
                m.handlers[op] = handler
            end

            outstr *= handler(netlist, inputs, outputs)
        end

        # move one level up in the DFG
        nodes, visited = traverse(m.dfg, nodes, visited)
    end

    return outstr
end
function generate(m::Module, f)
    # get runtime info and populate netlist
    netlist = Netlist()
    f(netlist)

    return generate(m, netlist)
end
generate(c::Tuple{Module, Function}, dut, args...) = generate(c[1], (netlist) -> c[2](dut, c[1], netlist, args...))


text:
# This file includes two kinds of codes
#   - Codes for backward compatibility
#   - Glue codes that might nolonger be necessary in the future

# patch for issue #110
if isdefined(Base, :ComposedFunction) # Julia >= 1.6.0-DEV.85
    # https://github.com/JuliaLang/julia/pull/37517
    _round(tform::ComposedFunction; kwargs...) = _round(tform.outer; kwargs...) ∘ _round(tform.inner; kwargs...)
end

@static if !isdefined(Base, :IdentityUnitRange)
    const IdentityUnitRange = Base.Slice
else
    using Base: IdentityUnitRange
end

@static if VERSION < v"1.1"
    @inline isnothing(x) = x === nothing
end

# FIXME: upstream https://github.com/JuliaGraphics/ColorVectorSpace.jl/issues/75
@inline _nan(::Type{HSV{Float16}}) = HSV{Float16}(NaN16,NaN16,NaN16)
@inline _nan(::Type{HSV{Float32}}) = HSV{Float32}(NaN32,NaN32,NaN32)
@inline _nan(::Type{HSV{Float64}}) = HSV{Float64}(NaN,NaN,NaN)
@inline _nan(::Type{T}) where {T} = nan(T)

if !hasmethod(Constant{Nearest}, ())
    # `Constant{Nearest}()` is not defined for Interpolations <= v0.13.2
    # https://github.com/JuliaMath/Interpolations.jl/pull/426
    construct_interpolation_type(::Type{T}) where T<:Union{Linear, Constant} = T()
    construct_interpolation_type(::Type{Constant{Nearest}}) = Constant()
else
    construct_interpolation_type(::Type{T}) where T<:Union{Linear, Constant} = T()
end


text:
"""
	mask(cmd0::String="", arg1=[]; kwargs...)

Reads (length,azimuth) pairs from file and plot a windmask diagram.

Full option list at [`psmask`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html)

Parameters
----------

- **I** : **inc** : -- Str or Number --

    Set a fixed azimuth projection for masks [Default uses track azimuth, but see -A].
	[`-I`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#i)
- $(GMT.opt_R)

- $(GMT.opt_B)
- **C** : **end_clip_path** : -- Bool or [] --

    Mark end of existing clip path. No input file is needed.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#C)
- **D** : **dump** : -- Str --

    Dump the (x,y) coordinates of each clipping polygon to one or more output files
    (or stdout if template is not given).
	[`-D`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#d)
- **F** : **oriented_polygons** : -- Str or [] --

    Force clip contours (polygons) to be oriented so that data points are to the left (-Fl [Default]) or right (-Fr)
	[`-F`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#f)
- **G** : **fill** : -- Number or Str --

    Set fill shade, color or pattern for positive and/or negative masks [Default is no fill].
	[`-G`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#g)
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- **L** : **node_grid** : -- Str --

    Save the internal grid with ones (data constraint) and zeros (no data) to the named nodegrid.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#l)
- **N** : **invert** : -- Bool or [] --

    Invert the sense of the test, i.e., clip regions where there is data coverage.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#n)
- $(GMT.opt_P)
- **Q** : **cut_number** : -- Number or Str --

    Do not dump polygons with less than cut number of points [Dumps all polygons].
	[`-Q`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#q)
- **S** : **search_radius** : -- Number or Str --

    Sets radius of influence. Grid nodes within radius of a data point are considered reliable.
	[`-S`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#s)
- **T** : **tiles** : -- Bool or [] --

    Plot tiles instead of clip polygons. Use -G to set tile color or pattern. Cannot be used with -D.
	[`-T`](http://gmt.soest.hawaii.edu/doc/latest/psmask.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_r)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function mask(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("psmask", cmd0, arg1)	# Speedy mode
	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/12c")
	cmd, = parse_JZ(cmd, d)
	cmd  = parse_UVXY(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_swap_xy(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R, opt_i, opt_bi, opt_di)

	cmd = add_opt(cmd, 'C', d, [:C :end_clip_path])
	cmd = add_opt(cmd, 'D', d, [:D :dump])
	cmd = add_opt(cmd, 'F', d, [:F :oriented_polygons])
	cmd = add_opt(cmd, 'G', d, [:G :fill])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'L', d, [:L :node_grid])
	cmd = add_opt(cmd, 'N', d, [:N :invert])
	cmd = add_opt(cmd, 'Q', d, [:Q :cut_number])
	cmd = add_opt(cmd, 'S', d, [:S :search_radius])
	cmd = add_opt(cmd, 'T', d, [:T :tiles])

	cmd = finish_PS(d, cmd, output, K, O)
	return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "psmask", arg1)
end

# ---------------------------------------------------------------------------------------------------
mask!(cmd0::String="", arg1=[]; K=true, O=true,  first=false, kw...) =
	mask(cmd0, arg1; K=K, O=O,  first=first, kw...)

mask(arg1=[], cmd0::String=""; K=false, O=false,  first=true, kw...) =
	mask(cmd0, arg1; K=K, O=O,  first=first, kw...)

mask!(arg1=[], cmd0::String=""; K=true, O=true,  first=false, kw...) =
	mask(cmd0, arg1; K=K, O=O,  first=first, kw...)

psmask  = mask			# Alias
psmask! = mask!			# Alias


text:
using Nemo

# _Fast_ composed sums and composed products of polynomials,
# using the algorithm described in:
# "Fast computation of special resultants"
# by Bostan, Flajolet, Salvy, and Schost

# derivative of polynomial
derivative(c::Vector) = c[2:end] .* (1:length(c)-1)

function polyinv(coeffs::Vector, n)
	R, x = Nemo.PowerSeriesRing(Nemo.FlintQQ, n, "x")
	a = R(map(Nemo.FlintQQ, coeffs), length(coeffs), n, 0)
	ai = inv(a)
	return Nemo.fmpq[coeff(ai,i) for i=0:n-1]
end

# compute newton power series of polynomial given with coefficients coeff,
# in base field R,x.
# See fig.1 in reference
function to_newton(coeffs::Vector{BigInt},n,R,x)
	# first, make monic.
	coeffs = coeffs//coeffs[end]

	d = length(coeffs)-1
	a_cfs = reverse(derivative(coeffs))
	b_cfs = reverse(coeffs)

	# initialize power series polynomials
	a = R(map(Nemo.FlintQQ, a_cfs))
	b = R(map(Nemo.FlintQQ, b_cfs))
	b0 = R(polyinv(b_cfs, n))

	c  = truncate(a*b0, d)

	r = R()
	x_power = R(1)
	x_d = x^d

	l = round(Int64,floor(n/d))
	for j = 0 : l
		r += c*x_power
		x_power *= x_d
		c = -mullow(shift_right(b*c,d),b0,d)
	end
	return r
end

to_array(p) = Rational{BigInt}[Rational(coeff(p,i)) for i=0:Nemo.degree(p)]

# tr: traces i.e. newton series
# This algorithm is based on the Leverrier-Faddeev algorithm
# see: http://math.stackexchange.com/questions/405822/what-is-the-fastest-way-to-find-the-characteristic-polynomial-of-a-matrix
function from_newton(tr::Vector{T}) where {T<:Number}
	# special case
	if tr==[1]
		return [0,1]
	end
	n = length(tr)
	c = Array{T}(UndefInitializer(),n)
	c[end] = one(T)
	for k = 1 : n-1
		next_c = -sum(tr[2:(k+1)].*c[end-k+1:end])/k
		c[end-k] = next_c
	end
	return c
end

# Hadamard (element-wise) product of two polynomials
function hadm(p,q,R)
	n = max(Nemo.degree(p),Nemo.degree(q))
	R([Nemo.coeff(p,i)*Nemo.coeff(q,i) for i=0:n])
end

# composed product of two polynomials, given as coeffs p and q
function composed_product(p::Vector{BigInt},q::Vector{BigInt})
	# compute newton series
	n = (length(p)-1)*(length(q)-1)+1
	R, x = Nemo.PolynomialRing(Nemo.FlintQQ, "x")
	a = to_newton(p,n,R,x)
	b = to_newton(q,n,R,x)

	# multiply newton series and invert
	pq = from_newton(to_array(hadm(a,b,R)))

	# convert to integer and return
	return map(numerator, pq*lcm(map(denominator, pq)))
end

# composed sum of two polynomials, given as coeffs p and q
function composed_sum(p::Vector{BigInt},q::Vector{BigInt})
	# compute newton series
	n = (length(p)-1)*(length(q)-1)+1
	R, x = Nemo.PolynomialRing(Nemo.FlintQQ, "x")
	a = to_newton(p,n,R,x)
	b = to_newton(q,n,R,x)

	# exp series
	ee  = R([Nemo.FlintQQ(1//factorial(BigInt(i))) for i=0:(n-1)])
	eei = R([Nemo.FlintQQ(   factorial(BigInt(i))) for i=0:(n-1)])

	# multiply newton series and invert
	m = truncate(hadm(a,ee,R)*hadm(b,ee,R),n)
	pq = from_newton(to_array(hadm(m,eei,R)))

	# convert to integer and return
	return map(numerator, pq*lcm(map(denominator, pq)))
end


text:
module DataDrivenDiffEq

using LinearAlgebra
using DiffEqBase
using ModelingToolkit
using QuadGK
using Statistics
using DSP
using FiniteDifferences, DataInterpolations
using Compat

abstract type AbstractBasis end;
abstract type AbstractKoopmanOperator end;

include("./optimizers/Optimize.jl")
using .Optimize

export set_threshold!, set_threshold
export STRRidge, ADMM, SR3

export ADM

export WeightedSum, WeightedExponentialSum, GoalProgramming
export weights

include("./basis.jl")
export Basis
export variables, jacobian, dynamics
export free_parameters, parameters, variables

include("./koopman/algorithms.jl")
export DMDPINV, DMDSVD, TOTALDMD

include("./koopman/koopman.jl")
export eigen, eigvals, eigvecs
export modes, frequencies
export is_discrete, is_continouos
export operator, generator
export inputmap, outputmap, updateable, isstable

include("./koopman/linearkoopman.jl")
export LinearKoopman, update!

include("./koopman/nonlinearkoopman.jl")
export NonlinearKoopman, reduce_basis

include("./koopman/exact_dmd.jl")
export DMD, gDMD

include("./koopman/dmdc.jl")
export DMDc, gDMDc

include("./koopman/extended_dmd.jl")
export EDMD, gEDMD

#include("./exact_dmd.jl")
#export ExactDMD
#export eigen, eigvals, eigvecs
#export modes, frequencies, isstable
#export dynamics, update!
#
#include("./extended_dmd.jl")
#export ExtendedDMD
#export dynamics, linear_dynamics
#export reduce_basis, update!

#include("./dmdc.jl")
#export DMDc
#export eigen, eigvals, eigvecs
#export get_dynamics, get_input_map, dynamics

include("./sindy/results.jl")
export SparseIdentificationResult
export print_equations
export get_coefficients, get_error, get_sparsity, get_aicc

include("./sindy/sindy.jl")
export SInDy
export sparse_regression, sparse_regression!

include("./sindy/isindy.jl")
export ISInDy

include("./utils.jl")
export AIC, AICC, BIC
export optimal_shrinkage, optimal_shrinkage!
export savitzky_golay
export burst_sampling, subsample

end # module


text:
@testset "Topostructures" begin
  s1 = ExplicitStructure(connect.([(1,2,3),(4,3,2)], Triangle))
  @test connectivities(s1) == connect.([(1,2,3),(4,3,2)], Triangle)

  s2 = convert(HalfEdgeStructure, s1)
  @test ncells(s2) == 2
  @test all(edgeoncell(s2, c).cell == c for c in 1:ncells(s2))
  @test edgeoncell(s2, 1).head ∈ 1:3
  @test edgeoncell(s2, 2).head ∈ 4:6

  s3 = convert(ExplicitStructure, s2)
  @test s3 == s1
end


text:
### evolving graphs with integer nodes and timestamps

g = EvolvingGraph()
add_edge!(g, 1, 2, 1)
add_edge!(g, 1, 3, 1)
add_edge!(g, 1, 4, 2)
add_edge!(g, 1, 2, 2)
add_edge!(g, 2, 1, 3)
add_edge!(g, 2, 3, 3)
add_edge!(g, 1, 2, 1)

@test is_directed(g)
ns = nodes(g)
@test Node(1, 1) in ns
@test num_nodes(g) == 4
@test num_edges(g) == 6
@test rand(1:3) in timestamps(g)
@test num_timestamps(g) == 3

g = EvolvingGraph(is_directed = false)
add_edge!(g, 1, 2, 1)
add_edge!(g, 1, 3, 1)
add_edge!(g, 1, 4, 2)
add_edge!(g, 1, 2, 2)
add_edge!(g, 2, 1, 3)
add_edge!(g, 2, 3, 3)
add_edge!(g, 1, 2, 1)
@test is_directed(g) == false
@test num_nodes(g) == 4
@test num_edges(g) == 12

#### general evolving graphs

aa = ['a', 'b', 'c', 'c', 'a']
bb = ['b', 'a', 'a', 'b', 'b']
tt = ["t1", "t2", "t3", "t4", "t5"]
gg = evolving_graph_from_arrays(aa, bb, tt, is_directed = false)
display(gg)
nodes(gg)
N, T = eltype(gg)
# @test forward_neighbors(gg, 'c', "t4") == [('b', "t4")]

@test num_nodes(gg) == 3
edges(gg)
edges(gg, "t1")
@test num_edges(gg) == 10
timestamps(gg)
@test num_timestamps(gg) == 5

# convert to matrix
@test adjacency_matrix(gg, "t2") == full(sparse_adjacency_matrix(gg, "t2"))

g2 = EvolvingGraph{Node{Int},Char}(is_directed = true)
add_edge!(g2, 1, 2, 'a')
add_edge!(g2, 2, 3, 'b')
@test num_nodes(g2) == 3
@test num_edges(g2) == 2
@test num_timestamps(g2) == 2

# remove edge
g = EvolvingGraph{Node{Int}, String}()
add_edge!(g, 1, 2, "t1")
add_edge!(g, 2, 3, "t2")
add_edge!(g, 4, 2, "t2")
add_edge!(g, 4, 2, "t1")
add_edge!(g, 2, 1, "t3")
v = EvolvingGraphs.find_node(g, 1)
@test has_active_node(g, 1, "t1")

@test num_edges(g) == 5

add_bunch_of_edges!(g, [(1,3, "t1"), (1,4, "t2")])

n = num_edges(g)
@test is_directed(g)

forward_neighbors(g, active_nodes(g)[1])

wg = EvolvingGraph{Node{String}, Int}(is_directed = false, is_weighted = true)
add_edge!(wg, "a", "b", 1, weight = 2.)
add_edge!(wg, "b", "c", 1, weight = 1.)
add_edge!(wg, "a", "e", 2, weight = 4.)
add_edge!(wg, "e", "f", 3, weight = 2.)
@test num_edges(wg) == 8
display(wg)

@test adjacency_matrix(wg, 1) == full(sparse_adjacency_matrix(wg, 1))

wg = EvolvingGraph(is_weighted = true)
add_edge!(wg, 1, 2 ,1, weight = 1)
add_edge!(wg, 2, 1, 2, weight = 2)


text:
# ---
# title: 1319. Number of Operations to Make Network Connected
# id: problem1319
# author: Indigo
# date: 2021-01-23
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Union Find
# link: <https://leetcode.com/problems/number-of-operations-to-make-network-connected/description/>
# hidden: true
# ---
#
# There are `n` computers numbered from `0` to `n-1` connected by ethernet
# cables `connections` forming a network where `connections[i] = [a, b]`
# represents a connection between computers `a` and `b`. Any computer can reach
# any other computer directly or indirectly through the network.
#
# Given an initial computer network `connections`. You can extract certain
# cables between two directly connected computers, and place them between any
# pair of disconnected computers to make them directly connected. Return the
# _minimum number of times_ you need to do this in order to make all the
# computers connected. If it's not possible, return -1.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/02/sample_1_1677.png)**
#
#
#
#     Input: n = 4, connections = [[0,1],[0,2],[1,2]]
#     Output: 1
#     Explanation: Remove cable between computer 1 and 2 and place between computers 1 and 3.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/02/sample_2_1677.png)**
#
#
#
#     Input: n = 6, connections = [[0,1],[0,2],[0,3],[1,2],[1,3]]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: n = 6, connections = [[0,1],[0,2],[0,3],[1,2]]
#     Output: -1
#     Explanation: There are not enough cables.
#
#
# **Example 4:**
#
#
#
#     Input: n = 5, connections = [[0,1],[0,2],[3,4],[2,3]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#   * `1 <= connections.length <= min(n*(n-1)/2, 10^5)`
#   * `connections[i].length == 2`
#   * `0 <= connections[i][0], connections[i][1] < n`
#   * `connections[i][0] != connections[i][1]`
#   * There are no repeated connections.
#   * No two computers are connected by more than one cable.
#
#
## @lc code=start
using LeetCode

function make_connected(n::Int, connections::Vector{Vector{Int}})
    if length(connections) < n - 1
        return -1
    end
    father = collect(1:n)
    find_root(u::Int)::Int = (father[u] == u) ? u : (father[u] = find_root(father[u]))
    issame_root(u::Int, v::Int)::Bool = find_root(u) == find_root(v)
    function merge(u::Int, v::Int)
        u_root = find_root(u)
        v_root = find_root(v)
        if u_root != v_root
            n -= 1
            father[u_root] = v_root
        end
    end
    for connection in connections
        merge(connection[1], connection[2])
    end
    return n - 1
end
## @lc code=end


text:
struct FooBar end

# default case: error out
FooBar(::T) where T = FooBar(T)
FooBar(T::Type) =
    error("Type $T doesn't implement the FooBar interface.")

add_foo_and_bar(x) = add_foo_and_bar(FooBar(x), x)
add_foo_and_bar(::FooBar, x) = foo(x) + bar(x)

register_foobar(T::Type) =
    if hasmethod(foo, Tuple{T}) && hasmethod(bar, Tuple{T})
        @eval FooBar(::Type{$T}) = FooBar()
    else
        error("Type $T must implement `foo` and `bar` methods")
    end

max_stars_repo_name: ninjaaron/dispatching-design-patterns

max_issues_repo_name: ninjaaron/dispatching-design-patterns

max_forks_repo_name: ninjaaron/dispatching-design-patterns


text:
## Tests relating to data integrity. ##

import Currencies.CurrencyData

@testset "Data" begin

# long symbols should be unique
alllongsymbols = values(CurrencyData.LONG_SYMBOL)
@test length(Set(alllongsymbols)) == length(alllongsymbols)

# symbols should all be registered currencies
registeredcurrencies = keys(CurrencyData.ISO4217)
for sym in keys(CurrencyData.SHORT_SYMBOL)
    @test sym ∈ registeredcurrencies
end
for sym in keys(CurrencyData.LONG_SYMBOL)
    @test sym ∈ registeredcurrencies
end

# locale keys should be registered currencies, and value :before or :after
for (sym, val) in CurrencyData.LOCAL_SYMBOL_LOCATION
    @test sym ∈ registeredcurrencies
    @test val ∈ (:before, :after)
end

end  # @testset "Data"


text:
import Gadfly
Gadfly.set_default_plot_size(20Gadfly.cm, 12Gadfly.cm)

# make energy conservation and velocity distribution plots

temps_1069 = readdlm("q2/temps-rho-0.75-temp-1.069.dat")
temps_1304 = readdlm("q2/temps-rho-0.75-temp-1.304.dat")

speeds_1069 = readdlm("q2/speeds-rho-0.75-temp-1.069.dat")
speeds_1304 = readdlm("q2/speeds-rho-0.75-temp-1.304.dat")

ke_1069 = readdlm("q2/ke-rho-0.75-temp-1.069.dat")
ke_1304 = readdlm("q2/ke-rho-0.75-temp-1.304.dat")

pe_1069 = readdlm("q2/pe-rho-0.75-temp-1.069.dat")
pe_1304 = readdlm("q2/pe-rho-0.75-temp-1.304.dat")

e_1069 = pe_1069 + ke_1069
e_1304 = pe_1304 + ke_1304

function mdplot(values::Array, title::String, name::String, N::Integer)
    p = Gadfly.plot(
        x=[10:10:N],
        y=values,
        Gadfly.Geom.line,
        Gadfly.Guide.xlabel("Simulation Step"),
        Gadfly.Guide.ylabel(name),
        Gadfly.Guide.title("$title vs. Simulation Step for Liquid Argon")
    )

    # save as svg and png
    Gadfly.draw(Gadfly.SVG("$name.svg", 20Gadfly.cm, 12Gadfly.cm), p)
    Gadfly.draw(Gadfly.PNG("$name.png", 12Gadfly.cm, 8Gadfly.cm), p)

    return p
end

# m = 1, k = 1 in our natural units
function maxboltzbuilder(T::Real)
    return v -> √(2/π*T^3) * v.^2 .* exp(-v.^2 / (2T))
end

function mdplotspeeds(speeds::Array, name::String, T::Real)
    xs = [0:0.01:1]
    p = Gadfly.plot(
        Gadfly.layer(
            x = speeds,
            Gadfly.Geom.histogram(bincount=100, density=true),
            color = ["Measured speed histogram for T = $T"]
        ),
###         Gadfly.layer(
###             x = xs,
###             y = maxboltzbuilder(T)(xs),
###             Gadfly.Geom.line,
###             color = ["Predicted by Maxwell-Boltzmann"]
###         ),
###         color = ["Legend"],
        Gadfly.Guide.title("Temperature Distribution for T=$T vs. Boltzmann Distribution"),
        Gadfly.Guide.xlabel("Velocity"),
        Gadfly.Guide.ylabel("Probability"),
        Gadfly.Guide.xticks(ticks=[0:0.1:1])
    )

    # save as svg and png
    Gadfly.draw(Gadfly.SVG("$name.svg", 20Gadfly.cm, 12Gadfly.cm), p)
    Gadfly.draw(Gadfly.PNG("$name.png", 12Gadfly.cm, 8Gadfly.cm), p)

    return p
end

speed_plots = Gadfly.Plot[]
push!(speed_plots, mdplotspeeds(speeds_1069[401:end], "temp-dist-1.069", 1.069));
push!(speed_plots, mdplotspeeds(speeds_1304[401:end], "temp-dist-1.304", 1.304));

temp_plots = Gadfly.Plot[]
push!(temp_plots, mdplot(temps_1069, "Temperatures (T = 1.069)", "temp-1.069", 20000))
push!(temp_plots, mdplot(temps_1304, "Temperatures (T = 1.304)", "temp-1.304", 20000))

ke_plots = Gadfly.Plot[]
push!(ke_plots, mdplot(ke_1069, "Kinetic Energies (T = 1.069)", "temp-1.069", 20000))
push!(ke_plots, mdplot(ke_1304, "Kinetic Energies (T = 1.304)", "temp-1.304", 20000))

pe_plots = Gadfly.Plot[]
push!(pe_plots, mdplot(pe_1069, "Potential Energies (T = 1.069)", "temp-1.069", 20000))
push!(pe_plots, mdplot(pe_1304, "Potential Energies (T = 1.304)", "temp-1.304", 20000))

e_plots = Gadfly.Plot[]
push!(e_plots, mdplot(e_1069, "Total Energies (T = 1.069)", "temp-1.069", 20000))
push!(e_plots, mdplot(e_1304, "Total Energies (T = 1.304)", "temp-1.304", 20000))


text:
using ArcadeLearningEnvironment
using Logging
using CGP
using Images


text:
"""
    invert(a)

Return a new nested container by reversing the order of the nested container `a`, for
example turning a dictionary of arrays into an array of dictionaries, such that
`a[i][j] === invert(a)[j][i]`.

Note that in order for the keys of the inner and outer structure to be known, the input
container `a` must not be empty.

# Examples

```julia
julia> invert([[1,2,3],[4,5,6]])
3-element Array{Array{Int64,1},1}:
 [1, 4]
 [2, 5]
 [3, 6]

julia> invert((a = [1, 2, 3], b = [2.0, 4.0, 6.0]))
3-element Array{NamedTuple{(:a, :b),Tuple{Int64,Float64}},1}:
 (a = 1, b = 2.0)
 (a = 2, b = 4.0)
 (a = 3, b = 6.0)
```
""" # Array-Array
@inline function invert(a::AbstractArray{T}) where {T <: AbstractArray}
    f = first(a)
    innersize = size(a)
    outersize = size(f)
    innerkeys = keys(a)
    outerkeys = keys(f)

    @boundscheck for x in a
        if size(x) != outersize
            error("keys don't match")
        end
    end

    out = Array{Array{eltype(T),length(innersize)}}(undef, outersize)
    @inbounds for i in outerkeys
        out[i] = Array{eltype(T)}(undef, innersize)
    end

    return _invert!(out, a, innerkeys, outerkeys)
end

@propagate_inbounds function invert!(out::AbstractArray, a::AbstractArray)
    innerkeys = keys(a)
    outerkeys = keys(first(a))

    @boundscheck for x in a
        if keys(x) != outerkeys
            error("keys don't match")
        end
    end

    @boundscheck if keys(out) != outerkeys
        error("keys don't match")
    end

    @boundscheck for x in out
        if keys(x) != innerkeys
            error("keys don't match")
        end
    end

    return _invert!(out, a, innerkeys, outerkeys)
end

# Note: keys are assumed verified already
function _invert!(out, a, innerkeys, outerkeys)
    @inbounds for i ∈ innerkeys
        tmp = a[i]
        for j ∈ outerkeys
            out[j][i] = tmp[j]
        end
    end
    return out
end

# Tuple-tuple
@inline function invert(a::NTuple{n, NTuple{m, Any}}) where {n, m}
    if @generated
        exprs = [:(tuple($([:(a[$j][$i]) for j in 1:n]...))) for i in 1:m]
        return :(tuple($(exprs...)))
    else
        ntuple(i -> ntuple(j -> a[j][i], Val(n)), Val(m))
    end
end

# Tuple-Array
@inline function invert(a::NTuple{n, AbstractArray}) where {n}
    arrayinds = keys(a[1])

    @boundscheck for x in a
        if keys(x) != arrayinds
            error("indices are not uniform")
        end
    end

    T = _eltypes(typeof(a))
    out = similar(first(a), T)

    @inbounds invert!(out, a)
end

@inline function invert!(out::AbstractArray{<:NTuple{n, Any}}, a::NTuple{n, AbstractArray}) where n
    @boundscheck for x in a
        if keys(x) != keys(out)
            error("indices do not match")
        end
    end

    if @generated
        return quote
            @inbounds for i in keys(out)
                out[i] = $(:(tuple($([:( a[$j][i] ) for j in 1:n]...))))
            end

            return out
        end
    else
        @inbounds for i in keys(out)
            out[i] = map(x -> @inbounds(x[i]), a)
        end

        return out
    end
end

# Note that T is a concrete type, so the AbstractArrays should have all their type parameters
@pure function _eltypes(::Type{T}) where {T <: Tuple{Vararg{AbstractArray}}}
    types = T.parameters
    eltypes = map(eltype, types)
    return Tuple{eltypes...}
end

struct Indexer{i}
end

Indexer(i::Int) = Indexer{i}()
(::Indexer{i})(x) where {i} = @inbounds x[i]

# Array-Tuple
@inline function invert(a::AbstractArray{<:NTuple{n, Any}}) where {n}
    if @generated
        exprs = [ :(map($(Indexer(i)), a)) for i in 1:n ]
        return :( tuple($(exprs...)) )
    else
        ntuple(i -> map(x -> x[i], a), Val(n))
    end
end

# NamedTuple-NamedTuple
@inline function invert(a::NamedTuple{names1, <:Tuple{Vararg{NamedTuple{names2}}}}) where {names1, names2}
    if @generated
        exprs = [:(NamedTuple{names1}(tuple($([:(getfield(getfield(a, $(QuoteNode(n1))), $(QuoteNode(n2)))) for n1 in names1]...)))) for n2 in names2]
        return :(NamedTuple{names2}(tuple($(exprs...))))
    else
        NamedTuple{names2}(ntuple(i -> NamedTuple{names1}(ntuple(j -> getfield(getfield(a, names1[j]), names2[i]), Val(length(names1)))), Val(length(names2))))
    end
end

# NamedTuple-Array
@inline function invert(a::NamedTuple{names, <:Tuple{Vararg{AbstractArray}}}) where {names}
    arrayinds = keys(a[1])

    @boundscheck for x in a
        if keys(x) != arrayinds
            error("indices are not uniform")
        end
    end

    T = _eltypes(typeof(a))
    out = similar(first(a), NamedTuple{names, T})

    @inbounds invert!(out, a)
end

@inline function invert!(out::AbstractArray{<:NamedTuple{names}}, a::NamedTuple{names, <:Tuple{Vararg{AbstractArray}}}) where names
    @boundscheck for x in a
        if keys(x) != keys(out)
            error("indices do not match")
        end
    end

    if @generated
        return quote
            @inbounds for i in keys(out)
                out[i] = $(:(NamedTuple{names}(tuple($([:( a[$(QuoteNode(n))][i] ) for n in names]...)))))
            end

            return out
        end
    else
        @inbounds for i in keys(out)
            out[i] = map(x -> @inbounds(x[i]), a)
        end

        return out
    end
end

# TODO check that the AbstractArrays have all their type parameters
@pure function _eltypes(::Type{NamedTuple{names, T}}) where {names, T <: Tuple{Vararg{AbstractArray}}}
    types = T.parameters
    eltypes = map(eltype, types)
    return Tuple{eltypes...}
end

# Array-NamedTuple
struct GetFielder{name}
end

GetFielder(name::Symbol) = GetFielder{name}()
(::GetFielder{name})(x) where {name} = getfield(x, name)

@inline function invert(a::AbstractArray{<:NamedTuple{names}}) where {names}
    if @generated
        exprs = [ :(map($(GetFielder(n)), a)) for n in names ]
        return :( NamedTuple{names}(tuple($(exprs...))) )
    else
        NamedTuple{names}(ntuple(i -> map(x -> getfield(x, names[i]), a), Val(n)))
    end
end

# NamedTuple-Tuple
@inline function invert(a::NamedTuple{names1, <:Tuple{Vararg{Any, n}}}) where {names1, n}
    if @generated
        exprs = [:(NamedTuple{names1}(tuple($([:(getfield(getfield(a, $(QuoteNode(n1))), $i)) for n1 in names1]...)))) for i in 1:n]
        return :(tuple($(exprs...)))
    else
        ntuple(i -> NamedTuple{names1}(ntuple(j -> getfield(getfield(a, names1[j]), i), Val(length(names1)))), Val(n))
    end
end

# Tuple-NamedTuple
@inline function invert(a::Tuple{Vararg{NamedTuple{names2}, n}}) where {n, names2}
    if @generated
        exprs = [:(tuple($([:(getfield(getfield(a, $i), $(QuoteNode(n2)))) for i in 1:n]...))) for n2 in names2]
        return :(NamedTuple{names2}(tuple($(exprs...))))
    else
        NamedTuple{names2}(ntuple(i -> ntuple(j -> getfield(getfield(a, j), names2[i]), Val(n)), Val(length(names2))))
    end
end

# Dict-Dict

# Dict-Array

# Array-Dict

# Dict-Tuple

# Tuple-Dict

# NamedTuple-Dict

# Dict-NamedTuple


text:
using LinearAlgebra
using Plots
using Random
using BenchmarkTools
using JuMP

function kernelf(X1, X2, l=1.0, sigma_f=1.0)

    # Isotropic squared exponential kernel. Computes
    # a covariance matrix from points in X1 and X2.

    #Σ = zeros(size(X1,1), size(X2,1))
    Σ = Array{Float64,2}(undef, size(X1,1), size(X2,1))
    for i in 1:size(Σ,1)
        for j in 1:size(Σ,2)
            Σ[i,j] = sigma_f^2*exp(-(X1[i]-X2[j])^2/(2*l^2))
        end
    end

    return Σ

end

function plot_gp(μ, Σ, X, X_train=nothing, Y_train=nothing)
    σ = 1.96*sqrt.(diag(Σ))
    p = plot([X X], [μ μ], fillrange=[μ+σ μ-σ], fillalpha=0.3, c=:blue, label="")
    plot!(X, μ, c=:blue, label="GP mean")
    if X_train!=nothing && Y_train!=nothing
        plot!(X_train, Y_train, seriestype="scatter", color="red", label="samples")
    end
    display(p)
end

function kernelTest()
    X1 = collect(range(-5, step=1, 5))
    X2 = collect(range(-5, step=1, 5)) .+ 0.1.*randn(size(X1))
    plot(X1, label="original data")
    plot!(X2, label="corrupted data")
    μ = zeros(size(X1))
    Σ = kernelf(X1, X1)
    plot_gp(μ, Σ, X1)
end

function posteriorPredictor(X_s, X_train, Y_train, l=1.0, sigma_f=1.0, sigma_y=1e-8, K_inv=nothing)
    K_s = kernelf(X_train, X_s, l, sigma_f)
    K_ss = kernelf(X_s, X_s, l, sigma_f) + 1e-8*I
    if K_inv == nothing
        K = kernelf(X_train, X_train, l, sigma_f) + sigma_y^2*I
        K_inv = inv(K)
    end

    μ_s = K_s'*K_inv*Y_train
    #mu_s = K_s.T.dot(K_inv).dot(Y_train)
    Σ_s = K_ss - K_s'*K_inv*K_s
    #cov_s = K_ss - K_s.T.dot(K_inv).dot(K_s)

    return μ_s, Σ_s, K_inv
end

function testPredictor()
    lb = 0;
    ub = 5;
    n = 5;
    sigma_y = 0.0
    X_train = [-4 -3 -2 -1 1 4]'
    X_train = rand(n) .*(ub-lb) .+lb
    Y_train = sin.(X_train) + sigma_y*randn(size(X_train,1))
    X = collect(range(lb, ub, length=n))
    μ_s, Σ_s = posteriorPredictor(X, X_train, Y_train, 1.0, 1.0, sigma_y)
    plot_gp(μ_s, Σ_s, X, X_train, Y_train)
    plot!(X, sin.(X), lw=:2, color=:black, label="true function")
end

function hyperParameterOptimization(X_train, Y_train, noise, naive=true)
    model = Model()
    @variable(theta[1:2] >= 1e-5) #hyperparam variables
    if naive
    end
end

function f(x)
    return 8*sin.(x) + 1/2*x.^2 + 1/8*x.^3 + 1/32*x.^4 - 1/128*x.^5
end

#kernelTest()
#testPredictor()

lb = 1;
ub = 7;
n = 100;
l = 1.0;
sigma_f = 1.0
sigma_y = 0.9
#X_train = [-4 -3 -2 -1 1 4]'
#X_train = rand(n).*(ub-lb) .+ lb
X_train = 1.1*randn(n) .+ (ub - lb)/2 .+ lb
Y_train = f(X_train) + sigma_y*randn(size(X_train,1))
X = collect(range(lb, ub, length=n))
if sigma_y != 0.0
    μ_s, Σ_s, K_inv = posteriorPredictor(X, X_train, Y_train, l, sigma_f, sigma_y)
else
    μ_s, Σ_s, K_inv = posteriorPredictor(X, X_train, Y_train, l, sigma_f)
end
println("Timing training + inference")
@time posteriorPredictor(X, X_train, Y_train, 1.0, 1.0, sigma_y)
println("Timing inference")
@time posteriorPredictor(X, X_train, Y_train, 1.0, 1.0, sigma_y, K_inv)
plot_gp(μ_s, Σ_s, X, X_train, Y_train)
plot!(X, f(X), lw=:2, color=:black, label="true function", xlim=(lb,ub))

max_stars_repo_name: gbarsih/Safe-Optimal-Rendezvous

max_issues_repo_name: gbarsih/Safe-Optimal-Rendezvous

max_forks_repo_name: gbarsih/Safe-Optimal-Rendezvous


text:
module Scrypt

using Nettle
using SIMD

include("data/Salsa512.jl")
include("data/ScryptParameters.jl")
include("util.jl")

function scrypt(parameters::ScryptParameters, key::Vector{UInt8}, salt::Vector{UInt8}, derivedkeylength)
    derivedkeylength > 0 || ArgumentError("Must be > 0.") |> throw

    buffer = pbkdf2_sha256_1(key, salt, bufferlength(parameters))
    parallelbuffer = reshape(reinterpret(Salsa512, buffer), (elementblockcount(parameters), parameters.p))

    for i ∈ 1:parameters.p
        element = @views reshape(parallelbuffer[:, i], elementblockcount(parameters))
        smix!(element, parameters)
    end

    derivedkey = pbkdf2_sha256_1(key, buffer, derivedkeylength)
end

const HASH_LENGTH = 256 ÷ 8

function pbkdf2_sha256_1(key, salt::Vector{UInt8}, derivedkeylength)
    blockcount = cld(derivedkeylength, HASH_LENGTH)
    lastblockbytes = derivedkeylength - (blockcount - 1) * HASH_LENGTH

    salt = [salt; zeros(UInt8, 4)]
    salttail = view(salt, length(salt) - 3:length(salt))

    derivedkey = Matrix{UInt8}(undef, HASH_LENGTH, blockcount)

    for i ∈ 1:blockcount
        salttail[:] = reinterpret(UInt8, [UInt32(i)]) |> reverse
        derivedkey[:, i] = digest("sha256", key, salt)
    end

    derivedkey = reshape(derivedkey, blockcount * HASH_LENGTH)[1:derivedkeylength]
    return derivedkey
end

function smix!(element::AbstractVector{Salsa512}, parameters::ScryptParameters)
    workingbuffer = prepare(element)
    shufflebuffer = valloc(Salsa512, length(workingbuffer))
    scryptblock, workingbuffer, shufflebuffer = fillscryptblock!(workingbuffer, shufflebuffer, parameters.r, parameters.N)
    workingbuffer = mixwithscryptblock!(workingbuffer, scryptblock, shufflebuffer, parameters.r, parameters.N)
    restore!(element, workingbuffer)
end

const SALSA_BLOCK_REORDER_INDEXES = [13;  2;  7; 12;  1;  6; 11; 16;  5; 10; 15;  4;  9; 14;  3;  8]

function prepare(src::AbstractVector{Salsa512})
    dest = valloc(Salsa512, length(src))
    si = 1:length(src)
    dj = [2:length(dest); 1]

    for (i, j) ∈ zip(si, dj)
        dest[j] = src[i]
        permute!(uint32view(dest, j), SALSA_BLOCK_REORDER_INDEXES)
    end

    return dest
end #permute! is no faster than explicit vectorization, even with a few extra allocations

function restore!(dest::AbstractVector{Salsa512}, src::AbstractVector{Salsa512})
    si = 1:length(src)
    dj = [length(dest); 1:length(dest) - 1]

    for (i, j) ∈ zip(si, dj)
        dest[j] = src[i]
        invpermute!(uint32view(dest, j), SALSA_BLOCK_REORDER_INDEXES)
    end
end

function fillscryptblock!(workingbuffer::AbstractVector{Salsa512}, shufflebuffer::AbstractVector{Salsa512}, r, N)
    scryptblock = reshape(valloc(Salsa512, 2r * N), (2r, N))
    for i ∈ 1:N
        scryptelement = view(scryptblock, :, i)
        previousblock = lastblock = load_store!(workingbuffer, scryptelement, 1)
        for j ∈ 2:2r
            block = load_store!(workingbuffer, scryptelement, j)
            block = mixblock_shuffle_store!(block, previousblock, shufflebuffer, shuffleposition(j, r))
            previousblock = block
        end
        mixblock_shuffle_store!(lastblock, previousblock, shufflebuffer, 1)
        workingbuffer, shufflebuffer = shufflebuffer, workingbuffer
    end
    return scryptblock, workingbuffer, shufflebuffer
end

shuffleposition(j, halfblockcount) = (j - 2) ÷ 2 + 2 + (iseven(j) ? 0 : halfblockcount)

uint32view(x, i) = @inbounds reinterpret(UInt32, view(x, i:i))
vloadsalsa(x, i) = @inbounds vloada(Vec{16, UInt32}, uint32view(x, i), 1)
vloadsalsant(x, i) = @inbounds vloadnt(Vec{16, UInt32}, uint32view(x, i), 1)
vstoresalsa(v, x, i) = @inbounds vstorea(v, uint32view(x, i), 1)
vstoresalsant(v, x, i) = @inbounds vstorent(v, uint32view(x, i), 1)

function load_store!(workingbuffer::AbstractVector{Salsa512}, scryptelement::AbstractVector{Salsa512}, i)
    block = vloadsalsa(workingbuffer, i)
    vstoresalsant(block, scryptelement, i)
    return block
end

function mixwithscryptblock!(workingbuffer::AbstractVector{Salsa512}, scryptblock, shufflebuffer::AbstractVector{Salsa512}, r, N)
    for i ∈ 1:N
        n = integerify(workingbuffer, N)
        scryptelement = reshape(view(scryptblock, :, n), 2r)

        for j ∈ 1:r # prefetch first half of the element
            vprefetchnt(scryptelement, j)
        end

        previousblock = lastblock = load_xor(workingbuffer, scryptelement, 1)
        for j ∈ 2:2r
            if j ≤ (r + 1) # prefetch one additional block through end
                vprefetchnt(scryptelement, r + j - 1)
            end

            block = load_xor(workingbuffer, scryptelement, j)
            block = mixblock_shuffle_store!(block, previousblock, shufflebuffer, shuffleposition(j, r))
            previousblock = block
        end
        mixblock_shuffle_store!(lastblock, previousblock, shufflebuffer, 1)
        workingbuffer, shufflebuffer = shufflebuffer, workingbuffer
    end
    return workingbuffer
end

integerify(x::AbstractVector{Salsa512}, N) = uint32view(x, 1)[5] % N + 1

function load_xor(workingbuffer::AbstractVector{Salsa512}, scryptelement::AbstractVector{Salsa512}, i)
    block = vloadsalsa(workingbuffer, i)
    block ⊻= vloadsalsant(scryptelement, i)
    return block
end

function mixblock_shuffle_store!(block, previousblock, shufflebuffer, i)
    block ⊻= previousblock
    block = salsa20(block, 8)
    vstoresalsa(block, shufflebuffer, i)
    return block
end

function salsa20(block, iterations)
    inputblock = block

    lines = splitblock(block)
    for i ∈ 1:iterations
        lines = salsamix(lines)
        lines = salsatranspose(lines)
    end
    block = joinlines(lines)

    block += inputblock
    return block
end

splitblock(block) = (shufflevector(block, Val((0,1,2,3))),
                     shufflevector(block, Val((4,5,6,7))),
                     shufflevector(block, Val((8,9,10,11))),
                     shufflevector(block, Val((12,13,14,15))))

joinlines(lines) = @inbounds shufflevector(shufflevector(lines[1], lines[2], Val((0, 1, 2, 3, 4, 5, 6, 7))),
                                           shufflevector(lines[3], lines[4], Val((0, 1, 2, 3, 4, 5, 6, 7))),
                                           Val((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)))

function salsamix(lines)
    line1, line2, line3, line4 = lines
    line3 = salsa(line1, line2, line3, 7)
    line4 = salsa(line2, line3, line4, 9)
    line1 = salsa(line3, line4, line1, 13)
    line2 = salsa(line4, line1, line2, 18)
    return (line1, line2, line3, line4)
end

function salsa(addend1, addend2, xor_operand, rotationmagnitude)
    sum = addend1 + addend2
    rot = (sum << rotationmagnitude) | (sum >>> (sizeof(UInt32) * 8 - rotationmagnitude))
    return xor_operand ⊻ rot
end

function salsatranspose(lines)
    toline3 = @inbounds shufflevector(lines[1], Val((1, 2, 3, 0)))
    line1 = @inbounds shufflevector(lines[3], Val((3, 0, 1, 2)))
    line3 = toline3
    line4 = @inbounds shufflevector(lines[4], Val((2, 3, 0, 1)))
    return @inbounds (line1, lines[2], line3, line4)
end

export scrypt
export ScryptParameters

end


text:
"""
    initialize_beta!(gcm{Poisson_Bernoulli_VCModel})

Initialize the linear regression parameters `β` using GLM.jl
"""
function initialize_beta!(gcm::Poisson_Bernoulli_VCModel{T, VD, VL}) where {T <: BlasReal, VD, VL}
    # form df
    Xstack = []
    Y1stack = zeros(length(gcm.data))
    Y2stack = zeros(length(gcm.data))
    @inbounds for i in 1:length(gcm.data)
        push!(Xstack, gcm.data[i].X[1, 1:Integer((gcm.p)/2)])
        Y1stack[i] = gcm.data[i].y[1]
        Y2stack[i] = gcm.data[i].y[2]
    end
    X = vcat(transpose(Xstack)...)

    poisson_glm = GLM.glm(X, Y1stack, gcm.vecd[1][1], gcm.veclink[1][1])
    bernoulli_glm = GLM.glm(X, Y2stack, gcm.vecd[1][2], gcm.veclink[1][2])
    copyto!(gcm.β, [poisson_glm.pp.beta0; bernoulli_glm.pp.beta0])
    nothing
end

"""
    initialize_beta!(gcm{GLMCopulaVCModel})

Initialize the linear regression parameters `β` using GLM.jl
"""
function initialize_beta!(gcm::Union{GLMCopulaVCModel{T, D, Link}, GLMCopulaARModel{T, D, Link}, GLMCopulaCSModel{T, D, Link}, NBCopulaVCModel{T, D, Link}}) where {T <: BlasReal, D, Link}
    # form df
    Xstack = []
    Ystack = []
    @inbounds for i in 1:length(gcm.data)
        push!(Xstack, gcm.data[i].X)
        push!(Ystack, gcm.data[i].y)
    end
    Xstack = [vcat(Xstack...)][1]
    Ystack = [vcat(Ystack...)][1]
    fit_glm = GLM.glm(Xstack, Ystack, gcm.d[1], gcm.link[1])
    copyto!(gcm.β, fit_glm.pp.beta0)
    nothing
end

"""
    initialize_model!(gcm)

Initialize the linear regression parameters `β` using Newton's Algorithm under Independence Assumption, update variance components using MM-Algorithm.
"""
function initialize_model!(
    gcm::GLMCopulaARModel{T, D, Link}) where {T <: BlasReal, D, Link}
    println("initializing β using Newton's Algorithm under Independence Assumption")
    initialize_beta!(gcm)
    fill!(gcm.τ, 1.0)
    fill!(gcm.θ, 1.0)
    println("initializing variance components using MM-Algorithm")
    update_θ!(gcm)
    copyto!(gcm.σ2, gcm.θ)
    copyto!(gcm.ρ, 0.2)
    nothing
end

function initialize_model!(
    gcm::GLMCopulaCSModel{T, D, Link}) where {T <: BlasReal, D, Link}
    println("initializing β using Newton's Algorithm under Independence Assumption")
    initialize_beta!(gcm)
    copyto!(gcm.ρ, 0.2)
    nothing
end

"""
    initialize_model!(gcm{GLMCopulaVCModel, Poisson_Bernoulli_VCModel, NBCopulaVCModel})

Initialize the linear regression parameters `β` using GLM.jl, and update variance components using MM-Algorithm.
"""
function initialize_model!(
    gcm::Union{GLMCopulaVCModel{T, D, Link}, Poisson_Bernoulli_VCModel{T, VD, VL}}) where {T <: BlasReal, D, Link,  VD, VL}
    println("initializing β using Newton's Algorithm under Independence Assumption")
    initialize_beta!(gcm)
    @show gcm.β
    fill!(gcm.τ, 1.0)
    println("initializing variance components using MM-Algorithm")
    fill!(gcm.θ, 1.0)
    update_θ!(gcm)
    if sum(gcm.θ) >= 20
      fill!(gcm.θ, 1.0)
    end
    @show gcm.θ
    nothing
end

"""
    initialize_model!(gcm)
Initialize the linear regression parameters `β` and `τ=σ0^{-2}` by the least
squares solution.
"""
function initialize_model!(
    gcm::GaussianCopulaVCModel{T}
    ) where T <: BlasReal
    # accumulate sufficient statistics X'y
    xty = zeros(T, gcm.p)
    for i in eachindex(gcm.data)
        BLAS.gemv!('T', one(T), gcm.data[i].X, gcm.data[i].y, one(T), xty)
    end
    # least square solution for β
    ldiv!(gcm.β, cholesky(Symmetric(gcm.XtX)), xty)
    @show gcm.β
    # accumulate residual sum of squares
    rss = zero(T)
    for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rss += abs2(norm(gcm.data[i].res))
    end
    println("initializing dispersion using residual sum of squares")
    gcm.τ[1] = gcm.ntotal / rss
    @show gcm.τ
    println("initializing variance components using MM-Algorithm")
    fill!(gcm.θ, 1.0)
    update_θ!(gcm)
    @show gcm.θ
    nothing
end

# march 6 2022 try this initialize model instead so that r doesn't do weird things
function initialize_model!(
    gcm::NBCopulaVCModel{T, D, Link}) where {T <: BlasReal, D, Link}
    println("initializing β using GLM.jl")
    initialize_beta!(gcm) # initialize beta using negative binomial glm from GLM.jl
    @show gcm.β
    fill!(gcm.θ, 1.0)
    fill!(gcm.τ, 1.0)
    println("initializing variance components using MM-Algorithm")
    update_θ!(gcm)
    if sum(gcm.θ) >= 20
      fill!(gcm.θ, 1.0)
    end
    @show gcm.θ
    println("initializing r using Newton update")
    fill!(gcm.r, 1)
    QuasiCopula.update_r!(gcm)
    nothing
end

# # code inspired from https://github.com/JuliaStats/GLM.jl/blob/master/src/negbinfit.jl
function initialize_model!(gcm::NBCopulaCSModel{T, D, Link}) where {T <: BlasReal, D, Link}
    # initial guess for r = 1
    fill!(gcm.r, 1)

    # fit a Poisson regression model to estimate μ, η, β, τ
    println("Initializing NegBin r to Poisson regression values")
    nsample = length(gcm.data)
    gcsPoisson = Vector{GLMCopulaCSObs{T, Poisson{T}, LogLink}}(undef, nsample)
    for (i, gc) in enumerate(gcm.data)
      gcsPoisson[i] = GLMCopulaCSObs(gc.y, gc.X, Poisson(), LogLink())
    end
    gcmPoisson = GLMCopulaCSModel(gcsPoisson)
    initialize_model!(gcmPoisson) # initialize beta using poisson glm from GLM.jl
    copyto!(gcm.β, gcmPoisson.β)
    copyto!(gcm.σ2, gcmPoisson.σ2)

    # update r using maximum likelihood with Newton's method
    for gc in gcm.data
      fill!(gcm.τ, 1.0)
      fill!(gc.∇β, 0)
      fill!(gc.Hβ, 0)
      fill!(gc.varμ, 1)
      fill!(gc.res, 0)
    end
    println("initializing r using Newton update")
    QuasiCopula.update_r!(gcm)
    copyto!(gcm.ρ, 0.2)
    nothing
end

function initialize_model!(gcm::NBCopulaARModel{T, D, Link}) where {T <: BlasReal, D, Link}
    # initial guess for r = 1
    fill!(gcm.r, 1)

    # fit a Poisson regression model to estimate μ, η, β, τ
    println("Initializing NegBin r to Poisson regression values")
    nsample = length(gcm.data)
    gcsPoisson = Vector{GLMCopulaARObs{T, Poisson{T}, LogLink}}(undef, nsample)
    for (i, gc) in enumerate(gcm.data)
      gcsPoisson[i] = GLMCopulaARObs(gc.y, gc.X, Poisson(), LogLink())
    end
    gcmPoisson = GLMCopulaARModel(gcsPoisson)
    initialize_model!(gcmPoisson) # initialize beta using poisson glm from GLM.jl
    copyto!(gcm.β, gcmPoisson.β)
    copyto!(gcm.σ2, gcmPoisson.σ2)

    # update r using maximum likelihood with Newton's method
    for gc in gcm.data
      fill!(gcm.τ, 1.0)
      fill!(gc.∇β, 0)
      fill!(gc.Hβ, 0)
      fill!(gc.varμ, 1)
      fill!(gc.res, 0)
    end
    println("initializing r using Newton update")
    QuasiCopula.update_r!(gcm)
    copyto!(gcm.ρ, 0.2)
    nothing
end


text:
""" Returns the correct container spec for the selected type of JuMP Model"""
function _container_spec(m::M, axs...) where M<:JuMP.AbstractModel
    return JuMP.Containers.DenseAxisArray{JuMP.variable_type(m)}(undef, axs...)
end

function add_var_container!(canonical::Canonical, var_name::Symbol, axs...)
    canonical.variables[var_name] = _container_spec(canonical.JuMPmodel, axs...)
    return canonical.variables[var_name]
end

function add_cons_container!(canonical::Canonical, cons_name::Symbol, axs...)
    canonical.constraints[cons_name] = JuMPConstraintArray(undef, axs...)
    return canonical.constraints[cons_name]
end

function _add_param_container!(canonical::Canonical, param_reference::UpdateRef, axs...)
    canonical.parameters[param_reference] = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, axs...)
    return canonical.parameters[param_reference]
end

function _middle_rename(original::Symbol, split_char::String, addition::String)

    parts = split(String(original),split_char)

    return Symbol(parts[1],"_",addition,"_",parts[2])

end

function _remove_underscore(original::Symbol)

    if !occursin("_", String(original))
        return original
    end

    parts = split(String(original),"_")

    return parts[1]
end

"Replaces the string in `char` with the string`replacement`"
function replace_chars(s::String, char::String, replacement::String)
    return replace(s, Regex("[$char]") => replacement)
end

"Removes the string `char` from the original string"
function remove_chars(s::String, char::String)
    return replace_chars(s::String, char::String, "")
end


text:
using FractionalDiffEq, Plots

function test(du, u, p, t)
    du[1] = 1-4*u[1]+u[1]^2*u[2]
    du[2] = 3*u[1]-u[1]^2*u[2]
end

alpha = [0.8, 0.7]
t0=0; T=100
y0=[1.2; 2.8]
h=0.01
param=[1 3]
prob = FODESystem(test, alpha, y0, (t0, T))
#(t, y) = solve(prob, h, PECE())
(t, y) = solve(prob, h, PIEX())
plot(y[1, :], y[2, :])
#=
plot(t, y[1, :])
plot!(t, y[2, :])
=#