text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/extrema.jl
# Defines univariate and bivariate max and min operators.
#############################################################################

# Univariate min/max operators
@inline deriv_max(x::Float64, c::Float64) = (x < c) ? 0.0 : 1.0
@inline function cv_max(x::Float64, xL::Float64, xU::Float64, a::Float64)
    xU <= a && (return a, 0.0)
    a <= xL && (return x, 1.0)
    term = max(0.0, (x - a)/(xU - a))
    val = a + (xU - a)*term^MC_DIFF_MU1
    dval = MC_DIFF_MU1T*term^MC_DIFF_MU
    return val, dval
end
@inline function cv_max_ns(x::Float64, xL::Float64, xU::Float64, a::Float64)
    (x <= a) && (return a, 0.0)
    return x, 1.0
end

@inline cc_max(x::Float64, xL::Float64, xU::Float64, a::Float64) = dline_seg(max, deriv_max, x, xL, xU, a)
@inline function max_kernel(x::MC{N, Diff}, c::Float64, z::Interval{Float64}) where N
    midcv, cv_id = mid3(x.cc, x.cv, x.Intv.lo)
    midcc, cc_id = mid3(x.cc, x.cv, x.Intv.hi)
    cv, dcv = cv_max(midcv, x.Intv.lo, x.Intv.hi, c)
    cc, dcc = cc_max(midcc, x.Intv.lo, x.Intv.hi, c)
    gcc1,gdcc1 = cc_max(x.cv, x.Intv.lo, x.Intv.hi, c)
    gcv1,gdcv1 = cv_max(x.cv, x.Intv.lo, x.Intv.hi, c)
    gcc2,gdcc2 = cc_max(x.cc, x.Intv.lo, x.Intv.hi, c)
    gcv2,gdcv2 = cv_max(x.cc, x.Intv.lo, x.Intv.hi, c)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N, Diff}(cv, cc, z, cv_grad, cc_grad, x.cnst)
end
@inline function max_kernel(x::MC{N, T}, c::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    midcv, cv_id = mid3(x.cc, x.cv, x.Intv.lo)
    midcc, cc_id = mid3(x.cc, x.cv, x.Intv.hi)
    cv, dcv = cv_max_ns(midcv, x.Intv.lo, x.Intv.hi, c)
    cc, dcc = cc_max(midcc, x.Intv.lo, x.Intv.hi, c)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, z, cv_grad, cc_grad, x.cnst)
end
@inline max(x::MC, c::Float64) = max_kernel(x, c, max(x.Intv, c))
@inline min_kernel(x::MC, c::Float64, z::Interval{Float64}) = -max(-x,-c)
@inline min(x::MC, c::Float64) = -max(-x,-c)

# Multivariate min/max operators
@inline function psil_max(x::Float64, y::Float64, lambda::Interval{Float64}, nu::Interval{Float64}, f1::MC{N, Diff}, f2::MC{N, Diff}) where N
    if nu.hi <= lambda.lo
        val = x
    elseif lambda.hi <= nu.lo
        val = y
    elseif nu.lo <= lambda.lo && lambda.lo < nu.hi
        val = x + (nu.hi - lambda.lo)*max(0.0, (y - x)/(nu.hi - lambda.lo))^MC_DIFF_MU1
    else
        val = y + (lambda.hi - nu.lo)*max(0.0, (x - y)/(lambda.hi - nu.lo))^MC_DIFF_MU1
    end
    if nu.hi <= lambda.lo
        grad_val = f1.cv_grad
    elseif lambda.hi <= nu.lo
        grad_val = f1.cv_grad
    else
        grad_val = max(0.0, psil_max_dx(x, y, lambda, nu))*f1.cv_grad +
                   min(0.0, psil_max_dx(x, y, lambda, nu))*f1.cc_grad +
                   max(0.0, psil_max_dy(x, y, lambda, nu))*f2.cv_grad +
                   min(0.0, psil_max_dy(x, y, lambda, nu))*f2.cc_grad
    end
    return val, grad_val
end

@inline function psil_max_dx(x::Float64, y::Float64, lambda::Interval{Float64}, nu::Interval{Float64})
    if in(lambda.lo, nu)
        return 1.0 - MC_DIFF_MU1*max(0.0, (y - x)/(nu.hi - lambda.lo))^MC_DIFF_MU
    end
    return MC_DIFF_MU1*max(0.0, (x - y)/(lambda.hi - nu.lo))^MC_DIFF_MU
end
@inline function psil_max_dy(x::Float64, y::Float64, lambda::Interval{Float64}, nu::Interval{Float64})
    if in(lambda.lo, nu)
        return MC_DIFF_MU1*max(0.0, (y - x)/(nu.hi - lambda.lo))^MC_DIFF_MU
    end
    return 1.0 - MC_DIFF_MU1*max(0.0, (x - y)/(lambda.hi - nu.lo))^MC_DIFF_MU
end
@inline function psir_max(x::Float64, y::Float64, xgrad::SVector{N,Float64}, ygrad::SVector{N,Float64},
                          lambda::Interval{Float64}, nu::Interval{Float64}) where N

    (nu.hi <= lambda.lo) && (return (x, xgrad))
    (lambda.hi <= nu.lo) && (return (y, ygrad))

    maxUU = max(lambda.hi, nu.hi)
    maxLU = max(lambda.lo, nu.hi)
    maxUL = max(lambda.hi, nu.lo)
    maxLL = max(lambda.lo, nu.lo)

    diamx = diam(lambda)
    diamy = diam(nu)
    delta = maxLL + maxUU - maxLU - maxUL
    theta_arg = max(0.0, (lambda.hi - x)diamx - (y - nu.lo)/diamy)
    thetar = delta*theta_arg^MC_DIFF_MU1

    coeff1 = maxUU - maxLU
    coeff2 = maxUU - maxUL
    val = maxUU - coeff1*(lambda.hi-x)/diamx - coeff2*(nu.hi - y)/diamy

    shared_term = -MC_DIFF_MU1T*delta*theta_arg^MC_DIFF_MU
    grad_val = ((coeff1 + shared_term)/diamx)*xgrad + ((coeff2 + shared_term)/diamy)*ygrad
    return val, grad_val
end
@inline function max_kernel(x::MC{N, Diff}, y::MC{N, Diff}, z::Interval{Float64}) where N

    if (y.Intv.hi <= x.Intv.lo) || (x.Intv.hi <= y.Intv.lo)
        cv, cv_grad = psil_max(x.cv, y.cv, x.Intv, y.Intv, x, y)
    elseif (y.Intv.lo <= x.Intv.lo) & (x.Intv.lo < y.Intv.hi)
        temp = mid3v(x.cv, x.cc, y.cv - (y.Intv.hi - x.Intv.lo)*MC_DIFF_DIV)
        cv, cv_grad = psil_max(temp, y.cv, x.Intv, y.Intv, x, y)
    else
        temp = mid3v(y.cv, y.cc, x.cv - (x.Intv.hi - y.Intv.lo)*MC_DIFF_DIV)
        cv, cv_grad = psil_max(x.cv, temp, x.Intv, y.Intv, x, y)
    end

    cc, cc_grad = psir_max(x.cc, y.cc, x.cc_grad, y.cv_grad, x.Intv, y.Intv)
    return MC{N, Diff}(cv, cc, z, cv_grad, cc_grad, (x.cnst && y.cnst))
end

@inline function max_kernel(x::MC{N, MV}, y::MC{N, MV}, z::Interval{Float64}) where N
    if x.Intv.hi <= y.Intv.lo
        cc = y.cc
        cc_grad = y.cnst ? zero(SVector{N,Float64}) : y.cc_grad
    elseif x.Intv.lo >= y.Intv.hi
        cc = x.cc
        cc_grad = x.cnst ? zero(SVector{N,Float64}) : x.cc_grad
    else
        maxLU = max(x.Intv.lo, y.Intv.hi)
        maxUL = max(x.Intv.hi, y.Intv.lo)

        m1a = isthin(x) ? 0.0 : (x.cc - lo(x))/diam(x)
        m1b = isthin(y) ? 0.0 : (y.cc - lo(y))/diam(y)
        maxLL = max(x.Intv.lo, y.Intv.lo)
        g1cc = maxLL + m1a*(maxUL - maxLL) + m1b*(maxLU - maxLL)

        m2a = isthin(x) ? 0.0 : (x.cc - hi(x))/diam(x)
        m2b = isthin(y) ? 0.0 : (y.cc - hi(y))/diam(y)
        maxUU = max(x.Intv.hi, y.Intv.hi)
        g2cc = maxUU + m2a*(maxUU - maxLU) + m2b*(maxUU - maxUL)

        if g1cc < g2cc
            cc = g1cc
            cc_grad = (maxUL - maxLL)*x.cc_grad + (maxLU - maxLL)*y.cc_grad
        else
            cc = g2cc
            cc_grad = (maxUU - maxLU)*x.cc_grad + (maxUU - maxUL)*y.cc_grad
        end
    end
    cv = max(x.cv, y.cv)
    cv_grad = (x.cv > y.cv) ? (x.cnst ? zero(SVector{N,Float64}) : x.cv_grad) :
                              (y.cnst ? zero(SVector{N,Float64}) : y.cv_grad)
    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    return MC{N, MV}(cv, cc, z, cv_grad, cc_grad, y.cnst ? x.cnst : (x.cnst ? y.cnst : (x.cnst || y.cnst)))
end

@inline function max_kernel(x::MC{N, NS}, y::MC{N, NS}, z::Interval{Float64}) where N
    if x.Intv.hi <= y.Intv.lo
        cc = y.cc
        cc_grad = y.cnst ? zero(SVector{N,Float64}) : y.cc_grad
    elseif x.Intv.lo >= y.Intv.hi
        cc = x.cc
        cc_grad = x.cnst ? zero(SVector{N,Float64}) : x.cc_grad
    else
        ccMC = 0.5*(x + y + abs(x - y))
        cc = ccMC.cc
        cc_grad = ccMC.cc_grad
    end
    cv = max(x.cv, y.cv)
    cv_grad = (x.cv > y.cv) ? (x.cnst ? zero(SVector{N,Float64}) : x.cv_grad) :
                              (y.cnst ? zero(SVector{N,Float64}) : y.cv_grad)
    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    return MC{N, NS}(cv, cc, z, cv_grad, cc_grad, y.cnst ? x.cnst : (x.cnst ? y.cnst : (x.cnst || y.cnst)))
end

@inline function min_kernel(x::MC{N, NS}, y::MC{N, NS}, z::Interval{Float64}) where N
    if x.Intv.hi <= y.Intv.lo
        cv = x.cv
        cv_grad = x.cnst ? zero(SVector{N,Float64}) : x.cv_grad
    elseif x.Intv.lo >= y.Intv.hi
        cv = y.cv
        cv_grad = y.cnst ? zero(SVector{N,Float64}) : y.cv_grad
    else
        cvMC = 0.5*(x + y - abs(x - y))
        cv = cvMC.cv
        cv_grad = cvMC.cv_grad
    end
    cc = min(x.cc, y.cc)
    cc_grad = (x.cv > y.cv) ? (x.cnst ? zero(SVector{N,Float64}) : x.cv_grad) :
                              (y.cnst ? zero(SVector{N,Float64}) : y.cv_grad)
    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    return MC{N, NS}(cv, cc, z, cv_grad, cc_grad, y.cnst ? x.cnst : (x.cnst ? y.cnst : (x.cnst || y.cnst)))
end

@inline max(x::MC, y::MC) = max_kernel(x, y, max(x.Intv, y.Intv))
@inline min_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{MV, Diff}} = -max(-x, -y)
@inline min(x::MC, y::MC) = min_kernel(x, y, min(x.Intv, y.Intv))

@inline function union(x::MC{N,T}, y::MC{N,T}) where {N, T <: RelaxTag}
    cv_MC = min(x, y)
    cc_MC = max(x, y)
    return MC{N, NS}(cv_MC.cv, cc_MC.cc, Interval(cv.Intv.lo, cc.Intv.hi),
                     cv_MC.cv_grad, cc_MC.cc_grad, x.cnst && y.cnst)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Documenter, Pitaya

makedocs(sitename="Pitaya",
         modules = [Pitaya],
         authors="Nick Haupka",
         pages=[
            "API reference" => "index.md"])

deploydocs(;
  repo="github.com/naustica/Pitaya.git",
  devbranch="main",
  devurl="latest",
)

ext: jl

max_stars_repo_name: naustica/Pitaya

max_issues_repo_name: naustica/Pitaya

max_forks_repo_name: naustica/Pitaya


text:
using Documenter
using POMDPTesting

makedocs(
    format =:html,
    sitename = "POMDPTesting.jl"
)

deploydocs(
    repo = "github.com/JuliaPOMDP/POMDPTesting.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
#!/usr/bin/env julia

#import PlmDCA
using PlmDCA
using LinearAlgebra
using DelimitedFiles
using ArgParse
using NPZ

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--input", "-i"
            help = "Input alignment"
            required = true
        "--output", "-o"
            help = "Output root"
            required = true
        "--theta", "-t"
            help = "Reweighting value, default=0.2 (reweighting at 80% sequence identity)"
            default = 0.2
            arg_type = Float64
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    input = parsed_args["input"]
    output = parsed_args["output"]
    theta = parsed_args["theta"]

    a = plmdca(input, verbose = false, theta = theta)

    h = a.htensor
    J = a.Jtensor
    scores = a.score

    output_scores = output*".scores"
    output_stream = open(output_scores, "w")
    for pred in scores
        println(output_stream, string(pred[1]) * "\t" * string(pred[2]) * "\t" * string(pred[3]))
    end
    close(output_stream)

    # J is too big! Reduce is size
    # switch to matrix form of J
    q,N = size(h)
    J_m = zeros(q*N,q*N)
    for i=1:N
    	for j = 1:N
    		J_m[(1+q*(i-1)):(q+q*(i-1)),(1+q*(j-1)):(q+q*(j-1))] = J[:,:,i,j]
    	end
    end

    # Remove lower_triangular_part (use the fact that J_m = J_m')
    J_m -= LowerTriangular(J_m)

    # Reduce precision to 8 digits
    J_m = round.(J_m,digits=8)

    # Write a out files
    J_name = output*".J"
    writedlm(J_name, J_m)

    h_name = output*".h"
    writedlm(h_name,h)

    # Save also in NPZ format
#     name_h = h_name*".npz"
#     name_J = J_name*".npz"
#
#     npzwrite(name_h, h)
#     npzwrite(name_J, J)
#
#     name_h_compressed = name_h*".tar.gz"
#     name_J_compressed = name_J*".tar.gz"
#     run(`tar -czvf $name_h_compressed $name_h`)
#     run(`tar -czvf $name_J_compressed $name_J`)
#     run(`rm $name_h`)
#     run(`rm $name_J`)
end

main()

max_stars_repo_name: juan-rodriguez-rivas/covmut

max_issues_repo_name: juan-rodriguez-rivas/covmut

max_forks_repo_name: juan-rodriguez-rivas/covmut


text:
# find hyperparameter values using gird search

include(pwd()*"/dlmodel.jl")

using Distributions
using ProgressMeter

N = 1000

dist = Uniform(0, 0.1)

rand(dist, 3)

hyperparameters = zeros(N,3)
model_fit_regularization = zeros(N,3+3+1+1)
model_fit_train = zeros(N,2000)
model_fit_val = zeros(N,2000)

# find hyperparameters using random search
@time begin
@showprogress for i in 1:N

    Random.seed!(1234) # set random numbers

    lambda = rand(dist,3)

    hyperparameters[i,:] = lambda

    w = Any[ xavier(n_hidden_1,n_input), zeros(n_hidden_1,1),
             xavier(n_hidden_2,n_hidden_1), zeros(n_hidden_2,1),
             xavier(n_out,n_hidden_2), zeros(n_out,1)]

    optim = optimizers(w, Adam)

    dropout_percentage = 0
    loss_train_vec, loss_val_vec = train(2000, dropout_percentage, lambda, w, optim)

    model_fit_train[i,:] = loss_train_vec
    model_fit_val[i,:] = loss_val_vec

    # calc predictions
    y_pred_train_proc = predict(w, x_train, 0)
    y_pred_val_proc = predict(w, x_val, 0)

    y_pred_train = class_nll(y_pred_train_proc,[1,2])
    y_pred_val = class_nll(y_pred_val_proc,[1,2])

    # loss
    model_fit_regularization[i,end-1] = loss(w,x_train, y_train, 0, lambda)
    model_fit_regularization[i,end] = loss(w,x_val, y_val, 0, lambda)

    # classification results
    model_fit_regularization[i,1:3] = classification_results(y_pred_train, y_train, [1,2])
    model_fit_regularization[i,4:6] = classification_results(y_pred_val, y_val, [1,2])

end
end

# add run w.o regularization

w = Any[ xavier(n_hidden_1,n_input), zeros(n_hidden_1,1),
         xavier(n_hidden_2,n_hidden_1), zeros(n_hidden_2,1),
         xavier(n_out,n_hidden_2), zeros(n_out,1)]

optim = optimizers(w, Adam)

# training
dropout_percentage = 0.
lambda = [0.0,0.0,0.0]
@time loss_train_vec, loss_val_vec = train(2000, dropout_percentage, lambda, w, optim)

# calc predictions
y_pred_train_proc = predict(w, x_train, 0)
y_pred_val_proc = predict(w, x_val, 0)

y_pred_train = class_nll(y_pred_train_proc,[1,2])
y_pred_val = class_nll(y_pred_val_proc,[1,2])

# loss
training_loss_no_reg = loss(w,x_train, y_train, 0, lambda)
val_loss_no_reg = loss(w,x_val, y_val, 0, lambda)

# classification results
training_classification_res_no_reg = classification_results(y_pred_train, y_train, [1,2])
val_classification_res_no_reg = classification_results(y_pred_val, y_val, [1,2])

# find best hyperparameter settings
(~, idx_best_training_loss) = findmin(model_fit_regularization[:,end-1])
(~, idx_best_val_loss) = findmin(model_fit_regularization[:,end])

(~, idx_best_training_tp) = findmax(model_fit_regularization[:,2])
(~, idx_best_val_tp) = findmax(model_fit_regularization[:,5])

(~, idx_best_training_tn) = findmax(model_fit_regularization[:,3])
(~, idx_best_val_tn) = findmax(model_fit_regularization[:,6])

(~, idx_best_training_acc) = findmax(model_fit_regularization[:,1])
(~, idx_best_val_acc) = findmax(model_fit_regularization[:,4])

@printf "Best lambda_1 %.10f\n" lambda_1[1]
@printf "Best lambda_2 %.10f\n" lambda_2[4]
@printf "Best lambda_3 %.10f\n" lambda_3[1]

max_stars_repo_name: SamuelWiqvist/dl-compute-project

max_issues_repo_name: SamuelWiqvist/dl-compute-project

max_forks_repo_name: SamuelWiqvist/dl-compute-project


text:
module NumCME
using DocStringExtensions
using ReusePatterns, SparseArrays, StaticArrays, LinearAlgebra, Symbolics
import ModelingToolkit
import Catalyst
import ForwardDiff

include("cmemodel/cmemodel.jl")
include("statespace/statespace.jl")
include("fspvector/fspvector.jl")
include("fspmatrix/fspmatrix.jl")
include("transientcme/transientcme.jl")
include("forwardsensfspmatrix/forwardsensfspmatrix.jl")
include("forwardsenscme/forwardsenscme.jl")
end


text:
# Test naive implementation of entropy calculations towards efficient calculation of entropy on a manifold

using DocStringExtensions
using KernelDensityEstimate

include(joinpath(dirname(@__FILE__), "circularEntropyUtils.jl"))

using Gadfly, Colors
using Distributions

# some test points to work with
pts = randn(100)

# p = kde!(pts)

# N = length(pts)
CV  = zeros(100)
global i_global = 0
BW= range(0.1, 1.0, length=100)
for bw in BW
    global i_global
    i_global += 1
    CV[i_global] = manifoldLooCrossValidation(pts, bw, own=true)
end

BW[findfirst(CV .== maximum(CV))]

getBW(kde!(pts))[1,1]

## do the same with

PD = zeros(100)

global i_global = 0
for bw in BW
    global i_global
    i_global += 1
    pd = kde!(pts, [bw;])
    PD[i_global] = KernelDensityEstimate.entropy(pd)
end

##

# plot(y=loo, Geom.point)

plot(
    layer(x=BW, y=CV, Geom.line, Theme(default_color=colorant"red")),
    layer(x=BW, y=-PD.+0.1, Geom.line, Theme(default_color=colorant"blue"))
)

plot(
    layer(x=BW, y=CV+PD, Geom.line, Theme(default_color=colorant"red")),
)

plot(rbf, -5.0, 5.0)

PP = kde!([0.0;], [1.0])

PP = kde!(pts, [1.0])

gg = (x)->PP([x;])[1]

plot(gg, -5.0, 5.0)

0

## characterize computation times

using BenchmarkTools

# function rbftimetest()

val1 = [0.0;]
val2 = [0.0;]

@btime rbf!(val1, 0.9, 0.0, 1.0);

nn = Normal()

# gg = (x::Float64) -> pdf(nn, x)  # even slower with additional allocation

@btime val2[1] = pdf(nn, 0.9);

    # return val1[1], val2[1]
# end
# rbftimetest()

# nn = Normal()
#
# plot(
# layer((x)->kde!(pts)([x;])[1], -5,5),
# layer((x)->pdf(nn,x), -5,5)
# )

## optimize the KDE bandwidth

using Optim

# initial testing values
lower = 0.001
upper = 10.0

pts = 10.0.+randn(100)

minEntropyLOOCV = (bw) -> -manifoldLooCrossValidation(pts, bw)

# TODO Compare Optim.GoldenSection vs KDE.golden
@time res = optimize(minEntropyLOOCV, lower, upper, GoldenSection(), x_tol=0.001)

@time kde!(pts)

res.minimizer

getBW(kde!(pts))[1,1]

#


text:
# bucket - name of cluster's s3 bucket
# path - path to write file to in bucket
# download_path - either http(s) link to a file or a local Path indicating the source of the file
function verify_file_in_s3(bucket, path, download_path)
    if !s3_exists(Banyan.get_aws_config(), bucket, path)
        if typeof(download_path) == String &&
           (startswith(download_path, "https://") || startswith(download_path, "http://"))
            download(
                download_path,
                S3Path("s3://$(bucket)/$(path)", config = Banyan.get_aws_config()),
            )
        else  # upload local file
            cp(
                Path(download_path),
                S3Path("s3://$(bucket)/$(path)", config = Banyan.get_aws_config()),
            )
        end
    end
end

######################################
# SETUP AND CLEANUP HELPER FUNCTIONS #
######################################

function write_df_to_csv_to_s3(df, filename, filepath, bucket_name, s3path)
    CSV.write(filename, df)
    verify_file_in_s3(bucket_name, s3path, filepath)
end

function write_df_to_parquet_to_s3(df, filename, filepath, bucket_name, s3path)
    Parquet.write_parquet(filename, df)
    verify_file_in_s3(bucket_name, s3path, filepath)
end

function write_df_to_arrow_to_s3(df, filename, filepath, bucket_name, s3path)
    Arrow.write(filename, df)
    verify_file_in_s3(bucket_name, s3path, filepath)
end

function get_local_path_tripdata(s3_path)
    filename = s3_path.segments[end]
    if occursin("csv", filename)
        return "tripdata.csv"
    elseif occursin("parquet", filename)
        return "tripdata.parquet"
    elseif occursin("arrow", filename)
        return "tripdata.arrow"
    end
end

function setup_basic_tests(bucket_name=get_cluster_s3_bucket_name())
    iris_s3_paths = [
        "iris_large.csv",
        "iris_large.parquet",
        "iris_large.arrow",
        "iris_large_dir.csv/",
        "iris_species_info.csv",
        "iris_species_info.parquet",
        "iris_species_info.arrow",
    ]
    bucket_contents = s3_list_keys(Banyan.get_aws_config(), bucket_name)
    to_be_downloaded = [
        iris_s3_path for iris_s3_path in iris_s3_paths if
        # TODO: Use the following when AWSS3.jl supports folders
        # !s3_exists(Banyan.get_aws_config(), bucket_name, iris_s3_path)
        !(iris_s3_path in bucket_contents)
    ]
    if !isempty(to_be_downloaded)
        @info "Downloading $to_be_downloaded"
        iris_download_path = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris.csv"
        iris_species_info_download_path = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris_species_info.csv"
        iris_local_path = download(iris_download_path)
        iris_species_info_local_path = download(iris_species_info_download_path)
        df = CSV.read(iris_local_path, DataFrames.DataFrame, stringtype=String)
        df_s = CSV.read(iris_species_info_local_path, DataFrames.DataFrame, stringtype=String)
        # Duplicate df six times and change the species names
        species_list = df[:, :species]
        df = reduce(vcat, [df, df, df, df, df, df])
        for i = 4:18
            append!(species_list, Base.fill("species_$(i)", 50))
        end
        df[:, :species] = species_list
        write_df_to_csv_to_s3(
            df,
            "iris_large.csv",
            p"iris_large.csv",
            bucket_name,
            "iris_large.csv",
        )
        write_df_to_parquet_to_s3(
            df,
            "iris_large.parquet",
            p"iris_large.parquet",
            bucket_name,
            "iris_large.parquet",
        )
        write_df_to_arrow_to_s3(
            df,
            "iris_large.arrow",
            p"iris_large.arrow",
            bucket_name,
            "iris_large.arrow",
        )

        # Write to dir
        df_shuffle = df[shuffle(1:nrow(df)), :]
        chunk_size = 100
        for i = 1:9
            write_df_to_csv_to_s3(
                df_shuffle[((i-1)*chunk_size+1):i*chunk_size, :],
                "iris_large_chunk.csv",
                p"iris_large_chunk.csv",
                bucket_name,
                "iris_large_dir.csv/iris_large_chunk$(i).csv",
            )
        end

        write_df_to_csv_to_s3(
            df_s,
            "iris_species_info.csv",
            p"iris_species_info.csv",
            bucket_name,
            "iris_species_info.csv",
        )
        write_df_to_parquet_to_s3(
            df_s,
            "iris_species_info.parquet",
            p"iris_species_info.parquet",
            bucket_name,
            "iris_species_info.parquet",
        )
        write_df_to_arrow_to_s3(
            df_s,
            "iris_species_info.arrow",
            p"iris_species_info.arrow",
            bucket_name,
            "iris_species_info.arrow",
        )
    end
end

function setup_empty_tests(bucket_name=get_cluster_s3_bucket_name())
    # Write empty dataframe
    empty_df = DataFrames.DataFrame()
    if !ispath(S3Path("s3://$bucket_name/empty_df.csv", config = Banyan.get_aws_config()))
        write_df_to_csv_to_s3(
            empty_df,
            "empty_df.csv",
            p"empty_df.csv",
            bucket_name,
            "empty_df.csv",
        )
    end
    if !ispath(S3Path("s3://$bucket_name/empty_df.arrow", config = Banyan.get_aws_config()))
        write_df_to_arrow_to_s3(
            empty_df,
            "empty_df.arrow",
            p"empty_df.arrow",
            bucket_name,
            "empty_df.arrow",
        )
    end

    # Write empty dataframe with two columns
    empty_df2 = DataFrames.DataFrame(x = [], y = [])
    if !ispath(S3Path("s3://$bucket_name/empty_df2.csv", config = Banyan.get_aws_config()))
        write_df_to_csv_to_s3(
            empty_df2,
            "empty_df2.csv",
            p"empty_df2.csv",
            bucket_name,
            "empty_df2.csv",
        )
    end
    if !ispath(S3Path("s3://$bucket_name/empty_df2.arrow", config = Banyan.get_aws_config()))
        write_df_to_arrow_to_s3(
            empty_df2,
            "empty_df2.arrow",
            p"empty_df2.arrow",
            bucket_name,
            "empty_df2.arrow",
        )
    end
end

global n_repeats = 10

function setup_stress_tests(bucket_name=get_cluster_s3_bucket_name())
    # Copy n_repeats of each of the four files into S3.
    global n_repeats
    for month in ["01", "02", "03", "04"]
        println("In setup_stress_tests on month=$month")
        # Get source download path and the list of dst S3 paths
        download_path = "https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2012-$(month).csv"
        dst_s3_paths = []
        dst_s3_paths_missing = []
        for filetype in ["csv", "parquet", "arrow"]
            for ncopy = 1:n_repeats
                dst_path = "s3://$(bucket_name)/tripdata_large_$(filetype).$(filetype)/tripdata_$(month)_copy$(ncopy).$(filetype)"
                dst_s3_path = S3Path(dst_path, config = Banyan.get_aws_config())
                append!(dst_s3_paths, dst_s3_path)
                if !isfile(dst_s3_path)
                    append!(dst_s3_paths_missing, dst_s3_path)
                end
            end
        end
        # If at least one file doesn't exist in s3, we need to download and write to s3
        if length(dst_s3_paths_missing) > 0
            local_path = download(download_path)
            println("Downloaded $local_path")
            df = CSV.read(local_path, DataFrames.DataFrame)
            println("Read $local_path into memory")

            # Write csv to disk if there is at least one missing csv file
            if any(p -> occursin("csv", p), dst_s3_paths_missing)
                CSV.write("tripdata.csv", df)
            end
            # Write parquet to disk if there is at least one missing parquet file
            if any(p -> occursin("parquet", p), dst_s3_paths_missing)
                Parquet.write_parquet("tripdata.parquet", df)
            end
            # Write arrow to disk if there is at least one missing arrow file
            if any(p -> occursin("arrow", p), dst_s3_paths_missing)
                Arrow.write("tripdata.arrow", df)
            end

            # Loop over missing files and upload to s3
            for s3_path in dst_s3_paths_missing
                cp(
                    Path(get_local_path_tripdata(s3_path)),
                    s3_path,
                    config = Banyan.get_aws_config(),
                )
            end
        end
    end
end

function cleanup_tests(bucket_name=get_cluster_s3_bucket_name())
    # Delete all temporary test files that are prepended with "test-tmp__"
    @show bucket_name
    for p in s3_list_keys(Banyan.get_aws_config(), bucket_name)
        if contains(string(p), "test-tmp_")
            # s3_path = S3Path("s3://$bucket_name/$p", config = Banyan.get_aws_config())
            rm(S3Path("s3://$bucket_name/$p", config = Banyan.get_aws_config()), recursive=true)
        end
    end
end

##############################
# HELPER FUNCTIONS FOR TESTS #
##############################

function read_file(path)
    if endswith(path, ".csv")
        return BanyanDataFrames.read_csv(path)
    elseif endswith(path, ".parquet")
        return BanyanDataFrames.read_parquet(path)
    elseif endswith(path, ".arrow")
        return BanyanDataFrames.read_arrow(path)
    else
        error("Invalid file format")
    end
end

function write_file(path, df)
    if endswith(path, ".csv")
        BanyanDataFrames.write_csv(df, path)
    elseif endswith(path, ".parquet")
        BanyanDataFrames.write_parquet(df, path)
    elseif endswith(path, ".arrow")
        BanyanDataFrames.write_arrow(df, path)
    else
        error("Invalid file format")
    end
end

function get_save_path(bucket, df_name, path)
    return "s3://$(bucket)/test-tmp_$(df_name)_$(split(path, "/")[end])"
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
# This tests the behaviour of multiple successive solves with problem data updates

using Test, LinearAlgebra, SparseArrays
using COSMO

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

T = Float64

@testset "Model modifications" begin

# for T in UnitTestFloats
    # calling optimize!() twice should yield the same results (with less iterations)
    q = [1; 1.];
    P = sparse([4. 1; 1 2]);
    A = [1. 1; 1 0; 0 1];
    l = [1.; 0; 0];
    u = [1; 0.7; 0.7];
    Aa = [-A; A]
    ba = [u; -l]
    constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);
    model = COSMO.Model();
    assemble!(model, P, q, constraint1, settings= COSMO.Settings(check_termination = 1));
    @test model.states.IS_ASSEMBLED = true
    res = COSMO.optimize!(model);
    res2 = COSMO.optimize!(model);
    @test abs(res.obj_val - res2.obj_val) <= 1e-3
    @test res2.iter <= res.iter

    # updating q should solve correctly without reassembling the KKT matrix
    model = COSMO.Model();
    assemble!(model, P, q, constraint1);
    res = COSMO.optimize!(model);
    # update q and compare to known solution
    q_new = [2.; 3]
    COSMO.update!(model, q = q_new)
    res = COSMO.optimize!(model);
    @test abs(res.obj_val - 3.5) < 1e-3
    @test norm(res.x - [0.5; 0.5], 2) < 1e-3

    # updating b should solve correctly without having to reassemble the KKT matrix, i.e. same factorisation time
    # min   x_1 + x_2
    # s.t.  x_1 >= 2 --> x* = [2, 3]
    #       x_2 >= 3
    q = [1; 1.];
    P = spzeros(2, 2);
    A = spdiagm(0 => ones(2))
    b = [-2.; -3.]
    model = COSMO.Model();
    assemble!(model, P, q, COSMO.Constraint(A, b, COSMO.Nonnegatives), settings = COSMO.Settings(verbose_timing = true, check_termination = 20, ));
    res = COSMO.optimize!(model);

    # update b to give: x_1 >= 0, x_2 >= -1. --> x*=[0, -1]
    @test model.states.KKT_FACTORED
    COSMO.update!(model, b = [0.; 1.])
    res2 = COSMO.optimize!(model);
    @test norm(res2.x - [0; -1.], 2) < 1e-4
    @test res.times.init_factor_time == res2.times.init_factor_time

    # warm-starting of x,y should yield correct s = b - Ax and lead to a fast second solve

    # test full emptying of the model
    COSMO.empty_model!(model)
    @test model.states.IS_OPTIMIZED == false
    @test model.states.IS_ASSEMBLED == false
    @test model.states.KKT_FACTORED == false
    @test model.states.IS_SCALED == false

    assemble!(model, P, q, COSMO.Constraint(A, b, COSMO.Nonnegatives), settings = COSMO.Settings(verbose_timing = true, check_termination = 20));
    @test model.states.IS_ASSEMBLED
    res3 = COSMO.optimize!(model);
    @test abs(res.obj_val - res3.obj_val) < 1e-3
    @test model.states.IS_OPTIMIZED
    @test model.states.KKT_FACTORED
    @test model.states.IS_SCALED
end

# end


text:
using RealSense

include("rs-utils.jl")
include("gl-utils.jl")

# initialize GLFW window
glfwWidth = 1280
glfwHeight = 720
window = create_window_context(glfwWidth, glfwHeight)

# prepare drawing via OpenGL
depth_window = GLfloat[-1.0,-1.0,   0.0,-1.0,   0.0, 1.0,
                        0.0, 1.0,  -1.0, 1.0,  -1.0,-1.0]
rgb_window = GLfloat[0.0,-1.0,   1.0,-1.0,   1.0, 1.0,
                     1.0, 1.0,   0.0, 1.0,   0.0,-1.0]
vertex_texcoords = GLfloat[1.0,1.0,  0.0,1.0,  0.0,0.0,
                           0.0,0.0,  1.0,0.0,  1.0,1.0]
# create Vertex Buffer Objects
depthVBO = Ref{GLuint}(0)
glGenBuffers(1, depthVBO)
glBindBuffer(GL_ARRAY_BUFFER, depthVBO[])
glBufferData(GL_ARRAY_BUFFER, sizeof(depth_window), depth_window, GL_STATIC_DRAW)
rgbVBO = Ref{GLuint}(0)
glGenBuffers(1, rgbVBO)
glBindBuffer(GL_ARRAY_BUFFER, rgbVBO[])
glBufferData(GL_ARRAY_BUFFER, sizeof(rgb_window), rgb_window, GL_STATIC_DRAW)
texcoordsVBO = Ref{GLuint}(0)
glGenBuffers(1, texcoordsVBO)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_texcoords), vertex_texcoords, GL_STATIC_DRAW)
# create Vertex Array Objects
depthVAO = Ref{GLuint}(0)
glGenVertexArrays(1, depthVAO)
glBindVertexArray(depthVAO[])
glBindBuffer(GL_ARRAY_BUFFER, depthVBO[])
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)
rgbVAO = Ref{GLuint}(0)
glGenVertexArrays(1, rgbVAO)
glBindVertexArray(rgbVAO[])
glBindBuffer(GL_ARRAY_BUFFER, rgbVBO[])
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)
# create shader program
vertSource= joinpath(@__DIR__, "capture.vert")
fragSource = joinpath(@__DIR__, "capture.frag")
shaderProgramID = createprogram(vertSource, fragSource)
# generate texture handle
depthTex = Ref{GLuint}(0)
glGenTextures(1, depthTex)
rgbTex = Ref{GLuint}(0)
glGenTextures(1, rgbTex)

# initialize RealSense camera
err = Ref{Ptr{rs2_error}}(0)
ctx = rs2_create_context(RS2_API_VERSION, err)
checkerror(err)
device_list = rs2_query_devices(ctx, err)
checkerror(err)
dev_count = rs2_get_device_count(device_list, err)
checkerror(err)
@info "There are $dev_count connected RealSense devices."
_pipeline = rs2_create_pipeline(ctx, err)
checkerror(err)
# start the pipeline streaming
@static if is_apple()
    # we sleep 1s on MacOS, otherwise it will hit issue #1586:
    # https://github.com/IntelRealSense/librealsense/issues/1586
    sleep(1)
end
pipeline_profile = rs2_pipeline_start(_pipeline, err)
@assert err[] == C_NULL "Failed to start pipeline!"

# create Depth-Colorizer processing block that can be used to quickly visualize the depth data
color_map = rs2_create_colorizer(err)
checkerror(err)
# create target output queue for the results
frame_queue = rs2_create_frame_queue(1, err)
checkerror(err)
# set the output of the processing block to be the queue
rs2_start_processing_queue(color_map, frame_queue, err)
checkerror(err)

# OpenGL rendering loop
while !GLFW.WindowShouldClose(window)
    updatefps!(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, glfwWidth, glfwHeight)
    glUseProgram(shaderProgramID)

    # wait next frame set
    frames = rs2_pipeline_wait_for_frames(_pipeline, 5000, err)
    checkerror(err)
    # get the number of frames embedded within the composite frame
    num_of_frames = rs2_embedded_frames_count(frames, err)
    checkerror(err)
    for i = 0:num_of_frames-1
        # extract frame
        frame = rs2_extract_frame(frames, i, err)
        checkerror(err)
        # extract stream profile
        profile = rs2_get_frame_stream_profile(frame, err)
        checkerror(err)
        stream, format = Ref{rs2_stream}(0), Ref{rs2_format}(0)
        index, unique_id, framerate = Ref{Cint}(0), Ref{Cint}(0), Ref{Cint}(0)
        rs2_get_stream_profile_data(profile, stream, format, index, unique_id, framerate, err)
        checkerror(err)
        # select stream
        if stream[] == RS2_STREAM_DEPTH
            # frame reference is being "moved" into processing, so unless you
            # don't need frame anymore you must add an extra rs2_frame_add_ref
            rs2_frame_add_ref(frame, err)
            checkerror(err)
            # invoke processing
            rs2_process_frame(color_map, frame, err)
            checkerror(err)
            colorized_frame = rs2_wait_for_frame(frame_queue, 5000, err)
            checkerror(err)
            # get data
            depth_data = Ptr{UInt8}(rs2_get_frame_data(colorized_frame, err))
            checkerror(err)
            width = rs2_get_frame_width(frame, err)
            height = rs2_get_frame_height(frame, err)
            glBindTexture(GL_TEXTURE_2D, depthTex[])
            glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, depth_data)
            mipmap()
            glBindVertexArray(depthVAO[])
            glDrawArrays(GL_TRIANGLES, 0, 6)
            rs2_release_frame(colorized_frame)
        elseif stream[] == RS2_STREAM_COLOR
            rgb_data = Ptr{UInt8}(rs2_get_frame_data(frame, err))
            checkerror(err)
            width = rs2_get_frame_width(frame, err)
            height = rs2_get_frame_height(frame, err)
            glBindTexture(GL_TEXTURE_2D, rgbTex[])
            glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, rgb_data)
            mipmap()
            glBindVertexArray(rgbVAO[])
            glDrawArrays(GL_TRIANGLES, 0, 6)
        end
        # release resource
        rs2_release_frame(frame)
    end
    rs2_release_frame(frames)
    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(window)
end

# stop the pipeline streaming
rs2_pipeline_stop(_pipeline, err)
checkerror(err)

# release rs2 resources
rs2_delete_frame_queue(frame_queue)
rs2_delete_processing_block(color_map)
rs2_delete_pipeline_profile(pipeline_profile)
rs2_delete_pipeline(_pipeline)
rs2_delete_device_list(device_list)
rs2_delete_context(ctx)

# release OpenGL resources
GLFW.DestroyWindow(window)


text:
"""
Checks that the matrix is swapable under constrained permutations, i.e. one of
the two conformations called c1 and c2.
"""
function constrained_swapable(A::Array{Int64, 2})
  @assert size(A) == (2, 2)
  const c1 = [1 0; 0 1]
  const c2 = [0 1; 1 0]
  return (A == c1) | (A == c2)
end

"""
Checks that the matrix is swapable under free permutations, i.e. it has
at least one non-zero element, and at least one zero element.

This check is necessary because if not, the free swaps will collect all ones in
a few rows / columns, which results in very blocky (and therefore biased)
matrices. As a consequence, all sub-matrices that are "constrained swapable" are
"free swapable" too, but the reciprocal is not true.

Because more sub-matrices are free swapable, and because there are higher order
constraints on the acceptable swapped matrices.

Also, this is a serious lot of documentation in a function that is not exported.
"""
function free_swapable(A::Array{Int64, 2})
  @assert size(A) == (2, 2)
  return (sum(A) >= 1) & (sum(A) <= 3)
end

"""
Performs a constrained swap: the marginals are conserved.
"""
function constrained_swap(A::Array{Int64, 2})
  @assert size(A) == (2, 2)
  const c1 = [1 0; 0 1]
  const c2 = [0 1; 1 0]
  if A == c2
    return c1
  else A == c1
    return c2
  end
end

"""
Perfoms a free swap: the matrix cells are shuffled. The marginals are not
(necessarily) conserved.
"""
function free_swap(A::Array{Int64, 2})
  @assert size(A) == (2, 2)
  return reshape(A[shuffle([1, 2, 3, 4])], (2, 2))
end

"""
Performs 30000 2x2 swaps of a matrix by preserving the marginal totals (degree of the
network).
"""
function null_preserve_marginals(A::Array{Int64, 2})
  X = copy(A)
  nswaps = 30000
  dswaps = 0 # Done swaps
  tswaps = 0 # Attempted swaps
  START_TIME = time()
  while dswaps < nswaps
    tswaps += 1
    rows = StatsBase.sample(1:size(X,1), 2, replace=false)
    cols = StatsBase.sample(1:size(X,2), 2, replace=false)
    # Only if the submatrix is swapable do we actually swap it
    if constrained_swapable(X[rows,cols])
      dswaps += 1
      X[rows,cols] = constrained_swap(X[rows,cols])
    end

    if time() - START_TIME >= 120
      START_TIME = time()
    end
  end

  return X
end

"""
Performs 30000 2x2 swaps of a matrix by preserving the marginal totals of ROWS only.
"""
function null_preserve_rows_marginals(A::Array{Int64, 2})
  X = copy(A)
  nswaps = 30000
  dswaps = 0 # Done swaps
  tswaps = 0 # Attempted swaps
  START_TIME = time()
  while dswaps < nswaps
    tswaps += 1
    rows = StatsBase.sample(1:size(X,1), 2, replace=false)
    cols = StatsBase.sample(1:size(X,2), 2, replace=false)
    # Only if the submatrix is swapable do we actually swap it
    if free_swapable(X[rows,cols])
      subX = free_swap(X[rows,cols])
      Y = copy(X) # NOTE this may not be optimal because it moves potentially large objects in memory
      Y[rows,cols] = subX
      if no_empty_rows(Y) & no_empty_columns(Y) & same_row_marginals(X, Y)
        X[rows,cols] = subX
        dswaps += 1
      end
    end

    if time() - START_TIME >= 120
      START_TIME = time()
    end
  end

  return X
end

"""
Calls `null_preserve_rows_marginals` on `A'`.
"""
function null_preserve_columns_marginals(A::Array{Int64, 2})
  return null_preserve_rows_marginals(A')'
end

"""
Performs 30000 2x2 swaps of a matrix by preserving the fill only.
"""
function null_preserve_fill(A::Array{Int64, 2})
  X = copy(A)

  nswaps = 30000
  dswaps = 0 # Done swaps
  tswaps = 0 # Attempted swaps
  START_TIME = time()
  while dswaps < nswaps
    tswaps += 1
    rows = StatsBase.sample(1:size(X,1), 2, replace=false)
    cols = StatsBase.sample(1:size(X,2), 2, replace=false)
    # Only if the submatrix is swapable do we actually swap it
    if free_swapable(X[rows,cols])
      subX = free_swap(X[rows,cols])
      Y = copy(X) # NOTE this may not be optimal because it moves potentially large objects in memory
      Y[rows,cols] = subX
      if no_empty_rows(Y) & no_empty_columns(Y)
        X[rows,cols] = subX # NOTE this avoids copying B back into A
        dswaps += 1
      end
    end

    if time() - START_TIME >= 120
      START_TIME = time()
    end
  end

  return X
end


text:
import ChainRulesCore.rrule
import Yota.YotaRuleConfig

@testset "rulesets" begin

    f = reshape
    args = [(rand(3, 4),), (12,)]
    test_rrule(YotaRuleConfig(), Core._apply_iterate, iterate, f, args...; check_inferred=false)

end

# broacastable non-primitive
sin_inc(x::Number) = sin(x) + 1

@testset "generic broadcasted" begin
    # see the discussion here:
    # https://github.com/JuliaDiff/ChainRules.jl/issues/531
    for f in [sin, sin_inc]
        xs = rand(2)

        # manually get pullbacks for each element and apply them to seed 1.0
        pbs = !isnothing(rrule(f, xs[1])) ?
                [rrule(f, x)[2] for x in xs] :       # just in case
                [rrule_via_ad(YotaRuleConfig(), f, x)[2] for x in xs]
        dxs = [pbs[1](1.0)[2], pbs[2](1.0)[2]]

        # use rrule for broadcasted
        _, bcast_pb = rrule(YotaRuleConfig(), Broadcast.broadcasted, f, xs)
        dxs_bcast = bcast_pb(ones(2))[end]

        @test all(dxs .== dxs_bcast)
    end
end


text:
# This file is a part of JuliaFEM/TetGen.
# License is MIT: see https://github.com/JuliaFEM/TetGen.jl/blob/master/LICENSE.md

module TetGen
using Cxx

cxx"""#define TETLIBRARY"""
const header_dir = joinpath(Pkg.dir("TetGen"),"deps/usr/include/")
const tetgen_so = joinpath(Pkg.dir("TetGen"),"deps/usr/lib/libtet.so")
addHeaderDir(header_dir, kind=C_System)
Libdl.dlopen(tetgen_so, Libdl.RTLD_GLOBAL)
cxxinclude("tetgen.h")

in_ = @cxxnew tetgenio()
out_ = @cxxnew tetgenio()
addin_ = @cxxnew tetgenio()
bgmin_ = @cxxnew tetgenio()

end # module


text:
export ZabbixUser, get_users

struct ZabbixUser
  id::Int
  name::String
  surname::String
  alias::String
  _type::Int
  theme::String
end
ZabbixUser(dict) =
  ZabbixUser(
    parse(Int, dict["userid"]),
    dict["name"],
    dict["surname"],
    dict["alias"],
    parse(Int, dict["type"]),
    dict["theme"],
  )

function showres(u::ZabbixUser; offset=0)
  print(repeat(" ", offset), "ZabbixUser: ", u.id, " => ")
  if u.name != ""
    if u.surname != "" && u.name != u.surname
      print(u.surname, ", ")
    end
    print(u.name)
  else
    print(u.alias)
  end
  println()
end

function get_users(z)
  # TODO: getAccess
  _users = @checkerr make_request(z, "user.get")
  users = ZabbixUser[]
  for u in _users
    push!(users, ZabbixUser(u))
  end
  users
end


text:
"createHDF5: Creates or opens HDF5 file."
function createHDF5(sol::Solution, file::AbstractString; overwrite=true)
    if overwrite
        flag = "w"
        get_config(:verbosity) > 1 ? @info("Creating HDF5 file ", file) : nothing
        isfile(file) ? @warn("Overwriting existing HDF5 file.") : nothing
    else
        flag = "r+"
        get_config(:verbosity) > 1 ? @info("Opening HDF5 file ", file) : nothing
    end

    # create or open HDF5 file
    h5 = h5open(file, flag)

    return h5
end

"save_attributes: Saves common attributes of Solution to HDF5 file."
function save_attributes(sol::Solution)
    save_attributes(sol, hdf5(sol))
end

"save_attributes: Saves attributes of Deterministic Solutions to HDF5 file."
function save_attributes(sol::DeterministicSolution, h5::HDF5.File)
    # save attributes
    attributes(h5)["ntime"] = ntime(sol)
    attributes(h5)["nsave"] = nsave(sol)
    attributes(h5)["nsamples"]  = nsamples(sol)
end

"write_to_hdf5: Wrapper for saving Solution to HDF5 file."
function write_to_hdf5(solution::Solution)
    write_to_hdf5(solution, hdf5(solution), offset(solution))
end

"Creates HDF5 file, writes solution to file, and closes file."
function write_to_hdf5(solution::Solution, file::AbstractString)
    h5 = create_hdf5(solution, file)
    write_to_hdf5(solution, h5)
    close(h5)
end

"Creates HDF5 file and initialises datasets for solution object."
function create_hdf5!(solution::Solution, file::AbstractString; kwargs...)
    solution.h5 = create_hdf5(solution, file; kwargs...)
end

"Creates HDF5 file and initialises datasets for deterministic solution object."
function create_hdf5(solution::DeterministicSolution, file::AbstractString)
    # create HDF5 file and save attributes and common parameters
    h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution, h5)

    # create dataset
    init_timeteps_in_hdf5(solution, h5)
    init_solution_in_hdf5(solution, h5)

    return h5
end

function init_timeteps_in_hdf5(solution::Solution{DT,TT}, h5::HDF5.File) where {DT,TT}
    t = create_dataset(h5, "t", TT, ((solution.nt+1,), (-1,)), chunk=(1,))
    t[1] = timesteps(solution)[0]
end

function init_solution_in_hdf5(solution::SolutionODE{DT,TT,1}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    q[1] = solution.q[0]
end

function init_solution_in_hdf5(solution::SolutionODE{AT,TT,1}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    q[elaxes..., 1] = solution.q[0]
end

function init_solution_in_hdf5(solution::SolutionODE{DT,TT,2}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    q[1,:] = solution.q[0,:]
end

function init_solution_in_hdf5(solution::SolutionODE{AT,TT,2}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin,begin])
    elaxes = axes(solution.q[begin,begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.q,2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end
end

function init_solution_in_hdf5(solution::SolutionPODE{DT,TT,1}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    p = create_dataset(h5, "p", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    q[1] = solution.q[0]
    p[1] = solution.p[0]
end

function init_solution_in_hdf5(solution::SolutionPODE{AT,TT,1}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    q[elaxes..., 1] = solution.q[0]

    elsize = size(solution.p[begin])
    elaxes = axes(solution.p[begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    p[elaxes..., 1] = solution.p[0]
end

function init_solution_in_hdf5(solution::SolutionPODE{DT,TT,2}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    p = create_dataset(h5, "p", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    q[1,:] = solution.q[0,:]
    p[1,:] = solution.p[0,:]
end

function init_solution_in_hdf5(solution::SolutionPODE{AT,TT,2}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin,begin])
    elaxes = axes(solution.q[begin,begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.q,2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end

    elsize = size(solution.p[begin,begin])
    elaxes = axes(solution.p[begin,begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.p,2)
        p[elaxes..., 1, k] = solution.p[0, k]
    end
end

function init_solution_in_hdf5(solution::SolutionDAE{DT,TT,1}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    λ = create_dataset(h5, "λ", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    q[1] = solution.q[0]
    λ[1] = solution.λ[0]
end

function init_solution_in_hdf5(solution::SolutionDAE{AT,TT,1}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    q[elaxes..., 1] = solution.q[0]

    elsize = size(solution.λ[begin])
    elaxes = axes(solution.λ[begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    λ[elaxes..., 1] = solution.λ[0]
end

function init_solution_in_hdf5(solution::SolutionDAE{DT,TT,2}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    λ = create_dataset(h5, "λ", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    q[1,:] = solution.q[0,:]
    λ[1,:] = solution.λ[0,:]
end

function init_solution_in_hdf5(solution::SolutionDAE{AT,TT,2}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin,begin])
    elaxes = axes(solution.q[begin,begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.q,2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end

    elsize = size(solution.λ[begin,begin])
    elaxes = axes(solution.λ[begin,begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.λ,2)
        λ[elaxes..., 1, k] = solution.λ[0, k]
    end
end

function init_solution_in_hdf5(solution::SolutionPDAE{DT,TT,1}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    p = create_dataset(h5, "p", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    λ = create_dataset(h5, "λ", DT, ((solution.nt+1,), (-1,)), chunk=(1,))
    q[1] = solution.q[0]
    p[1] = solution.p[0]
    λ[1] = solution.λ[0]
end

function init_solution_in_hdf5(solution::SolutionPDAE{AT,TT,1}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    q[elaxes..., 1] = solution.q[0]

    elsize = size(solution.p[begin])
    elaxes = axes(solution.p[begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    p[elaxes..., 1] = solution.p[0]

    elsize = size(solution.λ[begin])
    elaxes = axes(solution.λ[begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt+1), (elsize..., -1)), chunk=(elsize..., 1))
    λ[elaxes..., 1] = solution.λ[0]
end

function init_solution_in_hdf5(solution::SolutionPDAE{DT,TT,2}, h5::HDF5.File) where {DT <: Number, TT}
    q = create_dataset(h5, "q", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    p = create_dataset(h5, "p", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    λ = create_dataset(h5, "λ", DT, ((solution.nt+1, solution.ni), (-1, solution.ni)), chunk=(1,1))
    q[1,:] = solution.q[0,:]
    p[1,:] = solution.p[0,:]
    λ[1,:] = solution.λ[0,:]
end

function init_solution_in_hdf5(solution::SolutionPDAE{AT,TT,2}, h5::HDF5.File) where {DT, AT <: Array{DT}, TT}
    elsize = size(solution.q[begin,begin])
    elaxes = axes(solution.q[begin,begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.q,2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end

    elsize = size(solution.p[begin,begin])
    elaxes = axes(solution.p[begin,begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.p,2)
        p[elaxes..., 1, k] = solution.p[0, k]
    end

    elsize = size(solution.λ[begin,begin])
    elaxes = axes(solution.λ[begin,begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt+1, solution.ni), (elsize..., -1, solution.ni)), chunk=(elsize..., 1, 1))
    for k in axes(solution.λ,2)
        λ[elaxes..., 1, k] = solution.λ[0, k]
    end
end

function copy_timeteps_to_hdf5(sol::Solution, h5::HDF5.File, j1, j2, n1, n2)
    if size(h5["t"],1) < j2
        HDF5.set_extent_dims(h5["t"], (j2,))
    end
    h5["t"][j1:j2] = timesteps(sol)[n1:n2]
end

function copy_solution_to_hdf5(solution::Union{SolutionODE{DT,TT,1},SolutionDAE{DT,TT,1}}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (j2,))
    end
    h5["q"][j1:j2] = solution.q[n1:n2]
end

function copy_solution_to_hdf5(solution::Union{SolutionODE{AT,TT,1},SolutionDAE{AT,TT,1}}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin])
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["q"][elaxes..., j] = solution.q[n]
    end
end

function copy_solution_to_hdf5(solution::Union{SolutionODE{DT,TT,2},SolutionDAE{DT,TT,2}}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["q"],1) < j2
        HDF5.set_extent_dims(h5["q"], (j2, size(h5["q"],2)))
    end
    h5["q"][j1:j2, :] = solution.q[n1:n2, :]
end

function copy_solution_to_hdf5(solution::Union{SolutionODE{AT,TT,2},SolutionDAE{AT,TT,2}}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin,begin])
    if size(h5["q"])[end-1] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[begin:end-2]..., j2, size(h5["q"])[end]))
    end
    for k in 1:nsamples(solution.q)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["q"][elaxes..., j, k] = solution.q[n, k]
        end
    end
end

function copy_solution_to_hdf5(solution::Union{SolutionPODE{DT,TT,1},SolutionPDAE{DT,TT,1}}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (j2,))
    end
    if size(h5["p"])[end] < j2
        HDF5.set_extent_dims(h5["p"], (j2,))
    end
    h5["q"][j1:j2] = solution.q[n1:n2]
    h5["p"][j1:j2] = solution.p[n1:n2]
end

function copy_solution_to_hdf5(solution::Union{SolutionPODE{AT,TT,1},SolutionPDAE{AT,TT,1}}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin])
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["q"][elaxes..., j] = solution.q[n]
    end

    elaxes = axes(solution.p[begin])
    if size(h5["p"])[end] < j2
        HDF5.set_extent_dims(h5["p"], (size(h5["p"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["p"][elaxes..., j] = solution.p[n]
    end
end

function copy_solution_to_hdf5(solution::Union{SolutionPODE{DT,TT,2},SolutionPDAE{DT,TT,2}}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["p"],1) < j2
        HDF5.set_extent_dims(h5["p"], (j2, size(h5["p"],2)))
    end
    if size(h5["p"],1) < j2
        HDF5.set_extent_dims(h5["p"], (j2, size(h5["p"],2)))
    end
    h5["q"][j1:j2, :] = solution.q[n1:n2, :]
    h5["p"][j1:j2, :] = solution.p[n1:n2, :]
end

function copy_solution_to_hdf5(solution::Union{SolutionPODE{AT,TT,2},SolutionPDAE{AT,TT,2}}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin,begin])
    if size(h5["q"])[end-1] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[begin:end-2]..., j2, size(h5["q"])[end]))
    end
    for k in 1:nsamples(solution.q)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["q"][elaxes..., j, k] = solution.q[n, k]
        end
    end

    elaxes = axes(solution.p[begin,begin])
    if size(h5["p"])[end-1] < j2
        HDF5.set_extent_dims(h5["p"], (size(h5["p"])[begin:end-2]..., j2, size(h5["p"])[end]))
    end
    for k in 1:nsamples(solution.p)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["p"][elaxes..., j, k] = solution.p[n, k]
        end
    end
end

function copy_multiplier_to_hdf5(solution::Union{SolutionDAE{DT,TT,1},SolutionPDAE{DT,TT,1}}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["λ"])[end] < j2
        HDF5.set_extent_dims(h5["λ"], (j2,))
    end
    h5["λ"][j1:j2] = solution.λ[n1:n2]
end

function copy_multiplier_to_hdf5(solution::Union{SolutionDAE{AT,TT,1},SolutionPDAE{AT,TT,1}}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.λ[begin])
    if size(h5["λ"])[end] < j2
        HDF5.set_extent_dims(h5["λ"], (size(h5["λ"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["λ"][elaxes..., j] = solution.λ[n]
    end
end

function copy_multiplier_to_hdf5(solution::Union{SolutionDAE{DT,TT,2},SolutionPDAE{DT,TT,2}}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["λ"],1) < j2
        HDF5.set_extent_dims(h5["λ"], (j2, size(h5["λ"],2)))
    end
    h5["λ"][j1:j2, :] = solution.λ[n1:n2, :]
end

function copy_multiplier_to_hdf5(solution::Union{SolutionDAE{AT,TT,2},SolutionPDAE{AT,TT,2}}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.λ[begin,begin])
    if size(h5["λ"])[end-1] < j2
        HDF5.set_extent_dims(h5["λ"], (size(h5["λ"])[begin:end-2]..., j2, size(h5["λ"])[end]))
    end
    for k in 1:nsamples(solution.λ)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["λ"][elaxes..., j, k] = solution.λ[n, k]
        end
    end
end

"Append solution to HDF5 file."
function write_to_hdf5(solution::DeterministicSolution, h5::HDF5.File=hdf5(solution), offset=offset(solution))
    # set convenience variables and compute ranges
    j1 = offset+2
    j2 = offset+1+solution.nt

    # copy data from solution to HDF5 dataset
    copy_timeteps_to_hdf5(solution, h5, j1, j2, 1, solution.nt)
    copy_solution_to_hdf5(solution, h5, j1, j2, 1, solution.nt)

    return nothing
end

"Append solution to HDF5 file."
function write_to_hdf5(solution::Union{SolutionDAE,SolutionPDAE}, h5::HDF5.File=hdf5(solution), offset=offset(solution))
    # set convenience variables and compute ranges
    j1 = offset+2
    j2 = offset+1+solution.nt

    # copy data from solution to HDF5 dataset
    copy_timeteps_to_hdf5(solution, h5, j1, j2, 1, solution.nt)
    copy_solution_to_hdf5(solution, h5, j1, j2, 1, solution.nt)
    copy_multiplier_to_hdf5(solution, h5, j1, j2, 1, solution.nt)

    return nothing
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/RCPSSPScenario.jl")
    include("../src/components/SulphateForcing.jl")

    scenario = addrcpsspscenario(m, scenario)
    sulfemit = add_comp!(m, SulphateForcing)

    sulfemit[:pse_sulphatevsbase] = scenario[:pse_sulphatevsbase]

    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)

    forcing = m[:SulphateForcing,:fs_sulphateforcing]
    forcing_compare = readpagedata(m, "test/validationdata/$valdir/fs_sulfateforcing.csv")

    @test forcing ≈ forcing_compare rtol = 1e-3
end


text:
# ==================================
# Forward approximation
# ==================================

"""
    Forward{EM, SO, SI, IT, BT} <: AbstractApproximationModel

Forward approximation model.

### Fields

- `exp`     -- exponentiation method
- `setops`  -- set operations method
- `sih`     -- symmetric interval hull
- `inv`     -- (optional, default: `false`) if `true`, assume that the state matrix
               is invertible and use its inverse in the `Φ` functions
- `backend` -- (optional, default: `nothing`) used if the algorithm needs to apply
               concrete polyhedral computations

### Algorithm

The transformations are:

- ``Φ ← \\exp(Aδ)``,
- ``Ω_0 ← CH(\\mathcal{X}_0, Φ\\mathcal{X}_0 ⊕ δU(0) ⊕ E_ψ(U(0), δ) ⊕ E^+(\\mathcal{X}_0, δ))``,
- ``V(k) ← δU(k) ⊕ E_ψ(U(k), δ)``.

Here we allow ``U`` to be a sequence of time varying non-deterministic input sets.

For the definition of the sets ``E_ψ`` and ``E^+`` see [[FRE11]](@ref).
The `Backward` method uses ``E^-``.
"""
struct Forward{EM, SO, SI, IT, BT} <: AbstractApproximationModel
    exp::EM
    setops::SO
    sih::SI
    inv::IT
    backend::BT
end

hasbackend(alg::Forward) = !isnothing(alg.backend)

# convenience constructor using symbols
function Forward(; exp=BaseExp, setops=:lazy, sih=:concrete, inv=false, backend=nothing)
    return Forward(_alias(exp), _alias(setops), Val(sih), Val(inv), backend)
end

function Base.show(io::IO, alg::Forward)
    print(io, "`Forward` approximation model with: \n")
    print(io, "    - exponentiation method: $(alg.exp) \n")
    print(io, "    - set operations method: $(alg.setops)\n")
    print(io, "    - symmetric interval hull method: $(alg.sih)\n")
    print(io, "    - invertibility assumption: $(alg.inv)")
    print(io, "    - polyhedral computations backend: $(alg.backend)")
end

Base.show(io::IO, m::MIME"text/plain", alg::Forward) = print(io, alg)

# ------------------------------------------------------------
# Forward Approximation: Homogeneous case
# ------------------------------------------------------------

# if A == |A|, then Φ can be reused in the computation of Φ₂(|A|, δ)
function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ, alg::Forward)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)

    Φ = _exp(A, δ, alg.exp)
    A_abs = _elementwise_abs(A)
    Φcache = sum(A) == abs(sum(A)) ? Φ : nothing
    P2A_abs = _Φ₂(A_abs, δ, alg.exp, alg.inv, Φcache)

    E₊ = sih(P2A_abs * sih((A * A) * X0, alg.sih), alg.sih)
    Ω0 = ConvexHull(X0, Φ * X0 ⊕ E₊)
    Ω0 = _apply_setops(Ω0, alg)
    X = stateset(ivp)
    Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
    return InitialValueProblem(Sdis, Ω0)
end

function discretize(ivp::IVP{<:CLCS, Interval{N, IA.Interval{N}}}, δ, alg::Forward) where {N}
    A = state_matrix(ivp)
    @assert size(A, 1) == 1
    X0 = initial_state(ivp)

    a = A[1, 1]
    aδ = a * δ
    Φ = exp(aδ)
    A_abs = abs(a)

    # use inverse method
    @assert !iszero(a) "the given matrix should be invertible"

    # a_sqr = a * a
    #P2A_abs = (1/a_sqr) * (Φ - one(N) - aδ)
    #Einit = (P2A_abs * a_sqr) * RA._symmetric_interval_hull(X0).dat

    #P2A_abs = (1/a_sqr) * (Φ - one(N) - aδ)
    Einit = (Φ - one(N) - aδ) * _symmetric_interval_hull(X0).dat

    Ω0 = Interval(hull(X0.dat, Φ * X0.dat + Einit))
    X = stateset(ivp)
    # the system constructor creates a matrix
    Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
    return InitialValueProblem(Sdis, Ω0)
end

# ------------------------------------------------------------
# Forward Approximation: Inhomogeneous case
# ------------------------------------------------------------

# TODO : specialize, add option to compute the concrete linear map
function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ, alg::Forward)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)

    Φ = _exp(A, δ, alg.exp)
    A_abs = _elementwise_abs(A)
    Φcache = sum(A) == abs(sum(A)) ? Φ : nothing
    P2A_abs = _Φ₂(A_abs, δ, alg.exp, alg.inv, Φcache)

    Einit = sih(P2A_abs * sih((A * A) * X0, alg.sih), alg.sih)

    U = next_set(inputset(ivp), 1)
    Eψ0 = sih(P2A_abs * sih(A * U, alg.sih), alg.sih)

    Ud = δ*U ⊕ Eψ0
    In = IdentityMultiple(one(eltype(A)), size(A, 1))

    Ω0 = ConvexHull(X0, Φ * X0 ⊕ Ud ⊕ Einit)
    Ω0 = _apply_setops(Ω0, alg.setops)
    X = stateset(ivp)
    Sdis = ConstrainedLinearControlDiscreteSystem(Φ, In, X, Ud)
    return InitialValueProblem(Sdis, Ω0)
end


text:
#module SimplePlot

#using NMR.Scaling
using .SimpleGraphics
using Printf

export Plot, Axes, Frame, Show, ContourPlot

const FrameHdef = 700
const FrameWdef = 1000
const tickLength = 8
const FontSize = 20
const DefaultFont = "Helvetica"

const xAxesLabelAttr = GraphicsAttributes([
  "font-family"=>DefaultFont,
  "font-size"=>"$FontSize",
  "frame"=>"black",
  "text-anchor"=>"middle"])

const yAxesLabelAttr = GraphicsAttributes([
    "font-family"=>DefaultFont,
    "font-size"=>"$FontSize",
    "frame"=>"black",
    "text-anchor"=>"end"])

const AxesLabelAttr = GraphicsAttributes([
        "font-family"=>DefaultFont,
        "font-size"=>"$(1.2*FontSize)",
        "frame"=>"black",
        "text-anchor"=>"middle"])

const FrameAttr = GraphicsAttributes("fill"=>"none","stroke"=>"black","stroke-width"=>"1.0")
const TickAttr = GraphicsAttributes("stroke"=>"black","stroke-width"=>"0.7")
const PlotAttr = GraphicsAttributes("stroke"=>"#0000b0","stroke-width"=>"5.0","fill"=>"none")

const LineColours = ["#0000b0","#0b0000","#000b00"] ;

function Frame(xscale::NiceScale,yscale::NiceScale;
                Reverse=[false,false],
                xlabel="",
                ylabel="",
                FrameW=FrameWdef,
                FrameH=FrameHdef)
   g=Group()

   # create Frame
   SimpleGraphics.push!(g,Rectangle(-FrameW/2,-FrameH/2,FrameW,FrameH,FrameAttr))

   # create Frame Ticks
   for x in xscale
     lbl=@sprintf("%.4g",x)
     xtick=FrameW*(Reverse[1] ? -1 : 1)*scaled(x,xscale)
    #  println( (x,xtick) )
     SimpleGraphics.push!(g,Polygon(VertexList([xtick -FrameH/2 ; xtick -FrameH/2+tickLength]),TickAttr))
     SimpleGraphics.push!(g,TextElement(VertexList([xtick -FrameH/2-FontSize]),"$lbl",xAxesLabelAttr))
   end

   xp = -FrameH/2
   yp = -FrameW/2

   if xlabel != ""
     SimpleGraphics.push!(g,TextElement(VertexList([0. yp-3*FontSize]),xlabel,AxesLabelAttr))
   end
   if ylabel != ""
     SimpleGraphics.push!(g,TextElement(VertexList([xp-4.0*FontSize 0.0]),ylabel,AxesLabelAttr,-90.0))
   end

   for y in yscale
     lbl=@sprintf("%.4g",y)
     ytick=FrameH*(Reverse[2] ? -1 : 1)*(scaled(y,yscale))
     SimpleGraphics.push!(g,Polygon(VertexList([-FrameW/2 ytick ; -FrameW/2+tickLength ytick]),TickAttr))
     SimpleGraphics.push!(g,TextElement(VertexList([-FrameW/2-FontSize/2 ytick+2]),"$lbl",yAxesLabelAttr))
   end

   return g
end

function Axes(xscale::NiceScale,yscale::NiceScale;
              xAxis=true,
              yAxis=true,
              autoPosition=true,
              Position=[0. 0.],
              Reverse=[false,false],
              xlabel="",
              ylabel="",
              FrameW=FrameWdef,
              FrameH=FrameHdef)
  g=Group()

  # if x scale incorporates 0, put axis there.
  # if scale is purely negative, put axis at the top.
  # if scale is purely positive, put axis at the bottom.

  if autoPosition
    yAxisPos = 0.0
    if xscale.minPoint > 0 yAxisPos = xscale.minPoint end
    if xscale.maxPoint < 0 yAxisPos = xscale.maxPoint end

    xAxisPos = 0.0
    if yscale.minPoint > 0 xAxisPos = yscale.minPoint end
    if yscale.maxPoint < 0 xAxisPos = yscale.maxPoint end
  else
    (yAxisPos,xAxisPos) = Position
  end

  # create Axis Ticks

  if xAxis
    yp = FrameH*(scaled(xAxisPos,yscale))
    SimpleGraphics.push!(g,Polygon(VertexList([-FrameW/2 yp; FrameW/2 yp]),FrameAttr))

    for x=xscale.niceMin:xscale.tickSpacing:xscale.niceMax
      lbl=@sprintf("%.4g",x)
      xtick= FrameW*(Reverse[1] ? -1 : 1)*scaled(x,xscale)
      SimpleGraphics.push!(g,Polygon(VertexList([xtick yp; xtick yp+tickLength]),TickAttr))
      SimpleGraphics.push!(g,TextElement(VertexList([xtick yp-FontSize]),"$lbl",xAxesLabelAttr))
    end

    if xlabel != ""
      SimpleGraphics.push!(g,TextElement(VertexList([0. yp-2.5*FontSize]),xlabel,AxesLabelAttr))
    end

  end

  if yAxis
    xp = FrameW*(scaled(yAxisPos,xscale))
    SimpleGraphics.push!(g,Polygon(VertexList([xp -FrameH/2; xp FrameH/2]),FrameAttr))
    for y=yscale.niceMin:yscale.tickSpacing:yscale.niceMax
      lbl=@sprintf("%.4g",y)
      ytick= FrameH*(Reverse[2] ? -1 : 1)*scaled(y,yscale)
      SimpleGraphics.push!(g,Polygon(VertexList([xp ytick ; xp+tickLength ytick]),TickAttr))
      SimpleGraphics.push!(g,TextElement(VertexList([xp-FontSize/2 ytick+2]),"$lbl",yAxesLabelAttr))
    end

    if ylabel != ""
	    SimpleGraphics.push!(g,TextElement(VertexList([xp-4.0*FontSize 0.0]),ylabel,AxesLabelAttr,-90.0))
    end

  end

  return g
end

function Plot(xrange,pts::Array{T,2};
              g=Group(),
              frame=true,
              axes=false,
              show=true,
              yscale=NiceScale(minimum(pts[:,1]),maximum(pts[:,1])),
              xscale=NiceScale(minimum(xrange),maximum(xrange)),
              Reverse=[false,false],
              style=GraphicsAttributes([]),
              FrameW=FrameWdef,
	      FrameH=FrameHdef) where {T}

    revrs=Reverse
    if frame g=Frame(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end
    if axes  g=Axes(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end

    (xreverse,yreverse)=Reverse

    xs = [ FrameW*(xreverse ? -1 : 1)*scaled(x,xscale) for x in xrange]

    n,m=size(pts)
    for k=1:m
      lineCol=GraphicsAttributes("stroke"=>LineColours[(k-1)%length(LineColours)+1])
      ys = [ FrameH*(yreverse ? -1 : 1)*scaled(y,yscale)  for y in pts[:,k]]
      SimpleGraphics.push!(g,Polygon( VertexList([xs ys]), merge(PlotAttr,lineCol,style)))
    end

    if show
        Show(g)
    end

    return g
end

function Plot(xrange,pts;
              g=Group(),
              frame=true,
              axes=false,
              show=true,
              yscale=NiceScale(minimum(pts),maximum(pts)),
              xscale=NiceScale(minimum(xrange),maximum(xrange)),
              Reverse=[false,false],
              style=GraphicsAttributes([]),
              FrameW=FrameWdef,
              FrameH=FrameHdef)

    revrs=Reverse
    if frame g=Frame(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end
    if axes  g=Axes(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end

    (xreverse,yreverse)=Reverse

    xs = [ FrameW*(xreverse ? -1 : 1)*scaled(x,xscale) for x in xrange]
    ys = [ FrameH*(yreverse ? -1 : 1)*scaled(y,yscale)  for y in pts]

    SimpleGraphics.push!(g,Polygon( VertexList([xs ys]), merge(PlotAttr,style)))

    if show
        Show(g)
    end

    return g
end

#end # SimplePlot

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
# Copyright 2021 Adrian Hill
#
# Licensed under the ImageMagick License (the "License"); you may not use
# this file except in compliance with the License.  You may obtain a copy
# of the License at
#
#   https://imagemagick.org/script/license.php
#
# 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.

## Ordered dithering matrices from ImageMagic
"""
    IM_checks()

ImageMagick's Checkerboard 2x2 dither
"""
IM_checks(; kwargs...) = OrderedDither(CHECKS; kwargs...)
const CHECKS = [
    1 2
    2 1
]//3

## ImageMagic's Halftones - Angled 45 degrees
# Initially added to ImageMagick by Glenn Randers-Pehrson, IM v6.2.8-6,
# modified to be more symmetrical with intensity by Anthony, IM v6.2.9-7

"""
    IM_h4x4a()

ImageMagick's Halftone 4x4 - Angled 45 degrees
"""
IM_h4x4a(; kwargs...) = OrderedDither(H4X4A; kwargs...)
const H4X4A = [
    4 2 7 5
    3 1 8 6
    7 5 4 2
    8 6 3 1
]//9

"""
    IM_h6x6a()

ImageMagick's Halftone 6x6 - Angled 45 degrees
"""
IM_h6x6a(; kwargs...) = OrderedDither(H6X6A; kwargs...)
const H6X6A =
    [
        14 13 10 8 2 3
        16 18 12 7 1 4
        15 17 11 9 6 5
        8 2 3 14 13 10
        7 1 4 16 18 12
        9 6 5 15 17 11
    ]//19

"""
    IM_h8x8a()

ImageMagick's Halftone 8x8 - Angled 45 degrees
"""
IM_h8x8a(; kwargs...) = OrderedDither(H8X8A; kwargs...)
const H8X8A =
    [
        13 7 8 14 17 21 22 18
        6 1 3 9 28 31 29 23
        5 2 4 10 27 32 30 24
        16 12 11 15 20 26 25 19
        17 21 22 18 13 7 8 14
        28 31 29 23 6 1 3 9
        27 32 30 24 5 2 4 10
        20 26 25 19 16 12 11 15
    ]//33

# ImageMagic's Halftones - Orthogonally Aligned, or Un-angled
# Initially added by Anthony Thyssen, IM v6.2.9-5 using techniques from
# "Dithering & Halftoning" by Gernot Haffmann
# http://www.fho-emden.de/~hoffmann/hilb010101.pdf

"""
    IM_h4x4o()

ImageMagick's Halftone 4x4 - Orthogonally Aligned
"""
IM_h4x4o(; kwargs...) = OrderedDither(H4X4O; kwargs...)
const H4X4O = [
    7 13 11 4
    12 16 14 8
    10 15 6 2
    5 9 3 1
]//17

"""
    IM_h6x6o()

ImageMagick's Halftone 6x6 - Orthogonally Aligned
"""
IM_h6x6o(; kwargs...) = OrderedDither(H6X6O; kwargs...)
const H6X6O =
    [
        7 17 27 14 9 4
        21 29 33 31 18 11
        24 32 36 34 25 22
        19 30 35 28 20 10
        8 15 26 16 6 2
        5 13 23 12 3 1
    ]//37

"""
    IM_h8x8o()

ImageMagick's Halftone 8x8 - Orthogonally Aligned
"""
IM_h8x8o(; kwargs...) = OrderedDither(H8X8O; kwargs...)
const H8X8O =
    [
        7 21 33 43 36 19 9 4
        16 27 51 55 49 29 14 11
        31 47 57 61 59 45 35 23
        41 53 60 64 62 52 40 38
        37 44 58 63 56 46 30 22
        15 28 48 54 50 26 17 10
        8 18 34 42 32 20 6 2
        5 13 25 39 24 12 3 1
    ]//65

## ImageMagic's Halftones - Orthogonally Expanding Circle Patterns
# Added by Glenn Randers-Pehrson, 4 Nov 2010, ImageMagick 6.6.5-6
"""
    IM_c5x5()

ImageMagick's Halftone 5x5 - Orthogonally Expanding Circle Patterns
"""
IM_c5x5(; kwargs...) = OrderedDither(C5X5; kwargs...)
const C5X5 = [
    1 21 16 15 4
    5 17 20 19 14
    6 21 25 24 12
    7 18 22 23 11
    2 8 9 10 3
]//26

"""
    IM_c6x6()

ImageMagick's Halftone 6x6 - Orthogonally Expanding Circle Patterns
"""
IM_c6x6(; kwargs...) = OrderedDither(C6X6; kwargs...)
const C6X6 =
    [
        1 5 14 13 12 4
        6 22 28 27 21 11
        15 29 35 34 26 20
        16 30 36 33 25 19
        7 23 31 32 24 10
        2 8 17 18 9 3
    ]//37

"""
    IM_c7x7()

ImageMagick's Halftone 7x7 - Orthogonally Expanding Circle Patterns
"""
IM_c7x7(; kwargs...) = OrderedDither(C7X7; kwargs...)
const C7X7 =
    [
        3 9 18 28 17 8 2
        10 24 33 39 32 23 7
        19 34 44 48 43 31 16
        25 40 45 49 47 38 27
        20 35 41 46 42 29 15
        11 21 36 37 28 22 6
        4 12 13 26 14 5 1
    ]//50


text:
module BaJakob

    using SearchBall
    using GeometryTypes
    using Reactive

    function calc_myPlayers(config::SearchBall.Config, dt::Number,
        set_type::Symbol, moved_object::Union{SearchBall.Object{<:Number}, Nothing}=nothing)

        world_state = config.ws

        try
            # Initialization
            if set_type == :init
                config.printables["dist:opp-off"] = Signal("")
                config.printables["dist:opp-ball"] = Signal("")
                config.printables["dist:off-ball"] = Signal("")
            end

        catch error
        end

        return world_state.myPlayers
    end

    function calc_opponend(config::SearchBall.Config,
        dt::Number, set_type::Symbol,
        moved_object::Union{SearchBall.Object{<:Number}, Nothing}=nothing)

        world_state = config.ws
        if set_type == :init
            # Initialization
        end

        return world_state.opponents
    end

    function run_strategies(config::SearchBall.Config, dt::Number = 0.1,
        set_type::Symbol=:default,
        moved_object::Union{SearchBall.Object{<:Number}, Nothing}=nothing)

        return (
            calc_myPlayers(config, dt, set_type, moved_object),
            calc_opponend(config, dt, set_type, moved_object)
        )
    end

end


text:
using Test
using Random
using StatsBase
using Stheno
using StringDistances
using Distances
using Stheno: ew, pw
using LinearAlgebra
using BOMoD
using BOMoD: filterconstraint, Combination,  ComposedConstructConstraints,
    Mod,GroupMod, _word2vec, degree, laplacian,norm_laplacian, _len1,_creatGP

include("./BO_test/test_functions.jl")

@testset "BOMod" begin
    include("./design_test/Mod_test.jl")
    include("./design_test/Filter_test.jl")
    include("./design_test/Space_test.jl")
    include("./design_test/Constraint_test.jl")
    include("./design_test/Combination_test.jl")
    include("./design_test/Constructs_test.jl")
    include("./BO_test/Kernel_test.jl")
    include("./BO_test/graph_test.jl")
    include("./BO_test/GP_test.jl")
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
#!/usr/bin/env julia

using WriteVTK

using Printf: @sprintf

const FloatType = Float32
const vtk_filename_noext = "collection"

function update_point_data!(p, q, vec, time)
    E = exp(-time)
    for I in CartesianIndices(p)
        i, j, k = I[1], I[2], I[3]
        p[I] = E * (i * i + k)
        q[I] = E * k * sqrt(j)
        vec[1, i, j, k] = E * i
        vec[2, i, j, k] = E * j
        vec[3, i, j, k] = E * k
    end
    p, q, vec
end

function update_cell_data!(cdata, time)
    Nj = size(cdata, 2) + 1
    α = 3pi / (Nj - 2)
    E = exp(-time)
    for I in CartesianIndices(cdata)
        i, j, k = I[1], I[2], I[3]
        cdata[I] = E * (2i + 3k * sin(α * (j - 1)))
    end
    cdata
end

function main()
    # Define grid.
    Ni, Nj, Nk, Nt = 20, 30, 40, 4

    x = zeros(FloatType, Ni)
    y = zeros(FloatType, Nj)
    z = zeros(FloatType, Nk)

    [x[i] = i*i/Ni/Ni for i = 1:Ni]
    [y[j] = sqrt(j/Nj) for j = 1:Nj]
    [z[k] = k/Nk for k = 1:Nk]

    # Arrays for scalar and vector fields assigned to grid points.
    p = zeros(FloatType, Ni, Nj, Nk)
    q = zeros(FloatType, Ni, Nj, Nk)
    vec = zeros(FloatType, 3, Ni, Nj, Nk)

    # Scalar data assigned to grid cells.
    # Note that in structured grids, the cells are the hexahedra formed between
    # grid points.
    cdata = zeros(FloatType, Ni - 1, Nj - 1, Nk - 1)

    # Test extents (this is optional!!)
    ext = [0, Ni-1, 0, Nj-1, 0, Nk-1] .+ 42

    # Initialise pvd container file
    @time outfiles = paraview_collection(vtk_filename_noext) do pvd
        # Create files for each time-step and add them to the collection
        for it = 0:Nt-1
            vtk = vtk_grid(@sprintf("%s_%02i", vtk_filename_noext, it), x, y, z;
                           extent=ext)
            # Add data for current time-step
            update_point_data!(p, q, vec, it + 1)
            update_cell_data!(cdata, it + 1)
            vtk_point_data(vtk, p, "p_values")
            vtk_point_data(vtk, q, "q_values")
            vtk_point_data(vtk, vec, "myVector")
            vtk_cell_data(vtk, cdata, "myCellData")
            vtk_save(vtk)
            collection_add_timestep(pvd, vtk, Float64(it+1))
        end
    end

    # Create a copy of above pvd for reloading
    cp(vtk_filename_noext * ".pvd",
       vtk_filename_noext * "_reload.pvd",
       force=true)
    pvd_reload = paraview_collection(vtk_filename_noext * "_reload", append=true)

    # add a vtk file
    vtk_reload = vtk_grid("collection_reload", [1, 2, 3], [1, 2, 3])
    collection_add_timestep(pvd_reload, vtk_reload, 5.0)
    pvd_reload_files = vtk_save(pvd_reload)
    append!(outfiles, pvd_reload_files)

    println("Saved:  ", join(outfiles, "  "))

    return outfiles::Vector{String}
end

main()


text:
# Basic Syntax

function f(x, y)
  x + y
end

println(f(2, 3))

g = f
println(g(2, 3))

# Pass by Sharing Function Behavior
arr = Int64[1, 2, 3]
i = 2

function h(x, i)
  x[i] -= 1
  i -= 1
end

println(h(arr, i));
println(arr);
println(i);

# Returning from a function

function divider(x, y)
  if (y != 0)
    return x / y
  end
  return "Divide by zero error!"
end

println(divider(1, 2))
println(divider(1, 0))

# Return a specific type

function divider(x, y) :: Int32
  if (y != 0)
    return x / y
  end
  return -1
end

# return nothing

function divider(x, y)
  if (y != 0)
    return x / y
  end
  return nothing
end

println(divider(1, 0))
println(typeof(divider(1, 0)));

# multiple Return Values with Tuples

function fun(x, y)
  x += 1
  y += 1
  x, y
end

println(fun(1, 2))
println(typeof(fun(1, 2)))

a, b = fun(1, 2)

# Operators are actually functions!

1 + 2 + 3

+(1, 2, 3)
-(4, 5)

# Anonymous functions
function (x)
  x^2 + 3x - 2
end

x -> x^2 + 3x - 2

println(map(x -> x^2 + 3x - 2, [1, 2, 3]))

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
module Mesh

importall ImmutableArrays
import GL

type Texture
	diffuse::GL.Texture
	normal::GL.Texture
	width::Uint32
	height::Uint32
	flags::Uint32
end

type Light
	origin::Vector3{Float32}
	color::Vector3{Float32}
	power::Float32
	Light(origin, color, power) = new(convert(Vector3{Float32}, origin), convert(Vector3{Float32}, color), convert(Float32, power))
end

type Face
	indices::Vector{Uint16}
	ibo::GL.Buffer
	normal::Vector3{Float32}
	texture::Texture
	u_axis::Vector4{Float32}
	v_axis::Vector4{Float32}
	lights::Vector{Light}
end

export draw
function draw(face::Face)
	GL.DrawElements(GL.TRIANGLES, face.indices)
end

typealias Vertex Vector3{Float32}

function islit(face::Face, vertices::Vector{Vertex}, light::Light)

	# skip faces that have no triangles
	if length(face.indices) < 3
		return false
	end

	# find d component of face's plane
	v = vertices[face.indices[1]+1]
	d = -dot(face.normal, v)

	# is light radius within range of face's plane?
	dist = dot(face.normal, light.origin) + d
	if dist < 0 || dist > light.power
		return false
	end

	# project light onto face's plane
	projected_origin = light.origin - face.normal * dist
	projected_power = light.power * sin(acos(dist / light.power))

	# check if projected light radius intersects any triangles in face
	for tri = 0:floor(length(face.indices)/3)-1
		inside_triangle = true
		for edge = 0:2
			i1 = tri*3+edge+1
			i2 = tri*3+(edge+1)%3+1
			v1 = vertices[face.indices[i1]+1]
			v2 = vertices[face.indices[i2]+1]
			edge_normal = cross(face.normal, v1-v2)
			edge_normal /= norm(edge_normal)
			d = -dot(edge_normal, v1)
			dist = dot(edge_normal, projected_origin) + d
			if abs(dist) < projected_power
				return true
			end
			if dist > 0
				inside_triangle = false
			end
		end
		if inside_triangle
			return true
		end
	end

	return false
end

end


text:
include("io_lammps.jl")
include("io_xyz.jl")
include("io_json.jl")


text:
module DocumenterCitations
global indexer = 1
global tracker = Dict()
using Documenter
using Documenter.Anchors
using Documenter.Builder
using Documenter.Documents
using Documenter.Selectors
using Documenter.Utilities
using Documenter.Expanders

using Markdown
using Bibliography
using Bibliography: xnames, xyear, xlink, xtitle, xin
using Unicode

export CitationBibliography
struct CitationBibliography <: Documenter.Plugin
    bib::Dict
end

function CitationBibliography(filename::AbstractString="")
    # global indexer = 1
    filename == "" && return CitationBibliography(Dict())
    bf = import_bibtex(filename)
    return CitationBibliography(bf)
end

struct Example end

include("citations.jl")
include("bibliography.jl")

end

max_stars_repo_name: timknab/DocumenterCitations

max_issues_repo_name: timknab/DocumenterCitations

max_forks_repo_name: timknab/DocumenterCitations


text:
module SparseBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

import Compat: UTF8String, view

const SUITE = BenchmarkGroup(["array"])

#########
# index #
#########

# Note that some of the "logical" tests are commented
# out because they require resolution of JuliaLang/julia#14717.

g = addgroup!(SUITE, "index")

# vector #
#--------#

sizes = (10^3, 10^4, 10^5)

if VERSION < v"0.5.0-dev+763"
    spvecs = map(s -> samesprand(s, 1, inv(sqrt(s))), sizes)
    splogvecs = map(s -> samesprandbool(s, 1, 1e-5), sizes)
else
    spvecs = map(s -> samesprand(s, inv(sqrt(s))), sizes)
    splogvecs = map(s -> samesprandbool(s, 1e-5), sizes)
end

for (s, v, l) in zip(sizes, spvecs, splogvecs)
    g["spvec", "array",   s] = @benchmarkable getindex($v, $(samerand(1:s, s)))
    g["spvec", "integer", s] = @benchmarkable getindex($v, $(samerand(1:s)))
    g["spvec", "range",   s] = @benchmarkable getindex($v, $(1:s))
    g["spvec", "logical", s] = @benchmarkable getindex($v, $(samerand(Bool, s)))
    # g["spvec", "splogical", s, nnz(v), nnz(l)] = @benchmarkable getindex($v, $l)
end

# matrix #
#--------#

sizes = (10, 10^2, 10^3)
inds = map(s -> samerand(1:s), sizes)
matrices = map(s -> samesprand(s, s, inv(sqrt(s))), sizes)
vectors = map(s -> samerand(1:s, s), sizes)
logvecs = map(s -> samerand(Bool, s), sizes)
splogmats = map(s -> samesprandbool(s, s, 1e-5), sizes)

if VERSION < v"0.5.0-dev+763"
    splogvecs = map(s -> samesprandbool(s, 1, 1e-5), sizes)
else
    splogvecs = map(s -> samesprandbool(s, 1, 1e-5), sizes)
end

for (s, m, v, l, sl, c) in zip(sizes, matrices, vectors, logvecs, splogvecs, inds)
    g["spmat", "col", "array", s] = @benchmarkable getindex($m, $v, $c)
    g["spmat", "col", "range", s] = @benchmarkable getindex($m, $(1:s), $c)
    if isdefined(Base, :OneTo)
        g["spmat", "col", "OneTo", s] = @benchmarkable getindex($m, $(Base.OneTo(s)), $c)
    end
    g["spmat", "col", "logical", s] = @benchmarkable getindex($m, $l, $c)
    # g["spmat", "col", "splogical", s] = @benchmarkable getindex($m, $sl, $c)
end

for (s, m, v, l, sl, r) in zip(sizes, matrices, vectors, logvecs, splogvecs, inds)
    g["spmat", "row", "array", s] = @benchmarkable getindex($m, $r, $v)
    g["spmat", "row", "range", s] = @benchmarkable getindex($m, $r, $(1:s))
    if isdefined(Base, :OneTo)
        g["spmat", "row", "OneTo", s] = @benchmarkable getindex($m, $r, $(Base.OneTo(s)))
    end
    g["spmat", "row", "logical", s] = @benchmarkable getindex($m, $r, $l)
    # g["spmat", "row", "splogical", s] = @benchmarkable getindex($m, $r, $sl)
end

for (s, m, v, l, sl, i) in zip(sizes, matrices, vectors, logvecs, splogmats, inds)
    g["spmat", "array", s] = @benchmarkable getindex($m, $v, $v)
    g["spmat", "integer", s] = @benchmarkable getindex($m, $i, $i)
    g["spmat", "range", s] = @benchmarkable getindex($m, $(1:s), $(1:s))
    if isdefined(Base, :OneTo)
        g["spmat", "OneTo", s] = @benchmarkable getindex($m, $(Base.OneTo(s)), $(Base.OneTo(s)))
    end
    g["spmat", "logical", s] = @benchmarkable getindex($m, $l, $l)
    g["spmat", "splogical", s] = @benchmarkable getindex($m, $sl)
end

for b in values(g)
    b.params.time_tolerance = 0.3
end

######################
# transpose (#14631) #
######################

small_sqr = samesprand(600, 600, 0.01)
small_rct = samesprand(600, 400, 0.01)
large_sqr = samesprand(20000, 20000, 0.01)
large_rct = samesprand(20000, 10000, 0.01)

if VERSION >= v"0.7.0-DEV.1415"
    g = addgroup!(SUITE, "transpose", ["adjoint"])
else
    g = addgroup!(SUITE, "transpose", ["ctranspose"])
end

for m in (small_sqr, small_rct, large_sqr, large_rct)
    cm = m + m*im
    s = size(m)
    g["transpose", s] = @benchmarkable transpose($m)
    g["transpose!", s] = @benchmarkable transpose!($(m.'), $m)
    if VERSION >= v"0.7.0-DEV.1415"
        g["adjoint", s] = @benchmarkable adjoint($cm)
        g["adjoint!", s] = @benchmarkable adjoint!($(cm.'), $cm)
    else
        g["ctranspose", s] = @benchmarkable ctranspose($cm)
        g["ctranspose!", s] = @benchmarkable ctranspose!($(cm.'), $cm)
    end
end

for b in values(g)
    b.params.time_tolerance = 0.3
end

##############
# arithmetic #
##############

g = addgroup!(SUITE, "arithmetic")

# unary minus, julialang repo issue #19503 / fix #19530
g["unary minus", size(small_sqr)] = @benchmarkable -$small_sqr
g["unary minus", size(large_sqr)] = @benchmarkable -$large_sqr

for b in values(g)
    b.params.time_tolerance = 0.3
end

################
# constructors #
################
g = addgroup!(SUITE, "constructors")

const UPLO = VERSION >= v"0.7.0-DEV.884" ? :U : true
for s in sizes
    nz = floor(Int, 1e-4*s*s)
    I = samerand(1:s, nz)
    J = samerand(1:s, nz)
    V = randvec(nz)
    g["IV", s] = @benchmarkable sparsevec($I, $V)
    g["IJV", s] = @benchmarkable sparse($I, $J, $V)
    g["Diagonal", s] = @benchmarkable sparse($(Diagonal(randvec(s))))
    g["Bidiagonal", s] = @benchmarkable sparse($(Bidiagonal(randvec(s), randvec(s-1), UPLO)))
    g["Tridiagonal", s] = @benchmarkable sparse($(Tridiagonal(randvec(s-1), randvec(s), randvec(s-1))))
    g["SymTridiagonal", s] = @benchmarkable sparse($(SymTridiagonal(randvec(s), randvec(s-1))))
end

end # module


text:
# Note that indexing starts at 1 in Julia
for i in 1:5
    print(i, '\n')
end


text:
@testset "Coordinates" begin

    @testset "Create" begin
        @testset "Brain Vision" begin
            bv = BrainVision(0, 0, 0)
            @testset "Show" begin
                 show(bv)
            end
        end
        @testset "Talairach" begin
            tal = Talairach(68.3, -26.9, 8.3)
            @testset "Show" begin
                 show(tal)
            end
        end
        @testset "SPM" begin
            spm = SPM(68.3, -26.9, 8.3)
            @testset "Show" begin
                 show(spm)
            end
        end
        @testset "Unknown" begin
            uk = UnknownCoordinate(1.2, 2, 3.2)
            @testset "Show" begin
                 show(uk)
            end
        end
    end

    @testset "Convert" begin
        @testset "MNI -> Talairach" begin

            # Values from table IV in Lancaster et al 2007
            # TODO Find better points to translate

            mni = SPM(73.7, -26.0, 7.0)
            tal = Talairach(68.3, -26.9, 8.3)
            @test isapprox(euclidean(convert(Talairach, mni), tal), 0; atol=2)

            # As an electrode

            # Test as an electrode
            e = Electrode("test", SPM(73.7, -26.0, 7.0), Dict())
            e = conv_spm_mni2tal(e)

            @test isapprox(e.coordinate.x, tal.x; atol = 1.5)
            @test isapprox(e.coordinate.y, tal.y; atol = 1.5)
            @test isapprox(e.coordinate.z, tal.z; atol = 1.5)
            @test isa(e, EEG.Sensor) == true
            @test isa(e, EEG.Electrode) == true
            @test isa(e.coordinate, EEG.Talairach) == true

        end

        @testset "BrainVision -> Talairach" begin

            # TODO this just tests it runs, need to check values

            bv = BrainVision(0, 0, 0)
            tal = Talairach(128, 128, 128)
            @test isapprox(euclidean(convert(Talairach, bv), tal), 0; atol=1.5)

        end
    end

    @testset "Distances" begin

        @test euclidean(Talairach(0, 0, 0), Talairach(1, 1, 1)) == sqrt.(3)

        v = [0, 0, 0]
        @test euclidean(Talairach(1, 1, 1), v) == sqrt.(3)
        @test euclidean(v, Talairach(1, 1, 1)) == sqrt.(3)

    end
end


text:
function check_tolerance(::Type{Sphere1D}, arg0::jdouble)
    return jcall(Sphere1D, "checkTolerance", void, (jdouble,), arg0)
end

function get_dimension(obj::Sphere1D)
    return jcall(obj, "getDimension", jint, ())
end

function get_instance(::Type{Sphere1D})
    return jcall(Sphere1D, "getInstance", Sphere1D, ())
end

function get_sub_space(obj::Sphere1D)
    return jcall(obj, "getSubSpace", Space, ())
end