text:
"""Check whether elements of a vector are increasing."""
function isincreasing(x::Vector)
    is_increasing = true
    for i in 2:length(x)
        is_increasing &= x[i] > x[i-1] ? true : false
    end
    return is_increasing
end

"""Set all negative entries in an array to zero."""
function clip_negatives!(A::AbstractArray{T}) where T
    @inbounds for i in eachindex(A)
        A[i] = max(A[i],zero(T))
    end
end

"""Returns an integer `m >= n` of argument `n` with only prime factors 2 and 3
to obtain an efficiently fourier-transformable number of longitudes.
    m = 2^i * 3^j >= n, with i>=0 but j = 0,1
"""
function roundup_fft(n::Int)
    lz = leading_zeros(n)                   # determine scale of n

    # create a mask that's 1 for all but the two most significant figures of n
    # for finding the next largest integer of n with factors 2 and 3
    mask = (1 << (8*sizeof(n)-lz-2))-1

    # round up by adding mask as offset and then masking all insignificant bits
    n_roundedup = (n + mask) & ~mask
    return n_roundedup
end


text:
module TropicalNumbers

export Tropical, TropicalF64, TropicalF32, TropicalF16, CountingTropicalF16, CountingTropicalF32, CountingTropicalF64, content
export CountingTropical
export TropicalTypes

include("tropical.jl")
include("counting_tropical.jl")

const TropicalTypes{T} = Union{CountingTropical{T}, Tropical{T}}

# alias
for NBIT in [16, 32, 64]
    @eval const $(Symbol(:Tropical, :F, NBIT)) = Tropical{$(Symbol(:Float, NBIT))}
    @eval const $(Symbol(:CountingTropical, :F, NBIT)) = CountingTropical{$(Symbol(:Float, NBIT)),$(Symbol(:Float, NBIT))}
end

# alias
for T in [:Tropical, :CountingTropical]
    # defining constants like `TropicalF64`.
    for OP in [:>, :<, :(==), :>=, :<=, :isless]
        @eval Base.$OP(a::$T, b::$T) = $OP(a.n, b.n)
    end
    @eval begin
        content(x::$T) = x.n
        content(x::Type{$T{X}}) where X = X
        Base.isapprox(x::AbstractArray{<:$T}, y::AbstractArray{<:$T}; kwargs...) = all(isapprox.(x, y; kwargs...))
        Base.show(io::IO, ::MIME"text/plain", t::$T) = Base.show(io, t)
        Base.isnan(x::$T) = isnan(content(x))

        # this is for CUDA matmul
        Base.:(*)(a::$T, b::Bool) = b ? a : zero(a)
        Base.:(*)(b::Bool, a::$T) = b ? a : zero(a)
        Base.:(/)(a::$T, b::Bool) = b ? a : a / zero(a)
        Base.:(/)(b::Bool, a::$T) = b ? inv(a) : zero(a)
        Base.div(a::$T, b::Bool) = b ? a : a / zero(a)
        Base.div(b::Bool, a::$T) = b ? inv(a) : zero(a)
    end
end

end # module

ext: jl


text:
function prepare_variables(objective, approach, x0, ∇fz, B)
    z = x0
    x = copy(z)

    if isa(B, Nothing)  # didn't provide a B
        if isa(first(approach), GradientDescent)
            # We don't need to maintain a dense matrix for Gradient Descent
            B = I
        else
            # Construct a matrix on the correct form and of the correct type
            # with the content of I_{n,n}
            B = I + abs.(0*x*x')
        end
    end
    # first evaluation
    if isa(first(approach), Newton)
        fz, ∇fz, B = objective(x, ∇fz, B)
    else
        fz, ∇fz = objective(x, ∇fz)
    end
    fx = copy(fz)
    ∇fx = copy(∇fz)
    return x, fx, ∇fx, z, fz, ∇fz, B
end

function converged(z, ∇fz, g_tol)
    g_converged = norm(∇fz) < g_tol
    return g_converged || any(isnan.(z))
end

include("inplace_loop.jl")
include("outofplace_loop.jl")


text:
using DrWatson
@quickactivate "dummy_project"
using DrWatsonSim
using BSON
using Dates

function long_running_computation(p,output_path)
    sleep(p[:duration])
    result = p[:a]^p[:b]
    BSON.bson(output_path, Dict(:result=>result))
    return nothing
end

duration = [0.1, 1]
a = [1,3]
b = 1
parameter = @dict duration a b

if in_simulation_mode()
    # Using simid() here is actually better, because then the
    # no lookup it the index file is need which means, that the
    # database must not be loked.
    m = Metadata(simdir())
    println("$(simid()): Loaded metadata file")
    m["type"] = "Simple Computation"
    m["started at"] = Dates.now()
    println("$(simid()): Creating new files")
    m_new = Metadata(simdir("newfile"))
    m_new["extra"] = "This should be blocked"
    println("$(simid()): Done creating new file")
end

@runsync x->long_running_computation(x, simdir("output.bson")) dict_list(parameter) datadir("sims")

if in_simulation_mode()
    @info "Simulation done" simid()
    isfile(simdir("output_first_run.bson")) || mv(simdir("output.bson"),simdir("output_first_run.bson"))
end

if !in_simulation_mode()
    @info "All simulations done"
    for folder in readdir(datadir("sims"))
        path = datadir("sims", folder)
        @rerunsync x->long_running_computation(x, simdir("output.bson")) path
    end
end

max_stars_repo_name: sebastianpech/DrWatsonSim

max_issues_repo_name: sebastianpech/DrWatsonSim

max_forks_repo_name: sebastianpech/DrWatsonSim


text:
module Lattices

import Base: size, ndims, length, show, nameof, ==, axes, strides, checkbounds
using Base.Iterators
using LinearAlgebra

export AbstractLattice, WeightedLattice, CoordinateLattice

export neighbors, bonds, sites, translation_vectors, to_site_id, to_coordinate, nsites

export HyperCubic, Chain, Square, Cubic
export Triangular, Honeycomb, Kagome

export AbstractBoundary, Periodic, Open, Helical

export Coordinate

include("coordinate.jl")
include("boundaries.jl")
include("types.jl")
include("neighbors.jl")
include("sites.jl")
include("bonds.jl")
include("hypercubic.jl")
include("triangular.jl")

#######################
## Boundary Conditions
#######################

function _apply_boundary_conditions(::Val{true}, bcs::NTuple{M, Periodic}, dims::NTuple{M, Int}, site::Coordinate{N, Int})::Coordinate{N, Int} where {M, N}
    return Coordinate(apply_boundary.(bcs, dims, site.coordinates[1:end-1])..., site.coordinates[end])
end

function _apply_boundary_conditions(::Val{true}, bcs::NTuple{M, B}, dims::NTuple{M, Int}, site::Coordinate{N, Int})::Union{Coordinate{N, Int}, Nothing} where {M, N, B <: AbstractBoundary}
    coords = site.coordinates
    clipped = _apply_boundary_conditions(Val{false}(), bcs, dims, Coordinate(coords[1:end-1]))
    if clipped === nothing
        return nothing
    else
        return Coordinate(clipped.coordinates..., coords[end])
    end
end

function _apply_boundary_conditions(::Val{false}, bcs::NTuple{N, Periodic}, dims::NTuple{N, Int}, site::Coordinate{N, Int})::Coordinate{N, Int} where N
    return Coordinate(apply_boundary.(bcs, dims, site.coordinates)...)
end

function _apply_boundary_conditions(::Val{false}, bcs::NTuple{N, B}, dims::NTuple{N, Int}, site::Coordinate{N, Int})::Union{Coordinate{N, Int}, Nothing} where {N, B <: AbstractBoundary}
    coords = site.coordinates
    clipped = zeros(Int, N)
    for i in 1:N
        c = apply_boundary(bcs[i], dims[i], coords[i])
        if isnothing(c)
            return nothing
        else
            clipped[i] = c
        end
    end
    return Coordinate(clipped...)
end

for L in [:Triangular, :Honeycomb, :Kagome, :HyperCubic]

    @eval function apply_boundary_conditions(lattice::$L, site::Coordinate)
        return _apply_boundary_conditions(hasunitcell(lattice), lattice.bcs, lattice.dims, site)
    end

end

end


text:
using ImplicitEquations
using Base.Test

f(x,y) = y-x
g(x,y) = y+x

ImplicitEquations.GRAPH(f == 0, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f <= 0, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f <  g, -5, 5, -5, 5, 2^4, 2^4)

ImplicitEquations.GRAPH(f ≶  g, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≷  g, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(neq(f, g), -5, 5, -5, 5, 2^4, 2^4)

ImplicitEquations.GRAPH(eq(f, g), -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ⩵ g, -5, 5, -5, 5, 2^4, 2^4)


text:
using GeometricAlgebra
using StaticArrays

const ga = GeometricAlgebra

@basis "++-"

blade = Blade(1, 4)

include("blades.jl")

spacetime = Signature(3, 1)
ℝ⁴ = Signature(4, 0)
sig_111 = Signature(1, 1, 1)

include("signatures.jl")

include("basis.jl")

mv_1 = convert(Multivector, 1v1 + 1v2)

mv_2 = 1v1 + 1v12 + 1v123

A = 1.3v1 + 2.7v12
B = 0.7v3 + 1.1v123 + 3.05v

include("multivectors.jl")


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
matplotlib.use("agg") # or try "macosx"
using SparseArrays

trialnum = 1
include("chip_unstructured_geometry.jl")
include("plot_inverse_one_iter.jl")

k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

function k_exact(x, y)
    k_mold + 1000 * k_chip_ref * (x-0.49)^2 / (1 + x^2)
end

function k_nn(x, y)
    out = abs(fc(x, [20,20,20,1])) .+ k_mold
    squeeze(out)
end

nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

NT = 15    # number of iterations for Newton's method, 8 is good for m=400

xy = mesh.nodes
xy2 = zeros(mesh.nedge, 2)
for i = 1:mesh.nedge
    xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
end
xy = [xy;xy2]

x, y = xy[chip_fem_idx, 1], xy[chip_fem_idx, 2]
k_chip = k_nn(x, y)

k_fem = k_air * constant(ones(ndof))
k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
k_fem = scatter_update(k_fem, chip_fem_idx, k_chip)
kgauss = dof_to_gauss_points(k_fem, mesh)  # try change to include k on edges

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
# heat_source_fem[chip_fem_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)  # try change to include k on edges

B = constant(compute_interaction_matrix(mesh)) # function not exist
# TODO: add unstructured mesh version of compute_interaction_matrix

# compute F
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, mesh))
# xy = fem_nodes(m, n, h)
# x, y = xy[:, 1], xy[:, 2]
# k = @. k_nn(x, y); k=stack(k)
# kgauss = dof_to_gauss_points(k, mesh)
# LaplaceK = compute_fem_laplace_matrix1(kgauss, mesh)

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
     bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

function compute_residual(S)
    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]
    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    interaction = compute_interaction_term(p, mesh) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, mesh)
    f2 = compute_fem_source_term1(vgauss.*uy, mesh)
    f3 = -interaction[1:ndof]
    f4 = Laplace*u
    # f5 = -F1
    F = f1 + f2 + f3 + f4 #+ f5

    g1 = compute_fem_source_term1(ugauss.*vx, mesh)
    g2 = compute_fem_source_term1(vgauss.*vy, mesh)
    g3 = -interaction[ndof+1:end]
    g4 = Laplace*v
    # g5 = -F2
    T_gauss = dof_to_gauss_points(T, mesh)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, mesh)

    G = g1 + g2 + g3 + g4 + buoyance_term #+ g5

    H0 = -B * [u;v] # + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, mesh) * T - heat_source
    R = [F;G;H0;T0]
    # return R  # original
    # return[F;G;H0;constant(zeros(ndof))]  # fix T=0
    # return[F;G;constant(zeros(nelem));constant(zeros(ndof))]  # fix p=0 and T=0
end

function compute_jacobian(S)
    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]

    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    M1 = constant(compute_fem_mass_matrix1(ux, mesh))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, mesh))

    N1 = constant(compute_fem_mass_matrix1(vy, mesh))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, mesh))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, mesh))

    gradT = eval_grad_on_gauss_pts1(T, mesh)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, mesh))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, mesh))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = constant(compute_fem_mass_matrix1(-buoyance_coef * constant(ones(ndof)), mesh))
    T_mat = [SparseTensor(spzeros(ndof, ndof)); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(nelem, ndof))]

    # N = 2 * ndof + nelem
    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros(ndof, nelem)) M]]

    # return [J[1:2*ndof+nelem, 1:2*ndof+nelem]  SparseTensor(spzeros(2*ndof+nelem, ndof))
    # SparseTensor(spzeros(ndof,2*ndof+nelem)) spdiag(ones(ndof))] # fix T=0

    # return [J[1:2*ndof, 1:2*ndof]  SparseTensor(spzeros(2*ndof, nelem+ndof))
    # SparseTensor(spzeros(nelem+ndof,2*ndof)) spdiag(ones(nelem+ndof))]  # fix p=0 and T=0
end

function solve_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, mesh)
    residual = scatter_update(residual, bd, zeros(length(bd)))    # residual[bd] .= 0.0 in Tensorflow syntax

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(nelem+3*ndof))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, nelem+3*ndof))

S_computed = S[end, :]
S_data = matread("fn$trialnum/xz_chip_unstructured_data.mat")["V"]

# sample_size = 20
# idx = rand(1:ndof, sample_size)
# idx = [idx; ndof .+ idx; 2*ndof+nelem .+ idx] # observe velocity and temperature
# observed_data = S_data[idx]

# noise = false
# noise_level = 0.05
# if noise
#     noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
#     observed_data = observed_data .* noise_ratio
# end

loss = mean((S_computed .- S_data)^2)
loss = loss * 1e10
# ---------------------------------------------------
# create a session and run
max_iter = 1
sess = Session(); init(sess)

for k = 1:10000
    loss_ = BFGS!(sess, loss, max_iter)
    printstyled("[#iter $k] loss=$loss_\n", color=:green)
    matwrite("fn$trialnum/loss$k.mat", Dict("L"=>loss_))
    close("all"); semilogy(loss_); title("loss vs. iteration")
    savefig("fn$trialnum/nn_loss$k.png")
    plot_velo_pres_temp_cond(k)
    ADCME.save(sess, "fn$trialnum/nn_sess$k.mat")
end


text:
# NSGAII implementation
# source: https://doi.org/10.1007/3-540-45356-3_83

function NSGAII_update(P::Vector{Sol}, config::Config)

    debug && DEBUG_correct_evaluations(P, crowding = true)

    #------ Offsprings -------#

    N = length(P)

    # generate offsprings
    offsprings = Vector{Sol}(undef, N)
    idx_offspring = 0
    while idx_offspring < N

        # choose parents
        parent_1 = SELECTION_binary_tournament(P)
        parent_2 = SELECTION_binary_tournament(P)
        while parent_2 == parent_1
            parent_2 = SELECTION_binary_tournament(P)
        end

        # recombine parents
        offspring = config.crossover(parent_1, parent_2)

        # mutate offspring
        if rand() < config.mutation_rate
            MUTATION_flip!(offspring)
        end

        # repair offspring
        config.reparation!(offspring)

        # check that the offspring is not already in the population
        is_new = true
        for sol in P
            if sol.z == offspring.z
                is_new = false
                break
            end
        end

        if is_new
            # store offspring
            idx_offspring += 1
            offsprings[idx_offspring] = offspring
        end
    end

    #------ Create merged population R, sort it and evaluate it ------#

    R = UTIL_concatenate_solutions_vectors(P,offsprings)

    F = INDICATOR_fast_non_dominated_sort!(R)

    min_size_new_pop = min(2*N, Int(floor(length(F[1])*1.5))) # compute new adaptative size

    N = max(N, min_size_new_pop)

    new_P = Vector{Sol}(undef, N)
    size_new_P = 0
    idx_front = 1

    while size_new_P + length(F[idx_front]) <= N # we can the whole front

        # evaluate front (for next iteration)
        INDICATOR_crowding_distance!(F[idx_front])

        # add the whole front
        for sol_front in F[idx_front]
            size_new_P += 1
            new_P[size_new_P] = sol_front
        end

        idx_front += 1
    end

    if size_new_P < N # we still have space for a few solutions of the next front
        # evaluate front
        INDICATOR_crowding_distance!(F[idx_front])

        # sort front by decreasing crowding distance
        perm = sortperm(F[idx_front], by = sol -> sol.crowding, rev = true)

        idx_sol = 0
        while size_new_P < N
            size_new_P += 1
            idx_sol += 1
            new_P[size_new_P] = F[idx_front][perm[idx_sol]]
        end

    end

    debug && DEBUG_correct_evaluations(new_P, crowding = true)

    return new_P

end

function NSGAII_solve(prob::_MOMKP, config::Config; seeding_sols = Vector{Sol}(undef, 0), max_t = 10.0, verbose = true)

    st = time()

    # create population
    P = INIT_random_feasible_population(config.size_pop_init, prob, seeding_sols = seeding_sols)

    # evaluate population
    F = INDICATOR_fast_non_dominated_sort!(P)
    for front in F
        INDICATOR_crowding_distance!(front)
    end

    # iterate
    it = 0
    while !UTIL_termination_criteria_met(time()-st, max_t)
        it += 1
        P = NSGAII_update(P, config)
    end

    # evaluate final front
    F = INDICATOR_fast_non_dominated_sort!(P)
    INDICATOR_crowding_distance!(F[1])

    last_front = UTIL_remove_doublons_from_front(F[1])
    sort!(last_front, by = sol -> sol.z[1])

    verbose && println("total time spent : ",time() - st)
    verbose && println("total nb of iter : ",it)
    verbose && println("total nb of sols : ",length(last_front),"\n")

    return last_front

end


text:
push!(LOAD_PATH,"../src/")
using Documenter, LinearRegression, StatsModels, DataFrames

makedocs(sitename="LinearRegression.jl")

deploydocs(
    repo = "github.com/ericqu/LinearRegression.jl.git",
    push_preview = false,
    devbranch = "main",
    devurl = "dev",
)


text:
include("./ols.jl")

data = readcsv("data/gasoline.csv")

# Get the columns that correspond with our X values
# and create a new matrix to hold them
# syntax: data[row, column]
# to get a whole column do data[:, column]
# to get a whole row do data[row, :]
x = [data[:,2] data[:,3] data[:,4]]
y = data[:,6]

# Create a new OLS object
reg = ols.tols(y, x, "Octane Rating", ["Error", "Component 1", "Component 2", "Component 3"])

max_stars_repo_path: Module 3/Chapter05/mainols.jl

max_issues_repo_path: Module 3/Chapter05/mainols.jl

max_forks_repo_path: Module 3/Chapter05/mainols.jl


text:
# TODO: make more efficient
"""
    GaussNewton solver
"""
struct GNSolver <: LinearSolver end

function gn_solver(A)
    GNSolver()
end

function linear_solve!(solver::GNSolver, x::Vector{T}, A::SparseMatrixCSC{T,Int}, b::Vector{T}) where T
    ldiv!(x, lu!(Array(A' * A)), Array(A' * b)) # solve
end

function linear_solve!(solver::GNSolver, x::Matrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T
    ldiv!(x, lu!(A' * A), A' * B) # solve
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
##Tests for the Small-Data KP
include("KPNormal.jl")
using Distributions, Random, KP

#Default simulation of normal random variates
function simData!(o, xs, ys, zs)
	zs[:] = randn!(zs) ./ sqrt(o.vs) + o.thetas
	o.noise[:] = randn!(o.noise) ./ sqrt(o.vs)
	xs[:] = zs + o.noise
	ys[:] = zs - o.noise
end

#Default implementation
#Leaves cs, thetas and vs fixed in o
sim!(o, x, y, z) = simData!(o, x, y, z)

#Default experiments fix a single theta, cs, vs across all simulation
type DefaultExp
	cs
	thetas
	vs
	noise
end

#cs and vs fixed across all runs
function OddEvenExp(seed::Integer, n_max::Integer, odd_theta, odd_tau, frac_fit=.1)
	Random.seed!(seed)
	cs = 1./frac_fit * ones(n_max)

	thetas = ones(n_max)
	thetas[1:2:n_max] = odd_theta
	vs = ones(n_max)
	vs[1:2:n_max] = odd_tau

	noise = zeros(n_max)
	DefaultExp(cs, thetas, vs, noise)
end

function threePartExp(seed::Integer, n_max::Integer, theta_l, v_l, theta_h, v_h, frac_fit = .1)
	Random.seed!(seed)
	cs = 1./frac_fit * ones(n_max)
	thetas = ones(n_max)
	vs = ones(n_max)
	thetas[1:3:n_max] = theta_l
	vs[1:3:n_max] = v_l
	thetas[3:3:n_max] = theta_h
	vs[3:3:n_max] = v_h

	noise = zeros(n_max)
	DefaultExp(cs, thetas, vs, noise)
end

#Following simulates N rvs according to the given distributions
#AS N increases, this should approximately have gaussian likelihood
#dists should be mean zero.
type CLTExp
	cs
	vs
	thetas
	dists  #should be a mean zero r.v.
	N
end

#thetas /vs follow the three pt experiment
#width calculated accordingly
#Obsservations are uniform, centered on theta with width
function CLTExp(seed::Integer, width_max, n_max, N; frac_fit = .1)
	Random.seed!(seed)
	@assert rem(N, 2) == 0
	cs = 1./frac_fit * ones(n_max)
	thetas = ones(n_max)
	vs = ones(n_max)
	thetas[1:3:n_max] = .01
	vs[1:3:n_max] = .01
	thetas[3:3:n_max] = 1.5
	vs[3:3:n_max] = .1

	widths = sqrt(12 ./ vs)

	# widths = rand(n_max) * width_max + .1
	# vs = 12 ./ widths.^2
	# thetas = randn(n_max) ./ sqrt(tau0)

	dists = Array(Distributions.Uniform, n_max)
	for ix = 1:n_max
		dists[ix] = Uniform(-widths[ix]/2, widths[ix]/2 )
	end
	CLTExp(cs, vs, thetas, dists, N)
end

function sim!(o::CLTExp, xs, ys, zs)
	const n_max = length(o.cs)
	zetas = zeros(o.N)
	const half_N = round(Int, o.N/2)
	for ix = 1:n_max
		rand!(o.dists[ix], zetas)
		xs[ix] = mean(zetas[1:half_N]) * sqrt(half_N) + o.thetas[ix]
		ys[ix] = mean(zetas[(half_N + 1):o.N]) * sqrt(half_N) + o.thetas[ix]
		zs[ix] = mean(zetas) * sqrt(o.N) + o.thetas[ix]
	end
end

## The following constitutes a counter-example to tauxy/2 policy
type UniformLikelihoodExp
	cs
	vs
	thetas
	dist_odd
	dist_even
end

#even numbers have higher mean, lower std. dev
function UniformLikelihoodExp(seed::Integer, theta_even, theta_odd, width_even, width_odd, n_max, lamp)
	Random.seed!(seed)
	@assert rem(n_max, 2) == 0

	#compute the value of c and populate
	@assert theta_even + width_even < lamp < theta_odd + width_odd
	c = 4 * width_odd / (theta_odd + width_odd - lamp)
	println("C Value:\t", c)
	cs = c * ones(n_max)

	#check that we've picked values corresponding to the "interesting" case
	@assert theta_odd + width_odd > theta_even + width_even
	thetas = theta_odd * ones(n_max)
	thetas[2:2:n_max] = theta_even
	widths = width_odd * ones(n_max)
	widths[2:2:n_max] = width_even

	#create the distributions for later
	dist_odd = Uniform(theta_odd - width_odd, theta_odd + width_odd)
	dist_even = Uniform(theta_even - width_even, theta_even + width_even)
	UniformLikelihoodExp(cs, widths.^(-2),  thetas, dist_odd, dist_even)
end

##Specialized implementation for non-gaussian likelihood
function sim!(o::UniformLikelihoodExp, xs, ys, zs)
	n_max = length(o.cs)
	half_n_max = round(Int, n_max/2)

	xs[1:2:n_max] = rand(o.dist_odd, half_n_max)
	xs[2:2:n_max] = rand(o.dist_even, half_n_max)
	ys[1:2:n_max] = rand(o.dist_odd, half_n_max)
	ys[2:2:n_max] = rand(o.dist_even, half_n_max)
	zs[:] = .5 * (xs + ys);
end

#cs and vs fixed across all runs
#thetas drawn as gaussians tau0
type NormalBayesExp
	cs
	vs
	thetas
	noise
	tau0
	mu0
end

function NormalBayesExp(seed::Integer, tau0, n_max; mu0 = 0., frac_fit=.1, avg_tau=tau0)
	Random.seed!(seed)
	cs = rand(n_max) * 2./frac_fit
	vs = 2. * avg_tau * rand(n_max)
	NormalBayesExp(cs, vs, zeros(n_max), zeros(n_max), tau0, mu0)
end

#cs vs remain fixed but thetas change
function sim!(o::NormalBayesExp, x, y, z)
	randn!(o.thetas)
	o.thetas /= sqrt(o.tau0)
	o.thetas += o.mu0
	simData!(o, x, y, z)
end

#cs and vs fixed across all runs
#thetas drawn as gamma(alpha, beta)
type GammaExp
	cs
	thetas
	vs
	noise
	alpha
	beta
end

#thetas follow an exponential distribution
#Costs are uniform
function GammaExp(seed::Integer, alpha, beta, n_max, frac_fit = .1, avg_tau=2)
	Random.seed!(seed)
	cs = rand(n_max) * 2./frac_fit
	vs = 2. * avg_tau * rand(n_max)
	GammaExp(cs, zeros(Float64, n_max), vs, zeros(Float64, n_max), alpha, beta)
end

function sim!(o::GammaExp, x, y, z)
	#cs & vs are fixed
	rand!(Gamma(o.alpha, o.beta), o.thetas)
	simData!(o, x, y, z)
end

#cs and vs fixed across all runs
#thetas drawn as uniform [a, b]
type UniformExp
	cs
	thetas
	vs
	noise
	a
	b
end

function UniformExp(seed::Integer, a, b, n_max, frac_fit = .1, avg_vs=2)
	Random.seed!(seed)
	cs = rand(n_max) * 2./frac_fit
	vs = 2. * avg_vs * rand(n_max)
	UniformExp(cs, zeros(Float64, n_max), vs, zeros(Float64, n_max), a, b)
end

function sim!(o::UniformExp, x, y, z)
	#cs & vs are fixed
	rand!(Uniform(o.a, o.b), o.thetas)
	simData!(o, x, y, z)
end

#cs and vs fixed across all runs
#thetas drawn as beta[a, b]
type BetaExp
	cs
	thetas
	vs
	noise
	a
	b
end

function BetaExp(seed::Integer, a, b, n_max, frac_fit = .1, avg_tau=2)
	Random.seed!(seed)
	cs = rand(n_max) * 2./frac_fit
	vs = 2. * avg_tau * rand(n_max)
	BetaExp(cs, zeros(Float64, n_max), vs, zeros(Float64, n_max), a, b)
end

function sim!(o::BetaExp, x, y, z)
	#cs & vs are fixed
	rand!(Beta(o.a, o.b), o.thetas)
	simData!(o, x, y, z)
end

######################
####
#For increasing n, many simulations compare
#		#SAA/ Naive Z
#		#FullInfo Zstar
#		#Our method X
#		#Our method Avg X Y
#		#Oracle x method
#		#Oracle Z method
#		#Best randomized solution with n^k choices k = 2, 4, 6
#Compare on metrics
#		#Value wrt theta
#		#Time to compute
#		#optimal value of tau0
function test_harness(f, numRuns, o, n_grid)
	const n_max = maximum(n_grid)
	xs = zeros(Float64, n_max)
	ys = zeros(Float64, n_max)
	zs = zeros(Float64, n_max)

	#write a header
	writedlm(f,  ["Run" "n" "Method" "thetaVal" "time" "tau0"], ',')

	for iRun = 1:numRuns
		#generate the entire path up to n_max
		sim!(o, xs, ys, zs)

		for n in n_grid
			#Compute performance of each method

			#Naive Z
			#Same as the ZZ/SAA
			tic()
			qs = x(o.cs[1:n], zs[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "SAA" thetaval t 0.], ',')

			#fullInfo val
			tic()
			qs = x(o.cs[1:n], o.thetas[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "FullInfo" thetaval t 0.], ',')

			#The "Bayes" value.. only possible bc we if we we are in bayesian
			if :tau0 in fieldnames(o)
				tic()
				qs = x(o.cs[1:n], shrink(zs[1:n], o.vs[1:n], o.tau0))
				t = toc()
				thetaval = dot(o.thetas[1:n], qs)/n
				writedlm(f,  [iRun n "Bayes" thetaval t o.tau0], ',')
			end

			#Tau MLE
			tic()
			tauMLE, qs = x_MLE(o.cs[1:n], zs[1:n], o.vs[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "MLE" thetaval t tauMLE], ',')

			#Tau MM
			tic()
			tauMM, qs = x_MM(o.cs[1:n], zs[1:n], o.vs[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "MM" thetaval t tauMM], ',')

			#The half-sample heuristic
			tic()
			qs, vals, objs = best_x_tau(o.cs[1:n], xs[1:n], o.vs[1:n], ys[1:n])
			t = toc()

			tau_RS = vals[argmax(objs)]/2
			qs = x(o.cs[1:n], shrink(zs[1:n], o.vs[1:n], tau_RS))
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "Rescaled" thetaval t tau_RS], ',')

			#The idealized stein approach
			tic()
			qs, vals, objs = KP.x_stein_exact(o.cs[1:n], zs[1:n], o.vs[1:n], o.thetas[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "ExactStein" thetaval t vals[argmax(objs)]], ',')

			#The stein approach with various kernels
			#Box with the optimized rate, i.e. h_n = n^-1/6 and scaling, altKernel
			h = n^-.16666
			tic()
			qs, vals, objs = KP.x_stein_box(o.cs[1:n], zs[1:n], o.vs[1:n], h, KP.box, tau_step = .05)
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n_max
			writedlm(f,  [iRun n "Box" thetaval t vals[argmax(objs)]], ',')

			#gauss with MISE rate n^-1/5 with scaling, altKernel
			h = n^-.2
			tic()
			qs, vals, objs = KP.x_stein_box(o.cs[1:n], zs[1:n], o.vs[1:n], h, KP.gauss, tau_step=.05)
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "Gauss" thetaval t vals[argmax(objs)]], ',')

			#gauss with optimized rate n^-1/6 with scaling, altKernel
			h = n^-.16666
			tic()
			qs, vals, objs = KP.x_stein_box(o.cs[1:n], zs[1:n], o.vs[1:n], h, KP.gauss, tau_step=.05)
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "GaussO" thetaval t vals[argmax(objs)]], ',')

			#sinc with 1/log(n) with scaling, altKernel
			h = 1/log(n + 1)
			tic()
			qs, vals, objs = KP.x_stein_box(o.cs[1:n], zs[1:n], o.vs[1:n], h, sinc, tau_step = .1)
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "Sinc" thetaval t vals[argmax(objs)]], ',')

			#sinc with 1/log(n) with scaling, altKernel
			h = n^-.166666
			tic()
			qs, vals, objs = KP.x_stein_box(o.cs[1:n], zs[1:n], o.vs[1:n], h, sinc, tau_step = .1)
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "SincO" thetaval t vals[argmax(objs)]], ',')

			#Shrinkage via Cross-Val ("Ridge")
			tic()
			qs, tau_CV = x_HO_MSE(o.cs[1:n], xs[1:n], ys[1:n], o.vs[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "CV_Shrink" thetaval t tau_CV], ',')

			#weighted l2 regularization.  uses the oracle value for now
			tic()
			qs, mu = x_l2reg_CV(o.cs[1:n], zs[1:n], o.vs[1:n], o.thetas[1:n])[1:2]
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "OracleReg" thetaval t mu], ',')

			#weighted l2 regularization.  uses Cross-val to find value and then resolves
			tic()
			qs, mu = x_l2reg_CV(o.cs[1:n], xs[1:n], o.vs[1:n], ys[1:n])[1:2]
			t = toc()
			thetaval = dot(o.thetas[1:n], x_l2reg(o.cs[1:n], zs[1:n], o.vs[1:n], mu)[1])/n
			writedlm(f,  [iRun n "RegCV" thetaval t mu], ',')

			#Plug-in SURE estimation for L2
			tic()
			qs, tau = x_sure_MSE(o.cs[1:n], zs[1:n], o.vs[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			writedlm(f,  [iRun n "SURE" thetaval t tau], ',')

			# #The primal stein approach
			# #use the optimized rate, i.e. h_n = n^-1/6
			# h = n^-.16666
			# tic()
			# qs, vals, objs = KP.x_stein_primal(o.cs[1:n], zs[1:n], o.vs[1:n], h)
			# t = toc()
			# yval = dot(ys[1:n], qs)/n
			# thetaval = dot(o.thetas[1:n], qs)/n
			# writedlm(f,  [iRun n "PrimalStein" yval thetaval t vals[argmax(objs)]], ',')

			#Oracle Value
			tic()
			qs, vals, objs = best_x_tau(o.cs[1:n], zs[1:n], o.vs[1:n], o.thetas[1:n])
			t = toc()
			thetaval = dot(o.thetas[1:n], qs)/n
			if abs(thetaval - maximum(objs)) > 1e-5
				f2 = open("Mismatch_log.csv", "w")
				writecsv(f2, [n iRun])
				writecsv(f2, xs')
				writecsv(f2, ys')
				writecsv(f2, zs')
				writecsv(f2, o.vs[1:n]')
				writecsv(f2, o.thetas[1:n]')
				writecsv(f2, o.cs[1:n]')
				close(f2)
				exit()
			end

			@assert abs(thetaval - maximum(objs)) <= 1e-5 "Weird Mismatch? \t $thetaval \t $(maximum(objs))"
			writedlm(f,  [iRun n "OracleZ" thetaval t vals[argmax(objs)]], ',')

		end
		flush(f)
	end
end

#######################

### Bayesian set-up with a non-zero mean
### costs and vs are random, but fixed over runs
function test_Gaussian(file_out, numRuns, n_grid, seed, tau0, mu0, avg_tau)
	#build the sim object
	o = NormalBayesExp(seed, tau0, maximum(n_grid), mu0=mu0, avg_tau=avg_tau)

	#run the testharness
	#output files
	f = open("$(file_out)_$(tau0)_$(mu0)_$(avg_tau)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "$(file_out)_$(tau0)_$(mu0)_$(avg_tau)_$(seed).csv"
end

### The odd even set-up.
function test_OddEven(file_out, numRuns, n_grid, seed, odd_theta, odd_tau)
	#build the sim object
	o = OddEvenExp(seed, maximum(n_grid), odd_theta, odd_tau)

	#run the testharness
	#output files
	f = open("$(file_out)_$(odd_theta)_$(odd_tau)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "$(file_out)_$(odd_theta)_$(odd_tau)_$(seed).csv"
end

### The three parts experimental set-up
function test_threePart(file_out, numRuns, n_grid, seed, theta_l, v_l, theta_h, v_h)
	o = threePartExp(seed, maximum(n_grid), theta_l, v_l, theta_h, v_h)
	f = open("$(file_out)_$(theta_l)_$(v_l)_$(theta_h)_$(v_h)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "$(file_out)_$(theta_l)_$(v_l)_$(theta_h)_$(v_h)_$(seed).csv"
end

### The Gamma Test
function test_Gamma(file_out, numRuns, n_grid, seed, alpha, beta)
	o = GammaExp(seed, alpha, beta, maximum(n_grid))

	f = open("$(file_out)_Gamma_$(alpha)_$(beta)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "$(file_out)_Gamma_$(alpha)_$(beta)_$(seed).csv"
end

### The Uniform Test
function test_Uniform(file_out, numRuns, n_grid, seed, a, b)
	o = UniformExp(seed, a, b, maximum(n_grid))

	f = open("$(file_out)_Uniform_$(a)_$(b)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "$(file_out)_Uniform_$(a)_$(b)_$(seed).csv"
end

### The Beta Test
function test_Beta(file_out, numRuns, n_grid, seed, a, b)
	o = BetaExp(seed, a, b, maximum(n_grid))

	f = open("$(file_out)_Beta_$(a)_$(b)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "$(file_out)_Beta_$(a)_$(b)_$(seed).csv"
end

### The Beta Test
function test_UniformLikelihood(numRuns, n_grid, seed, theta_high, width_low, width_high, lamp)
	o = UniformLikelihoodExp(seed, theta_high, 1., width_high, width_low, maximum(n_grid), lamp)
	f = open("UniformLikelihood_$(theta_high)_$(width_high)_$(width_low)_$(lamp)_$(seed).csv", "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return "UniformLikelihood_$(theta_high)_$(width_high)_$(width_low)_$(lamp)_$(seed).csv"
end

function test_CLTExp(file_out, numRuns, n_grid, seed, N, width_max)
	o = CLTExp(seed, width_max, maximum(n_grid), N)
	tag = "$(file_out)CLTExp_$(N)_$(width_max)_$(seed).csv"
	f = open(tag, "w")
	test_harness(f, numRuns, o, n_grid)
	close(f)
	return tag
end

#######
# A by-hand test for bandwidths and reg parameters
# o is assumed pre-initialized
function test_bandwidth(file_out, numRuns, o)
	const n = length(o.thetas)
	xs = zeros(Float64, n)
	ys = zeros(Float64, n)
	zs = zeros(Float64, n)

	f = open(file_out, "w")
	writedlm(f,  ["Run" "n" "Method" "bandwidth" "thetaVal" "tau0"], ',')
	bandwidths = linspace(-1, -.001, 30)
	#regs = linspace(.05, 5, 15)

	for iRun = 1:numRuns
		sim!(o, xs, ys, zs)

		for h in bandwidths
			qs, vals, objs = KP.x_stein_box(o.cs, zs, o.vs, n^h, KP.box, tau_step = .1)
			thetaval = dot(o.thetas, qs)/n
			writedlm(f,  [iRun n "Box" h thetaval vals[argmax(objs)]], ',')

			qs, vals, objs = KP.x_stein_box(o.cs, zs, o.vs, n^h, KP.gauss, tau_step = .1)
			thetaval = dot(o.thetas, qs)/n
			writedlm(f,  [iRun n "Gauss" h thetaval vals[argmax(objs)]], ',')

			qs, vals, objs = KP.x_stein_box(o.cs, zs, o.vs, n^h, sinc, tau_step = .1)
			thetaval = dot(o.thetas, qs)/n
			writedlm(f,  [iRun n "Sinc" h thetaval vals[argmax(objs)]], ',')
		end
	end
	close(f)
end

#########
n_grid = [2^i for i = 8:17]

#run some small examples to pre-compile for optimization
test_threePart("./temp/tempThreePart", 5, [100, 150], 876, .01, .01, 1.5, .1)
test_CLTExp("./temp/tempCLT", 5, [100, 150], 86, 10, 2)
# test_Gaussian("./temp/temp_Gaussian", 5, [100, 150], 87, 3, 1, 3)
# test_OddEven("./temp/temp_OddEven", 5,[100, 150], 87, 2, 2)

# test_Gamma("./temp/temp_Gamma", 5, [100, 150], 87, 1., 1.)
# test_Uniform("./temp/temp_Uniform", 5, [100, 150], 87, 1, 2)
# test_Beta("./temp/temp_Beta", 5, [100, 150], 87, .5, .5)
# test_bandwidth("./temp/tempbandwidth", 10, NormalBayesExp(8675309, 3, 100))

# #The counterexample
# n_grid = [2^i for i = 8:20]
# test_UniformLikelihood(2, [100, 150], 8675309, 10., 15., 1., 14.1)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Documenter, FMI

makedocs(sitename="FMI.jl",
          pages= Any[
                "Home" => "index.md"
                "Examples" => [
                    "Parameterize" => "parameterize.md"
                    "Simulate CoSimulation" => "simulateCS.md"
                    "Simulate ModelExchange" => "simulateME.md"
                ]
                "Library Functions" => Any[
                    "FMI library functions" => "library.md",
                    "FMU2 functions" => "fmu2.md"
                ]
         ]
         )

deploydocs(repo = "github.com/ThummeTo/FMI.jl.git", devbranch = "main")


text:
using Plots
sphere(x) = sum(x.^2)

dpi = 600
x = y = range(-3, stop=3, length=200)

l = Vector{Float64}()
for i in x
   for j in y
       v = [i j]
       append!(l, sphere(v))
   end
end

surface(x,y,l, dpi=dpi)
png("./sphere_surface.png")
contour(x,y,l, dpi=dpi, fill = true)
png("./sphere_contour.png")


text:
module SpectralLibraries

abstract type AbstractSpectralLibrary end

include("phoenix.jl")

end # module


text:
"""
problem.jl
"""

"""$(TYPEDEF) Trajectory Optimization Problem.
Contains the full definition of a trajectory optimization problem, including:
* dynamics model (`Model`)
* objective (`Objective`)
* constraints (`ConstraintSet`)
* initial and final states
* Primal variables (state and control trajectories)
* Discretization information: knot points (`N`), time step (`dt`), and total time (`tf`)

# Constructors:
```julia
Problem(model, obj, constraints, x0, xf, Z, N, tf) # defaults to RK3 integration
Problem{Q}(model, obj, constraints, x0, xf, Z, N, tf) where Q<:QuadratureRule
Problem(model, obj, xf, tf; x0, constraints, N, X0, U0, dt, integration)
Problem{Q}(prob::Problem)  # change integration
```
where `Z` is a trajectory (Vector of `KnotPoint`s)

# Arguments
* `model`: Dynamics model. Can be either `Discrete` or `Continuous`
* `obj`: Objective
* `X0`: Initial state trajectory. If omitted it will be initialized with NaNs, to be later overwritten by the solver.
* `U0`: Initial control trajectory. If omitted it will be initialized with zeros.
* `x0`: Initial state. Defaults to zeros.
* `xf`: Final state. Defaults to zeros.
* `dt`: Time step
* `tf`: Final time. Set to zero to specify a time penalized problem.
* `N`: Number of knot points. Defaults to 51, unless specified by `dt` and `tf`.
* `integration`: One of the defined integration types to discretize the continuous dynamics model.
Both `X0` and `U0` can be either a `Matrix` or a `Vector{Vector}`, but must be the same.
At least 2 of `dt`, `tf`, and `N` need to be specified (or just 1 of `dt` and `tf`).
"""
struct Problem{IR,T,Tm,To,Tx,Tix,Tu,Tiu,Tt,TE,TM,TMd,TV}
    # problem info
    n::Int
    m::Int
    N::Int
    model::Tm
    obj::To
    convals::Vector{Vector{ConVal}}
    X::Vector{Tx}
    X_tmp::Vector{Tx}
    ix::Tix
    U::Vector{Tu}
    U_tmp::Vector{Tu}
    iu::Tiu
    ts::Tt
    E::TE
    M::TM
    Md::TMd
    V::TV
end

function Problem(::Type{IR}, model::Tm, obj::To, constraints::ConstraintList,
                 X::Vector{Tx}, U::Vector{Tu}, ts::Tt, N::Int, M::TM, Md::TMd, V::TV) where {
                     IR<:QuadratureRule,Tm<:AbstractModel,To<:AbstractObjective,
                     Tx<:AbstractVector,Tu<:AbstractVector,Tt<:AbstractVector,
                     TM,TMd,TV}
    n, m = size(model)
    # allocate shared resources
    T = eltype(X[1])
    X_tmp = [V(zeros(T, n)) for k = 1:N+1] # 1 extra
    U_tmp = [V(zeros(T, m)) for k = 1:N+1] # 2 extra
    E = QuadraticCost(M(zeros(T, n, n)), M(zeros(T, m, m)), M(zeros(T, m, n)), V(zeros(T, n)),
                      V(zeros(T, m)), zero(T); checks=false)
    # construct indices into concatenated state and controls
    ix = V(1:n)
    iu = V((1:m) .+ n)
    # initial condition constraint for direct solve
    initial_state_constraint = GoalConstraint(copy(X[1]), V(1:n), n, m, M, V; direct=true)
    add_constraint!(constraints, initial_state_constraint, 1:1)
    # dynamics constraint for direct solve
    dynamics_constraint = DynamicsConstraint(IR, model, ts, ix, iu, n, m, M, V)
    add_constraint!(constraints, dynamics_constraint, 2:N-1)
    # create convals from constraints
    convals = convals_from_constraint_list(constraints)
    # put it all together
    Tix = typeof(ix)
    Tiu = typeof(iu)
    T = eltype(X[1])
    TE = typeof(E)
    return Problem{IR,T,Tm,To,Tx,Tix,Tu,Tiu,Tt,TE,TM,TMd,TV}(
        n, m, N, model, obj, convals, X, X_tmp, ix, U, U_tmp, iu, ts, E, M, Md, V
    )
end

"$(TYPEDSIGNATURES)
Get number of states, controls, and knot points"
Base.size(prob::Problem) = size(prob.model)..., prob.N


text:
function main()
    N, K = parse.(Int, split(readline()))

    f(N, K) = min(K-1, 2N + 1 - K) * (2 ≤ K ≤ 2N)

    ans = 0
    for x in 2:2N
        ans += f(N, x) * f(N, x - K)
    end

    println(ans)

end

main()


text:
### Serial "ordinary" Monte Carlo runner

immutable SerialMCBaseRunner <: SerialMCRunner
  burnin::Int
  thinning::Int
  nsteps::Int
  r::Range{Int}
  storegradlogtarget::Bool # Indicates whether to save the gradient of the log-target in the cases it is available

  function SerialMCBaseRunner(r::Range{Int}, s::Bool=false)
    burnin = first(r)-1
    thinning = r.step
    nsteps = last(r)
    @assert burnin >= 0 "Number of burn-in iterations should be non-negative."
    @assert thinning >= 1 "Thinning should be >= 1."
    @assert nsteps > burnin "Total number of MCMC iterations should be greater than number of burn-in iterations."
    new(burnin, thinning, nsteps, r, s)
  end
end

SerialMCBaseRunner(r::UnitRange{Int}, s::Bool=false) = SerialMCBaseRunner(first(r):1:last(r), s)

SerialMCBaseRunner(; burnin::Int=0, thinning::Int=1, nsteps::Int=100, storegradlogtarget::Bool=false) =
  SerialMCBaseRunner((burnin+1):thinning:nsteps, storegradlogtarget)

typealias SerialMC SerialMCBaseRunner

function run(m::MCModel, s::MCSampler, r::SerialMC, t::MCTuner=VanillaMCTuner(), job::Symbol=:task)
  tic()

  mcjob = MCJob(m, s, r; tuner=t, job=job)

  # Pre-allocation for storing results
  mcchain::MCChain = MCChain(m.size, length(r.r); storegradlogtarget=r.storegradlogtarget)
  ds = Dict{Any, Any}("step" => collect(r.r))

  # Sampling loop
  i::Int = 1
  for j in 1:r.nsteps
    mcstate = mcjob.receive(1)
    if in(j, r.r)
      mcchain.samples[i, :] = mcstate.successive.sample
      mcchain.logtargets[i] = mcstate.successive.logtarget

      if r.storegradlogtarget
        mcchain.gradlogtargets[i, :] = mcstate.successive.gradlogtarget
      end

      # Save diagnostics
      for (k,v) in mcstate.diagnostics
        # If diagnostics name not seen before, create column
        if !haskey(ds, k)
          ds[k] = Array(typeof(v), length(ds["step"]))
        end

        ds[k][i] = v
      end

      i += 1
    end
  end

  mcchain.diagnostics, mcchain.runtime = ds, toq()
  mcchain
end

function resume!(m::MCModel, s::MCSampler, r::SerialMC, c::MCChain, t::MCTuner=VanillaMCTuner(), j::Symbol=:task;
  nsteps::Int=100)
  m.init = vec(c.samples[end, :])
  mcrunner::SerialMC = SerialMC(burnin=0, thinning=r.thinning, nsteps=nsteps, storegradlogtarget=r.storegradlogtarget)
  run(m, s, mcrunner, t, j)
end

resume(m::MCModel, s::MCSampler, r::SerialMC, c::MCChain, t::MCTuner=VanillaMCTuner(), j::Symbol=:task;
  nsteps::Int=100) =
  resume!(deepcopy(m), s, r, c, t, j; nsteps=nsteps)


text:
module BoundaryTriangulationsTests

using Gridap

model = CartesianDiscreteModel(partition=(3,4,2))

tags = [23,24,25]
grid = BoundaryGrid(model,tags)

trian = Triangulation(grid)
descr = grid.descriptor

btrian = BoundaryTriangulation(trian,descr)
test_triangulation(btrian)

tags = [23,24,25]
btrian = BoundaryTriangulation(model,tags)

tags = "boundary"
btrian = BoundaryTriangulation(model,tags)

tags = ["physical_tag_23","physical_tag_24","physical_tag_25"]
btrian = BoundaryTriangulation(model,tags)

btrian = BoundaryTriangulation(model)

end # module


text:
"""
    range32(start, stop, step)

StepRangeLen but for Float32.
"""
function range32(start::Real, stop::Real, step::Real)
  res = range(Float32(start), stop=Float32(stop), step = Float32(step))
  return StepRangeLen{Float32,Float32,Float32}(res)
end

max_stars_repo_name: williamjsdavis/ThermoChroNN

max_issues_repo_name: williamjsdavis/ThermoChroNN

max_forks_repo_name: williamjsdavis/ThermoChroNN


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################
"""
Example considering a multi-state Markov model for a coal fired power generating
unit taken from:

> A Multi-State Markov Model for a Short-Term Reliability Analysis of a Power
  Generating Unit by A. Lisnianski, D. Elmakias, and H. Ben Haim (2012)

A unit's available generating capacity can take any real value in [0,gⁿᵒᵐ].
This is approximated by a discrete state continuous time stochastic process
Gₐ(t) as follows:
- Two states to describe the limits of the generating capacity:
    - State 1, where g₁ = 0
    - State N, where gₙ = gⁿᵒᵐ
- N-2 additional states seperated with an 'averaged' output gᵢ within an
  interval Δg = gⁿᵒᵐ/(N-2).

If the random times between transitions are ignored and only the transition time
instances are of interest, the stochastic process Gₐ(t) reduces to a Markov
model Gₐₘ(t), with approximated transition intensities âᵢⱼ from state i to
state j≠i.

The case study in the paper considers a four-state model to represent a coal
fired power generating unit, with the input given in Table 1[^1].

Table 1: States' output and transition intensities.
| g₁  = 0 MW        | g₂  = 247 MW      | g₃  = 482 MW      | g₄  = 576 MW      |
|-------------------|-------------------|-------------------|-------------------|
| â₁₂ ≈ 0.08000 1/h | â₂₁ ≈ 0.02912 1/h | â₃₁ ≈ 0.00000 1/h | â₄₁ ≈ 0.00015 1/h |
| â₁₃ ≈ 0.01333 1/h | â₂₃ ≈ 0.32353 1/h | â₃₂ ≈ 0.02885 1/h | â₄₂ ≈ 0.00010 1/h |
| â₁₄ ≈ 0.00000 1/h | â₂₄ ≈ 0.02912 1/h | â₃₄ ≈ 0.35577 1/h | â₄₃ ≈ 0.00069 1/h |

[^1]: The data in the table is an approximation. Exact results are obtained by
using the fraction of kᵢⱼ over Tᵢ as given in the paper.
"""

# load pkgs
using Unitful
using MultiStateSystems

# initialize the state-transition diagram corresponding to the coal fired power
# generating unit.
stdᶜᶠ = STD()

# add the states to the std
add_states!(stdᶜᶠ, power = [0u"MW", 247u"MW", 482u"MW", 576u"MW"],
                   init  = [0.0, 0.0, 0.0, 1.0])

# add the transitions to the std
add_transitions!(stdᶜᶠ, states = [(1,2),(1,3)],
                        rate = [0.08000u"1/hr",0.01333u"1/hr"])
add_transitions!(stdᶜᶠ, states = [(2,1),(2,3),(2,4)],
                        rate = [0.02941u"1/hr",0.32353u"1/hr",0.02941u"1/hr"])
add_transitions!(stdᶜᶠ, states = [(3,2),(3,4)],
                        rate = [0.02885u"1/hr",0.35577u"1/hr"])
add_transitions!(stdᶜᶠ, states = [(4,1),(4,2),(4,3)],
                        rate = [0.00015u"1/hr",0.00010u"1/hr",0.00069u"1/hr"])

# return the std
return stdᶜᶠ


text:
using Redis
import DataStructures: OrderedSet

if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

include(joinpath(dirname(@__FILE__),"client_tests.jl"))
include(joinpath(dirname(@__FILE__),"redis_tests.jl"))


text:
module SorbetHelperLinkedGateBlockActivatorTrigger

using ..Ahorn, Maple

@mapdef Trigger "SorbetHelper/LinkedGateBlockActivatorTrigger" LinkedGateBlockActivatorTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, linkTag::String="", flag::String="", inverted::Bool=false)

const placements = Ahorn.PlacementDict(
    "Linked Gate Block Activator (Sorbet Helper)" => Ahorn.EntityPlacement(
        LinkedGateBlockActivatorTrigger,
        "rectangle"
    )
)

Ahorn.editingOrder(entity::LinkedGateBlockActivatorTrigger) = String["x", "y", "width", "height", "linkTag", "flag", "inverted"]

end

max_stars_repo_name: Kitty-Cats/SorbetHelper

max_issues_repo_name: Kitty-Cats/SorbetHelper

max_forks_repo_name: Kitty-Cats/SorbetHelper


text:
using ParallelAccelerator

#ParallelAccelerator.ParallelIR.set_debug_level(3)

@acc function find_chg(k,m,W,Wp)
	W_tmp = [abs(W[i,j] - Wp[i,j]) for i in 1:m, j in 1:k]
	s = [sum(W_tmp[:,j]) for j in 1:k]
	chg = maximum(s)
end

function main(m::Int, k::Int)
    W  = Array{Float64}(m, k)
    Wp = Array{Float64}(m, k)
    fill!(W, 3)
    fill!(Wp, 5)
    chg = find_chg(k,m,W,Wp)

    println("chg = ", chg)

end

main(100, 100)


text:
using StatsBase
using Random, Test

Random.seed!(1234)

#### weighted sample with replacement

function check_wsample_wrep(a::AbstractArray, vrgn, wv::AbstractWeights, ptol::Real; ordered::Bool=false)
    K = length(wv)
    (vmin, vmax) = vrgn
    (amin, amax) = extrema(a)
    @test vmin <= amin <= amax <= vmax
    p0 = values(wv) ./ sum(wv)
    if ordered
        @test issorted(a)
        if ptol > 0
            @test isapprox(proportions(a, vmin:vmax), p0, atol=ptol)
        end
    else
        @test !issorted(a)
        ncols = size(a,2)
        if ncols == 1
            @test isapprox(proportions(a, vmin:vmax), p0, atol=ptol)
        else
            for j = 1:ncols
                aj = view(a, :, j)
                @test isapprox(proportions(aj, vmin:vmax), p0, atol=ptol)
            end
        end
    end
end

import StatsBase: direct_sample!, alias_sample!

n = 10^5
wv = weights([0.2, 0.8, 0.4, 0.6])

a = direct_sample!(4:7, wv, zeros(Int, n, 3))
check_wsample_wrep(a, (4, 7), wv, 5.0e-3; ordered=false)
test_rng_use(direct_sample!, 4:7, wv, zeros(Int, 100))

a = alias_sample!(4:7, wv, zeros(Int, n, 3))
check_wsample_wrep(a, (4, 7), wv, 5.0e-3; ordered=false)

a = sample(4:7, wv, n; ordered=false)
check_wsample_wrep(a, (4, 7), wv, 5.0e-3; ordered=false)

a = sample(4:7, wv, n; ordered=true)
check_wsample_wrep(a, (4, 7), wv, 5.0e-3; ordered=true)

#### weighted sampling without replacement

function check_wsample_norep(a::AbstractArray, vrgn, wv::AbstractWeights, ptol::Real; ordered::Bool=false)
    # each column of a for one run

    vmin, vmax = vrgn
    (amin, amax) = extrema(a)
    @test vmin <= amin <= amax <= vmax
    n = vmax - vmin + 1

    for j = 1:size(a,2)
        aj = view(a,:,j)
        @assert allunique(aj)
        if ordered
            @assert issorted(aj)
        end
    end

    if ptol > 0
        p0 = values(wv) ./ sum(wv)
        @test isapprox(proportions(a[1,:], vmin:vmax), p0, atol=ptol)
    end
end

import StatsBase: naive_wsample_norep!, efraimidis_a_wsample_norep!,
                  efraimidis_ares_wsample_norep!, efraimidis_aexpj_wsample_norep!

n = 10^5
wv = weights([0.2, 0.8, 0.4, 0.6])

a = zeros(Int, 3, n)
for j = 1:n
    naive_wsample_norep!(4:7, wv, view(a,:,j))
end
check_wsample_norep(a, (4, 7), wv, 5.0e-3; ordered=false)
test_rng_use(naive_wsample_norep!, 4:7, wv, zeros(Int, 2))

a = zeros(Int, 3, n)
for j = 1:n
    efraimidis_a_wsample_norep!(4:7, wv, view(a,:,j))
end
check_wsample_norep(a, (4, 7), wv, 5.0e-3; ordered=false)
test_rng_use(efraimidis_a_wsample_norep!, 4:7, wv, zeros(Int, 2))

a = zeros(Int, 3, n)
for j = 1:n
    efraimidis_ares_wsample_norep!(4:7, wv, view(a,:,j))
end
check_wsample_norep(a, (4, 7), wv, 5.0e-3; ordered=false)
test_rng_use(efraimidis_ares_wsample_norep!, 4:7, wv, zeros(Int, 2))

a = zeros(Int, 3, n)
for j = 1:n
    efraimidis_aexpj_wsample_norep!(4:7, wv, view(a,:,j))
end
check_wsample_norep(a, (4, 7), wv, 5.0e-3; ordered=false)
test_rng_use(efraimidis_aexpj_wsample_norep!, 4:7, wv, zeros(Int, 2))

a = sample(4:7, wv, 3; replace=false, ordered=false)
check_wsample_norep(a, (4, 7), wv, -1; ordered=false)

a = sample(4:7, wv, 3; replace=false, ordered=true)
check_wsample_norep(a, (4, 7), wv, -1; ordered=true)


text:
using Plots, ApproxFun

x = Fun()
D = Derivative()

L=D^10 + cosh(x)*D^8 + x^3*D^6 + x^4*D^4 + cos(x)*D^2 + x^2
d = ChebyshevInterval()
B = [Dirichlet(d) ;
     Neumann(d)   ;
     [Evaluation(d,leftendpoint,k) for k=2:4]... ;
     [Evaluation(d,rightendpoint,k) for k=2:4]...]

u = [B; L] \ [ [0.,0.], [1.,1.], zeros(6)..., exp(x)]

plot(u)

max_stars_repo_path: ODEs/Tenth order ODE.jl

max_issues_repo_path: ODEs/Tenth order ODE.jl

max_forks_repo_path: ODEs/Tenth order ODE.jl


text:
#Switch argument order? Convention is buffer first, but that conflicts with julia praxis
#This is not be identical to vi(m)
#Symmetric words don't support unicode, we have to do this another way
const Word = r"\S+|^\s+$|^$"
const word = r"[\pL\pN_]+|[\pS\pP]+|^\s+$|^$"

crs(s::String, r::Regex) = length(matchall(r, s))
crs(b::Buffer, r::Regex, n::Integer) = crs(line(b, n), r)

#It should be possible to make this faster
#it may be necessary, num arg word movement should be very fast
function next_pos(s::String, r::Regex, n=1)
    c = collect(eachmatch(r, s))
    if length(c)>=n
        if length(s)>1
            return ind2chr(s, c[n].offset)
        else
            return 1
        end
    else
        return -1
    end
end

next_pos(b::Buffer, r::Regex, n=1) = next_pos(unirange(line(b), max(x(b),1):width(b)),
                                              r, n)
function next_pos_naive(b::Buffer, r::Regex, n=1)
    s = unirange(line(b), max(x(b),1):width(b))
    if next_pos(s, r, 1) == 1
        n += 1
    end
    m = next_pos(s, r, n)
    if m!=-1
        return pos(b).+[0,m-1]
    else
        lf = length(collect(eachmatch(r, s)))+1
        wc = cumsum(crs.(b.text[y(b)+1:end], r))
        yoffset = findfirst(x->x>n, lf+wc)
        mc = n - lf
        if yoffset>1
            mc -= wc[yoffset-1]
        end
        if mc==0
            xoffset = next_pos(line(b, y(b)+yoffset), r, 1)
        else
            #Pretty clumsy solution
            if next_pos(line(b, y(b)+yoffset), r, 1)!=1
                xoffset = next_pos(line(b, y(b)+yoffset), r, mc)
            else
                xoffset = next_pos(line(b, y(b)+yoffset), r, mc+1)
            end
        end
        return [y(b)+yoffset, xoffset]
    end
end

function prev_pos_naive(b::Buffer, r::Regex, n=1)
    s = unirange(line(b), min(x(b),width(b)):-1:1)
    if next_pos(s, r, 1) == 1
        n += 1
    end
    m = next_pos(s, r, n)
    if m!=-1
        return pos(b).-[0,m-1]
    else
        lf = length(collect(eachmatch(r, s)))+1
        wc = cumsum(crs.(b.text[y(b)-1:-1:1], r))
        yoffset = findfirst(x->x>n, lf+wc)
        mc = n - lf
        if yoffset>1
            mc -= wc[yoffset-1]
        end
        if mc==0
            xoffset = next_pos(reverse(line(b, y(b)-yoffset)), r, 1)
        else
            #Pretty clumsy solution
            if next_pos(line(b, y(b)+yoffset), r, 1)!=1
                xoffset = next_pos(reverse(line(b, y(b)-yoffset)), r, mc)
            else
                xoffset = next_pos(reverse(line(b, y(b)-yoffset)), r, mc+1)
            end
        end
        return [y(b)-yoffset, length(line(b, y(b)-yoffset)) - xoffset+1]
    end
end

function eom(s::String, r::Regex)
    match(r, s).offset + length(match(r, s).match) - 1
end
eom(b::Buffer, r::Regex) = eom(unirange(line(b), max(x(b),1):width(b)), r)
bom(b::Buffer, r::Regex) = -eom(unirange(line(b), min(x(b),width(b)):-1:1), r)

function move_to_nth(b::Buffer, r::Regex, n)
    b.cursor.pos .= next_pos_naive(b, r, n)
end

function back_to_nth(b::Buffer, r::Regex, n)
    b.cursor.pos .= prev_pos_naive(b, r, n)
end

function move_to_nth_eom(b::Buffer, r::Regex, n)
    if eom(b, r)>1
        b.cursor.pos[2] += eom(b, r) - 1
        n -= 1
    end
    b.cursor.pos .= next_pos_naive(b, r::Regex, n)
    b.cursor.pos[2] += eom(b, r) - 1
end

function back_to_nth_eom(b::Buffer, r::Regex, n)
    if -bom(b, r)>1
        b.cursor.pos[2] += bom(b, r) + 1
        n -= 1
    end
    b.cursor.pos .= prev_pos_naive(b, r::Regex, n)
    b.cursor.pos[2] += bom(b, r) + 1
end

function move_ops(r::Regex)
    crs_r(b::Buffer, n) = crs(b, r, n)
    crs_r(s::String) = crs(s, r)
    next_pos_r(s::String, n=1) =  next_pos(s, r, n)
    next_pos_naive_r(b::Buffer, n=1)  = next_pos_naive(b, r, n)
    prev_pos_naive_r(b::Buffer, n=1)  = prev_pos_naive(b, r, n)
    eom_r(o::Union{Buffer, String}) = eom(o, r)
    bom_r(o::Union{Buffer, String}) = bom(o, r)
    move_to_nth_r(b::Buffer, n=parse_n(b)) = move_to_nth(b, r, n)
    move_to_nth_eom_r(b::Buffer, n=parse_n(b)) = move_to_nth_eom(b, r, n)
    back_to_nth_r(b::Buffer, n=parse_n(b)) = back_to_nth(b, r, n)
    back_to_nth_eom_r(b::Buffer, n=parse_n(b)) = back_to_nth_eom(b, r, n)
    [crs_r,
     next_pos_r,
     next_pos_naive_r,
     prev_pos_naive_r,
     eom_r,
     bom_r,
     move_to_nth_r,
     move_to_nth_eom_r,
     back_to_nth_r,
     back_to_nth_eom_r,
    ]
end

(cwords,
 nextword_pos,
 nextword_pos_naive,
 prevword_pos_naive,
 eow,
 bow,
 move_to_nth_word,
 move_to_nth_eow,
 back_to_nth_word,
 back_to_nth_eow,
) = move_ops(word)

(cWords,
 nextWord_pos,
 nextWord_pos_naive,
 prevWord_pos_naive,
 eoW,
 boW,
 move_to_nth_Word,
 move_to_nth_eoW,
 back_to_nth_Word,
 back_to_nth_eoW,
) = move_ops(Word)


text:
include("attributes.jl")
include("executors.jl")
include("penaltyterm.jl")


text:
"""
    copywholeraster(source::AbstractDataset, dest::AbstractDataset; <keyword arguments>)

Copy all dataset raster data.

This function copies the complete raster contents of one dataset to another
similarly configured dataset. The source and destination dataset must have the
same number of bands, and the same width and height. The bands do not have to
have the same data type.

Currently the only `options` supported are : `\"INTERLEAVE=PIXEL\"` to
force pixel interleaved operation and `\"COMPRESSED=YES\"` to force alignment on
target dataset block sizes to achieve best compression. More options may be
supported in the future.

### Additional Remarks
This function is primarily intended to support implementation of driver
specific `createcopy()` functions. It implements efficient copying, in
particular \"chunking\" the copy in substantial blocks and, if appropriate,
performing the transfer in a pixel interleaved fashion.
"""
function copywholeraster(
        source::AbstractDataset,
        dest::AbstractDataset;
        options                 = StringList(C_NULL),
        progressfunc::Function  = GDAL.gdaldummyprogress,
        progressdata::Any       = C_NULL
    )
    result = GDAL.gdaldatasetcopywholeraster(source.ptr, dest.ptr, options,
        @cplprogress(progressfunc), progressdata)
    @cplerr result "Failed to copy whole raster"
end

"""
    unsafe_copy(dataset::AbstractDataset; [filename, [driver, [<keyword arguments>]]])

Create a copy of a dataset.

This method will attempt to create a copy of a raster dataset with the
indicated filename, and in this drivers format. Band number, size, type,
projection, geotransform and so forth are all to be copied from the
provided template dataset.

### Parameters
* `dataset`       the dataset being duplicated.

### Keyword Arguments
* `filename`      the filename for the new dataset. UTF-8 encoded.
* `driver`        the driver to use for creating the new dataset
* `strict`        ``true`` if the copy must be strictly equivalent, or more
normally ``false`` if the copy may adapt as needed for the output format.
* `options`       additional format dependent options controlling creation
of the output file. `The APPEND_SUBDATASET=YES` option can be specified to
avoid prior destruction of existing dataset.

### Returns
a pointer to the newly created dataset (may be read-only access).

### Additional Remarks
Note: many sequential write once formats (such as JPEG and PNG) don't implement
the `Create()` method but do implement this `CreateCopy()` method. If the
driver doesn't implement `CreateCopy()`, but does implement `Create()` then the
default `CreateCopy()` mechanism built on calling `Create()` will be used.

It is intended that `CreateCopy()` will often be used with a source dataset
which is a virtual dataset allowing configuration of band types, and other
information without actually duplicating raster data (see the VRT driver). This
is what is done by the gdal_translate utility for example.

This function will validate the creation option list passed to the driver
with the `GDALValidateCreationOptions()` method. This check can be disabled by
defining the configuration option `GDAL_VALIDATE_CREATION_OPTIONS=NO`.

After you have finished working with the returned dataset, it is required to
close it with `GDALClose()`. This does not only close the file handle, but also
ensures that all the data and metadata has been written to the dataset
(`GDALFlushCache()` is not sufficient for that purpose).

In some situations, the new dataset can be created in another process through
the GDAL API Proxy mechanism.
"""
function unsafe_copy(
        dataset::AbstractDataset;
        filename::AbstractString    = string("/vsimem/$(gensym())"),
        driver::Driver              = getdriver(dataset),
        strict::Bool                = false,
        options                     = StringList(C_NULL),
        progressfunc::Function      = GDAL.gdaldummyprogress,
        progressdata                = C_NULL
    )
    return Dataset(GDAL.gdalcreatecopy(
        driver.ptr,
        filename,
        dataset.ptr,
        strict,
        options,
        @cplprogress(progressfunc),
        progressdata
    ))
end

"""
    copy(dataset::AbstractDataset; [filename, [driver, [<keyword arguments>]]])

Create a copy of a dataset.

This method will attempt to create a copy of a raster dataset with the
indicated filename, and in this drivers format. Band number, size, type,
projection, geotransform and so forth are all to be copied from the
provided template dataset.

### Parameters
* `dataset`       the dataset being duplicated.

### Keyword Arguments
* `filename`      the filename for the new dataset. UTF-8 encoded.
* `driver`        the driver to use for creating the new dataset
* `strict`        ``true`` if the copy must be strictly equivalent, or more
    normally ``false`` if the copy may adapt as needed for the output format.
* `options`       additional format dependent options controlling creation
of the output file. `The APPEND_SUBDATASET=YES` option can be specified to
avoid prior destruction of existing dataset.

### Example
```
dataset = ArchGDAL.copy(originaldataset)
# work with dataset from here
```
or
```
ArchGDAL.copy(originaldataset) do dataset
    # work with dataset from here
end
```

### Returns
The newly created dataset.
"""
function copy(
        dataset::AbstractDataset;
        filename::AbstractString    = string("/vsimem/$(gensym())"),
        driver::Driver              = getdriver(dataset),
        strict::Bool                = false,
        options                     = StringList(C_NULL),
        progressfunc::Function      = GDAL.gdaldummyprogress,
        progressdata                = C_NULL
    )
    return IDataset(GDAL.gdalcreatecopy(
        driver.ptr,
        filename,
        dataset.ptr,
        strict,
        options,
        @cplprogress(progressfunc),
        progressdata
    ))
end

"""
    write(dataset::AbstractDataset, filename::AbstractString; kwargs...)

Writes the dataset to the designated filename.
"""
function write(dataset::AbstractDataset, filename::AbstractString; kwargs...)
    return destroy(unsafe_copy(dataset, filename = filename; kwargs...))
end

"""
    unsafe_create(filename::AbstractString; driver, width, height, nbands, dtype, options)

Create a new dataset.

What argument values are legal for particular drivers is driver specific, and
there is no way to query in advance to establish legal values.

That function will try to validate the creation option list passed to the driver
with the GDALValidateCreationOptions() method. This check can be disabled by
defining the configuration option GDAL_VALIDATE_CREATION_OPTIONS=NO.

After you have finished working with the returned dataset, it is required to
close it with GDALClose(). This does not only close the file handle, but also
ensures that all the data and metadata has been written to the dataset
(GDALFlushCache() is not sufficient for that purpose).

In some situations, the new dataset can be created in another process through
the GDAL API Proxy mechanism.

In GDAL 2, the arguments nXSize, nYSize and nBands can be passed to 0 when
creating a vector-only dataset for a compatible driver.
"""
function unsafe_create(
        filename::AbstractString;
        driver::Driver              = identifydriver(filename),
        width::Integer              = 0,
        height::Integer             = 0,
        nbands::Integer             = 0,
        dtype::DataType             = Any,
        options                     = StringList(C_NULL)
    )
    result = GDAL.gdalcreate(driver.ptr, filename, width, height, nbands,
        _GDALTYPE[dtype], options)
    return Dataset(result)
end

function unsafe_create(
        driver::Driver;
        filename::AbstractString    = string("/vsimem/$(gensym())"),
        width::Integer              = 0,
        height::Integer             = 0,
        nbands::Integer             = 0,
        dtype::DataType             = Any,
        options                     = StringList(C_NULL)
    )
    result = GDAL.gdalcreate(driver.ptr, filename, width, height, nbands,
        _GDALTYPE[dtype], options)
    return Dataset(result)
end

"""
    create(filename::AbstractString; driver, width, height, nbands, dtype, options)

Create a new dataset.

### Parameters
* `filename`       the filename for the dataset being created.

### Keyword Arguments
* `driver`        the driver to use for creating the new dataset
* `options`       additional format dependent options controlling creation
of the output file. `The APPEND_SUBDATASET=YES` option can be specified to
avoid prior destruction of existing dataset.
* `width`, `height`, `nbands`, `dtype`: only for raster datasets.

### Example
```
dataset = ArchGDAL.create(filename; ...)
# work with raster dataset from here
```
or
```
ArchGDAL.create(filename; ...) do dataset
    # work with vector dataset from here
end
```

### Returns
The newly created dataset.
"""
function create(
        filename::AbstractString;
        driver::Driver              = identifydriver(filename),
        width::Integer              = 0,
        height::Integer             = 0,
        nbands::Integer             = 0,
        dtype::DataType             = Any,
        options                     = StringList(C_NULL)
    )
    result = GDAL.gdalcreate(driver.ptr, filename, width, height, nbands,
        _GDALTYPE[dtype], options)
    return IDataset(result)
end

function create(
        driver::Driver;
        filename::AbstractString    = string("/vsimem/$(gensym())"),
        width::Integer              = 0,
        height::Integer             = 0,
        nbands::Integer             = 0,
        dtype::DataType             = Any,
        options                     = StringList(C_NULL)
    )
    result = GDAL.gdalcreate(driver.ptr, filename, width, height, nbands,
        _GDALTYPE[dtype], options)
    return IDataset(result)
end

"""
    unsafe_read(filename; flags=OF_ReadOnly, alloweddrivers, options, siblingfiles)

Open a raster file as a GDALDataset.

This function will try to open the passed file, or virtual dataset name by
invoking the Open method of each registered `GDALDriver` in turn. The first
successful open will result in a returned dataset. If all drivers fail then
`NULL` is returned and an error is issued.

### Parameters
* `filename`  the name of the file to access. In the case of exotic drivers
this may not refer to a physical file, but instead contain information for the
driver on how to access a dataset. It should be in UTF-8 encoding.
* `flags`     a combination of `GDAL_OF_*` flags (listed below) that may be
              combined through the logical `|` operator.

    - Driver kind: GDAL_OF_RASTER for raster drivers, GDAL_OF_VECTOR for vector
                   drivers. If none of the value is specified, both are implied.
    - Access mode: `OF_ReadOnly` (exclusive) or `OF_Update`.
    - Shared mode: `GDAL_OF_SHARED`. If set, it allows the sharing of
                   GDALDataset handles for a dataset with other callers that
                   have set GDAL_OF_SHARED. In particular, GDALOpenEx() will
                   consult its list of currently open and shared GDALDataset's,
                   and if the GetDescription() name for one exactly matches the
                   pszFilename passed to GDALOpenEx() it will be referenced and
                   returned, if GDALOpenEx() is called from the same thread.
    - Verbose error: GDAL_OF_VERBOSE_ERROR. If set, a failed attempt to open the
                   file will lead to an error message to be reported.
* `options`: additional format dependent options.

### Additional Remarks
Several recommendations:

* If you open a dataset object with `GA_Update` access, it is not recommended
to open a new dataset on the same underlying file.
* The returned dataset should only be accessed by one thread at a time. To use
it from different threads, you must add all necessary code (mutexes, etc.) to
avoid concurrent use of the object. (Some drivers, such as GeoTIFF, maintain
internal state variables that are updated each time a new block is read,
preventing concurrent use.)
* In order to reduce the need for searches through the operating system file
system machinery, it is possible to give an optional list of files with the
papszSiblingFiles parameter. This is the list of all files at the same level in
the file system as the target file, including the target file. The filenames
must not include any path components, are essentially just the output of
VSIReadDir() on the parent directory. If the target object does not have
filesystem semantics then the file list should be NULL.

In some situations (dealing with unverified data), the datasets can be opened
in another process through the GDAL API Proxy mechanism.

For drivers supporting the VSI virtual file API, it is possible to open a file
in a `.zip` archive (see `VSIInstallZipFileHandler()`), a `.tar/.tar.gz/.tgz`
archive (see `VSIInstallTarFileHandler()`), or a HTTP / FTP server
(see `VSIInstallCurlFileHandler()`)
"""
function unsafe_read(
        filename::AbstractString;
        flags           = OF_ReadOnly,
        alloweddrivers  = StringList(C_NULL),
        options         = StringList(C_NULL),
        siblingfiles    = StringList(C_NULL)
    )
    result = GDAL.gdalopenex(filename, Int(flags), alloweddrivers, options,
        siblingfiles)
    return Dataset(result)
end

"""
    read(filename; flags=OF_ReadOnly, alloweddrivers, options, siblingfiles)

Open a raster file

### Parameters
* `filename`: the filename of the dataset to be read.

### Keyword Arguments
* `flags`: a combination of `OF_*` flags (listed below) that may be
    combined through the logical `|` operator. It defaults to `OF_ReadOnly`.
    - Driver kind: `OF_Raster` for raster drivers, `OF_Vector` for vector
        drivers. If none of the value is specified, both are implied.
    - Access mode: `OF_ReadOnly` (exclusive) or `OF_Update`.
    - Shared mode: `OF_Shared`. If set, it allows the sharing of handles for a
        dataset with other callers that have set `OF_Shared`.
    - Verbose error: `OF_Verbose_Error`. If set, a failed attempt to open the
        file will lead to an error message to be reported.
* `options`: additional format dependent options.

### Example
```
dataset = ArchGDAL.read("point.shp")
# work with dataset from here
```
or
```
ArchGDAL.read("point.shp") do dataset
    # work with dataset from here
end
```

### Returns
The corresponding dataset.
"""
function read(
        filename::AbstractString;
        flags           = OF_ReadOnly | OF_Verbose_Error,
        alloweddrivers  = StringList(C_NULL),
        options         = StringList(C_NULL),
        siblingfiles    = StringList(C_NULL)
    )
    result = GDAL.gdalopenex(filename, Int(flags), alloweddrivers, options,
        siblingfiles)
    return IDataset(result)
end

unsafe_update(filename::AbstractString; flags = OF_Update, kwargs...) =
    unsafe_read(filename; flags = OF_Update | flags, kwargs...)

"""
    width(dataset::AbstractDataset)

Fetch raster width in pixels.
"""
width(dataset::AbstractDataset) = GDAL.gdalgetrasterxsize(dataset.ptr)

"""
    height(dataset::AbstractDataset)

Fetch raster height in pixels.
"""
height(dataset::AbstractDataset) = GDAL.gdalgetrasterysize(dataset.ptr)

"""
    nraster(dataset::AbstractDataset)

Fetch the number of raster bands on this dataset.
"""
nraster(dataset::AbstractDataset) = GDAL.gdalgetrastercount(dataset.ptr)

"""
    nlayer(dataset::AbstractDataset)

Fetch the number of feature layers on this dataset.
"""
nlayer(dataset::AbstractDataset) = GDAL.gdaldatasetgetlayercount(dataset.ptr)

"""
    getdriver(dataset::AbstractDataset)

Fetch the driver that the dataset was created with
"""
getdriver(dataset::AbstractDataset) =
    Driver(GDAL.gdalgetdatasetdriver(dataset.ptr))

"""
    filelist(dataset::AbstractDataset)

Fetch files forming dataset.

Returns a list of files believed to be part of this dataset. If it returns an
empty list of files it means there is believed to be no local file system files
associated with the dataset (for instance a virtual dataset). The returned file
list is owned by the caller and should be deallocated with `CSLDestroy()`.

The returned filenames will normally be relative or absolute paths depending on
the path used to originally open the dataset. The strings will be UTF-8 encoded
"""
filelist(dataset::AbstractDataset) = GDAL.gdalgetfilelist(dataset.ptr)

"""
    getlayer(dataset::AbstractDataset, i::Integer)

Fetch the layer at index `i` (between `0` and `nlayer(dataset)-1`)

The returned layer remains owned by the GDALDataset and should not be deleted by
the application.
"""
getlayer(dataset::AbstractDataset, i::Integer) =
    IFeatureLayer(GDAL.gdaldatasetgetlayer(dataset.ptr, i), ownedby = dataset)

unsafe_getlayer(dataset::AbstractDataset, i::Integer) =
    FeatureLayer(GDAL.gdaldatasetgetlayer(dataset.ptr, i))

"""
    getlayer(dataset::AbstractDataset, name::AbstractString)

Fetch the feature layer corresponding to the given name.

The returned layer remains owned by the GDALDataset and should not be deleted by
the application.
"""
function getlayer(dataset::AbstractDataset, name::AbstractString)
    return IFeatureLayer(
        GDAL.gdaldatasetgetlayerbyname(dataset.ptr, name),
        ownedby = dataset
    )
end

unsafe_getlayer(dataset::AbstractDataset, name::AbstractString) =
    FeatureLayer(GDAL.gdaldatasetgetlayerbyname(dataset.ptr, name))

"""
    deletelayer!(dataset::AbstractDataset, i::Integer)

Delete the indicated layer (at index i; between `0` to `nlayer()-1`)

### Returns
`OGRERR_NONE` on success, or `OGRERR_UNSUPPORTED_OPERATION` if deleting layers
is not supported for this dataset.
"""
function deletelayer!(dataset::AbstractDataset, i::Integer)
    result = GDAL.gdaldatasetdeletelayer(dataset.ptr, i)
    @ogrerr result "Failed to delete layer"
    return dataset
end

"""
    testcapability(dataset::AbstractDataset, capability::AbstractString)

Test if capability is available. `true` if capability available otherwise `false`.

One of the following dataset capability names can be passed into this function,
and a `true` or `false` value will be returned indicating whether or not the
capability is available for this object.

* `ODsCCreateLayer`: True if this datasource can create new layers.
* `ODsCDeleteLayer`: True if this datasource can delete existing layers.
* `ODsCCreateGeomFieldAfterCreateLayer`: True if the layers of this datasource
        support CreateGeomField() just after layer creation.
* `ODsCCurveGeometries`: True if this datasource supports curve geometries.
* `ODsCTransactions`: True if this datasource supports (efficient) transactions.
* `ODsCEmulatedTransactions`: True if this datasource supports transactions
        through emulation.

The #define macro forms of the capability names should be used in preference to
the strings themselves to avoid misspelling.

### Parameters
* `dataset`: the dataset handle.
* `capability`: the capability to test.
"""
testcapability(dataset::AbstractDataset, capability::AbstractString) =
    Bool(GDAL.gdaldatasettestcapability(dataset.ptr, capability))

function listcapability(
        dataset::AbstractDataset,
        capabilities = (GDAL.ODsCCreateLayer,
                        GDAL.ODsCDeleteLayer,
                        GDAL.ODsCCreateGeomFieldAfterCreateLayer,
                        GDAL.ODsCCurveGeometries,
                        GDAL.ODsCTransactions,
                        GDAL.ODsCEmulatedTransactions)
    )
    return Dict{String, Bool}(
        c => testcapability(dataset, c) for c in capabilities
    )
end

"""
    unsafe_executesql(dataset::AbstractDataset, query::AbstractString; dialect, spatialfilter)

Execute an SQL statement against the data store.

The result of an SQL query is either NULL for statements that are in error, or
that have no results set, or an OGRLayer pointer representing a results set from
the query. Note that this OGRLayer is in addition to the layers in the data
store and must be destroyed with ReleaseResultSet() before the dataset is closed
(destroyed).

For more information on the SQL dialect supported internally by OGR review the
OGR SQL document. Some drivers (i.e. Oracle and PostGIS) pass the SQL directly
through to the underlying RDBMS.

Starting with OGR 1.10, the SQLITE dialect can also be used.

### Parameters
* `dataset`: the dataset handle.
* `query`: the SQL statement to execute.
* `spatialfilter`: geometry which represents a spatial filter. Can be NULL.
* `dialect`: allows control of the statement dialect. If set to NULL, the
    OGR SQL engine will be used, except for RDBMS drivers that will use their
    dedicated SQL engine, unless OGRSQL is explicitly passed as the dialect.
    Starting with OGR 1.10, the SQLITE dialect can also be used.

### Returns
an OGRLayer containing the results of the query.
Deallocate with ReleaseResultSet().
"""
function unsafe_executesql(
        dataset::AbstractDataset,
        query::AbstractString;
        dialect::AbstractString = "",
        spatialfilter::Geometry = Geometry(GDALGeometry(C_NULL))
    )
    return FeatureLayer(GDALFeatureLayer(GDAL.gdaldatasetexecutesql(
        dataset.ptr,
        query,
        spatialfilter.ptr,
        dialect
    )))
end

"""
    releaseresultset(dataset::AbstractDataset, layer::FeatureLayer)

Release results of ExecuteSQL().

This function should only be used to deallocate OGRLayers resulting from an
ExecuteSQL() call on the same GDALDataset. Failure to deallocate a results set
before destroying the GDALDataset may cause errors.

### Parameters
* `dataset`: the dataset handle.
* `layer`: the result of a previous ExecuteSQL() call.
"""
function releaseresultset(dataset::AbstractDataset, layer::FeatureLayer)
    GDAL.gdaldatasetreleaseresultset(dataset.ptr, layer.ptr)
    destroy(layer)
end

"""
    getband(dataset::AbstractDataset, i::Integer)

Fetch a band object for a dataset from its index.
"""
getband(dataset::AbstractDataset, i::Integer) =
    IRasterBand(GDAL.gdalgetrasterband(dataset.ptr, i), ownedby = dataset)

unsafe_getband(dataset::AbstractDataset, i::Integer) =
    RasterBand(GDAL.gdalgetrasterband(dataset.ptr, i))

"""
    getgeotransform!(dataset::AbstractDataset, transform::Vector{Cdouble})

Fetch the affine transformation coefficients.

Fetches the coefficients for transforming between pixel/line (P,L) raster
space, and projection coordinates (Xp,Yp) space.

```julia
   Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
   Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];
```

In a north up image, `padfTransform[1]` is the pixel width, and
`padfTransform[5]` is the pixel height. The upper left corner of the upper left
pixel is at position `(padfTransform[0],padfTransform[3])`.

The default transform is `(0,1,0,0,0,1)` and should be returned even when a
`CE_Failure` error is returned, such as for formats that don't support
transformation to projection coordinates.

### Parameters
* `buffer`   a six double buffer into which the transformation will be placed.

### Returns
`CE_None` on success, or `CE_Failure` if no transform can be fetched.
"""
function getgeotransform!(dataset::AbstractDataset, transform::Vector{Cdouble})
    @assert length(transform) == 6
    result = GDAL.gdalgetgeotransform(dataset.ptr, pointer(transform))
    @cplerr result "Failed to get geotransform"
    return transform
end

getgeotransform(dataset::AbstractDataset) =
    getgeotransform!(dataset, Array{Cdouble}(undef, 6))

"""
    setgeotransform!(dataset::AbstractDataset, transform::Vector{Cdouble})

Set the affine transformation coefficients.
"""
function setgeotransform!(dataset::AbstractDataset, transform::Vector{Cdouble})
    @assert length(transform) == 6
    result = GDAL.gdalsetgeotransform(dataset.ptr, pointer(transform))
    @cplerr result "Failed to transform raster dataset"
    return dataset
end

"""
    ngcp(dataset::AbstractDataset)

Get number of GCPs for this dataset. Zero if there are none.
"""
ngcp(dataset::AbstractDataset) = GDAL.gdalgetgcpcount(dataset.ptr)

"""
    getproj(dataset::AbstractDataset)

Fetch the projection definition string for this dataset in OpenGIS WKT format.

It should be suitable for use with the OGRSpatialReference class. When a
projection definition is not available an empty (but not `NULL`) string is
returned.
"""
getproj(dataset::AbstractDataset) = GDAL.gdalgetprojectionref(dataset.ptr)

"""
    setproj!(dataset::AbstractDataset, projstring::AbstractString)

Set the projection reference string for this dataset.
"""
function setproj!(dataset::AbstractDataset, projstring::AbstractString)
    result = GDAL.gdalsetprojection(dataset.ptr, projstring)
    @cplerr result "Could not set projection"
    return dataset
end

"""
    buildoverviews!(dataset::AbstractDataset, overviewlist::Vector{Cint}; bandlist, resampling="NEAREST",
        progressfunc, progressdata)
        Build raster overview(s).

If the operation is unsupported for the indicated dataset, then CE_Failure is
returned, and CPLGetLastErrorNo() will return CPLE_NotSupported.

### Parameters
* `overviewlist` overview decimation factors to build.

### Keyword Parameters
* `panBandList`  list of band numbers. Must be in Cint (default = all)
* `sampling`     one of "NEAREST" (default), "GAUSS","CUBIC","AVERAGE","MODE",
                 "AVERAGE_MAGPHASE" or "NONE" controlling the downsampling
                 method applied.
* `progressfunc` a function to call to report progress, or `NULL`.
* `progressdata` application data to pass to the progress function.
"""
function buildoverviews!(
        dataset::AbstractDataset,
        overviewlist::Vector{Cint};
        bandlist::Vector{Cint}     = Cint[],
        resampling::AbstractString = "NEAREST",
        progressfunc::Function     = GDAL.gdaldummyprogress,
        progressdata               = C_NULL
    )
    result = GDAL.gdalbuildoverviews(
        dataset.ptr,
        resampling,
        length(overviewlist),
        overviewlist,
        length(bandlist),
        bandlist,
        @cplprogress(progressfunc),
        progressdata
    )
    @cplerr result "Failed to build overviews"
    return dataset
end

function destroy(dataset::AbstractDataset)
    GDAL.gdalclose(dataset.ptr)
    dataset.ptr = C_NULL
end


text:
using Statistics

export BVHAccelerator

"""
    PrimitiveInfo

A type holding information about a primitive (its index in an array, its centre point and
its bounding box). These are used when constructing the BVH accelerator tree structure.

"""
struct PrimitiveInfo
    number::Int64
    centroid::Point3
    BBox::BoundingBox
end

"""
    LinearBVHNode

A node of the flattened version of the BVH three. These nodes are used in the actual
ray-tracer acceleration.

# Fields

- `BBox::BoundingBox`: Bounding box of the primitive(s) contained in this node.
- `offset::Int64`: If this is a leaf node, contains the index of the corresponding primitive
  in the `primitives` array of the `BVHAccelerator`.
- `leaf::Bool`: Whether this is a leaf or branch node.
- `axis::Int64`: The axis along which the primitives are split if this is a branch node.
  Knowing this relative to the ray direction is relevant in the traversal of the tree.

"""
struct LinearBVHNode
    BBox::BoundingBox
    offset::Int64
    leaf::Bool
    axis::Int64
end

"""
    BVHBuildNode

A mutable type representing a node in the Bounding Volume Hierarchy tree. This is used when
creating the tree, which is then flattened into a more efficient structure using the
immutable `LinearBVHNode` types.

"""
mutable struct BVHBuildNode
    leaf::Bool
    split_axis::Int64
    first_offset::Int64
    BBox::BoundingBox
    childA::BVHBuildNode
    childB::BVHBuildNode
    BVHBuildNode() = new()
end

"""
    init_leaf!(node::BVHBuildNode, idx::Integer, BB::BoundingBox)

Initialize the `BVHBuildNode` as a leaf node of the BVH tree.

"""
function init_leaf!(node::BVHBuildNode, idx::Integer, BB::BoundingBox)
    node.leaf = true
    node.first_offset = idx
    node.BBox = BB
    node.split_axis = -1
    return node
end

"""
    init_interior!(node::BVHBuildNode, idx::Integer, BB::BoundingBox)

Initialize the `BVHBuildNode` as an interior node of the BVH tree.

"""

function init_interior!(node::BVHBuildNode, axis::Integer, childA::BVHBuildNode, childB::BVHBuildNode)
    node.leaf = false
    node.split_axis = axis
    node.BBox = BoundingBox(childA.BBox, childB.BBox)
    node.childA = childA
    node.childB = childB
    node.first_offset = -1
    return node
end

"""
    BVH_build_recursive(build_data, startN, endN, primitives, ordered_primitives,
                        total_nodes)

Recursively construct the BVH tree using `build_data` and `primitives`. This will populate
the array `ordered_primitives` with the primitives in the same order as to the tree leaf
nodes which refer to them.

"""
function BVH_build_recursive(build_data::Array{PrimitiveInfo,1}, startN::Integer,
                             endN::Integer, primitives::Array{T,1}, ordered_primitives::Array{T,1},
                             total_nodes::Integer) where T <: Primitive

    total_nodes += 1
    node = BVHBuildNode()

    if startN < 1 || endN < 1
        error()
    end

    if startN == endN
        # Only one primitive in the range, create leaf node
        n = build_data[startN].number
        push!(ordered_primitives, primitives[n])
        offset = length(ordered_primitives)
        init_leaf!(node, offset, build_data[startN].BBox)
    else
        # More than one primitive in the range, make partition and call recursively
        CBox = BoundingBox([build_data[i].centroid for i = startN:endN])
        dim = max_extent(CBox)
        pmid = mean(CBox.pMin, CBox.pMax)
        cmp(x,y) = x.centroid[dim] < y.centroid[dim]
        build_data[startN:endN] = sort(build_data[startN:endN] ; lt=cmp)
        mid = ceil(Int64, (startN+endN) / 2)

        # TODO: FIX THIS BETTER! THIS HAS SOME SUBTLE BUG
        #for i = startN:endN
        #    if (pmid[dim] <= build_data[i].centroid[dim])
        #        mid = i
        #        break
        #    end
        #end

        childA, total_nodes = BVH_build_recursive(build_data, startN, mid-1, primitives,
                                                  ordered_primitives, total_nodes)
        childB, total_nodes = BVH_build_recursive(build_data, mid, endN, primitives,
                                                  ordered_primitives, total_nodes)
        init_interior!(node, dim, childA, childB)
    end
    return node, total_nodes
end

"""
    flatten_BVH!(node::BVHBuildNode, linear::Array{LinearBVHNode,1}, idx::Integer)

Recursively traverse the BVH tree starting from `node`, and produce the flattened version
of the tree in the array `linear`.

"""
function flatten_BVH!(node::BVHBuildNode, linear::Array{LinearBVHNode,1}, idx::Integer)
    idx0 = idx
    if node.leaf
        new_node = LinearBVHNode(node.BBox, node.first_offset, true, -1)
    else
        idx1 = flatten_BVH!(node.childA, linear, idx+1)
        idx2 = flatten_BVH!(node.childB, linear, idx1+1)
        new_node = LinearBVHNode(node.BBox, idx1+1, false, node.split_axis)
        idx = idx2
    end
    linear[idx0] = new_node
    return idx
end

# LEFT: your number is my number, plus 1. Also, give me the biggest number in your subtree.
# RIGHT: your number is the biggest number in the LEFT subtree, plus 1
# I got to give upwards the biggest number in my subtree. But that's just the highest number
# in my right subtree.

"""
    BVHAccelerator{P<:Primitive} <: Aggregate

The Bounding Volume Hierarchy accelerator keeps the list of primitives, and a list of
`LinearBVHNode` objects, which contain all the information about the tree structure.

"""
struct BVHAccelerator{P<:GeometricPrimitive} <: Aggregate
    primitives::Array{P,1}
    nodes::Array{LinearBVHNode,1}
end

"""
    world_bounds(B::BVHAccelerator)

Return the bounding box containing all the primitives in the Accelerator.

"""
world_bounds(B::BVHAccelerator) = B.nodes[1].BBox

"""
    BVHAccelerator(prims::Array{T,1}) where T<:Primitive

Construct a `BVHAccelerator` from the list of primitives.

"""
function BVHAccelerator(prims::Array{T,1}) where T<:Primitive
    # Refine primitives to their fully intersectable components
    primtype = eltype(prims)
    primitives = primtype[]
    for p in prims
        fully_refine!(p, primitives)
    end

    # Extract relevant info from primitives and push into an array
    build_data = PrimitiveInfo[]
    for (i, obj) in enumerate(primitives)
        BB = world_bounds(obj)
        c = mean(BB.pMin, BB.pMax)
        info = PrimitiveInfo(i, c, BB)
        push!(build_data, info)
    end

    # Generate binary tree structure of bounding boxes
    ordered = primtype[]
    root_node, total_nodes = BVH_build_recursive(build_data, 1, length(primitives), primitives,
                                                 ordered, 0)

    # Transform the tree structure into a more efficient form
    linear = Array{LinearBVHNode}(undef,total_nodes)
    flatten_BVH!(root_node, linear, 1)
    BVHAccelerator(ordered, linear)
end

# Two arrays that are reused by the `intersect` and `intersectP` methods to
# avoid reallocating them every time those functions are called (which is very
# often). A significant optimization.
const TODO_ARRAY = zeros(Int64, (64, Threads.nthreads()))
const dir_is_neg = zeros(Bool, (3, Threads.nthreads()))

IntersectionTypes = Union{
                    Nothing,
                    Intersection{Lambert{SampledSpectrum}},
                    Intersection{SpecularDiffuse{SampledSpectrum}},
                    Intersection{AshikhminShirleySingle{SampledSpectrum}}
                    }

"""
    intersect(ray::Ray, BVH::BVHAccelerator)

Find the nearest intersection of the ray and the primitives in the `BVHAccelerator`. If you
only need to know whether a ray intersects something or not, use [`intersectP`](@ref)
instead, since `intersect` must find _all possible_ intersections to determine the nearest
one and is therefore much slower.

"""
function intersect(ray::Ray, BVH::BVHAccelerator)
    if length(BVH.nodes) == 0
        return nothing
    end
    thread_id = Threads.threadid()
    nodeN = 1
    todo_offset = 0
    todo = TODO_ARRAY
    for i = 1:64
        todo[i,thread_id] = 0
    end
    for i = 1:3
        dir_is_neg[i,thread_id] = ray.direction[i] < 0.0
    end

    tmin = Inf
    best_isect = nothing

    while true
        # Check the node given by the index nodeN
        node = BVH.nodes[nodeN]
        if intersectP(ray, node.BBox)
            # We hit the bounding box of this node.
            if node.leaf
                # This is a leaf node.
                # Check intersection with the primitive in the node and keep best intersection.
                isect = intersect(ray, BVH.primitives[node.offset])::IntersectionTypes
                best_isect = update_isect(isect, best_isect)::IntersectionTypes

                # If there's nothing in the todo queue, we exit the loop.
                if todo_offset == 0
                    break
                end
                # Othewise, set nodeN based on current todo value and shift todo queue left.
                nodeN = todo[todo_offset,thread_id]
                todo_offset -= 1
            else
                # This is not a leaf node.
                # We point nodeN to either the left or right child node, depending on the ray
                # direction. We shift right in the todo queue and put the other node there.
                todo_offset += 1
                if dir_is_neg[node.axis,thread_id]
                    # Put left child in todo queue and set nodeN to right child.
                    todo[todo_offset,thread_id] = nodeN + 1
                    nodeN = node.offset
                else
                    # Put right child in todo queue and set nodeN to left child.
                    todo[todo_offset,thread_id] = node.offset
                    nodeN = nodeN + 1
                end
            end
        else
            # We did not hit a bounding box.
            # If there's nothing in the todo queue, we exit.
            if todo_offset == 0
                break
            end
            # Otherwise, set nodeN based on current todo value and shift left in the todo queue.
            nodeN = todo[todo_offset,thread_id]
            todo_offset -= 1

        end
    end
    return best_isect
end

"""
    intersectP(ray::Ray, BVH::BVHAccelerator)

Returns `true` if the given ray intersects any primitive in the `BVHAccelerator`, else
returns false. This is much faster to do than [`intersect`](@ref), since it's enough to
find any intersection, not specifically the nearest one.

Possibly the most often called function in the entire ray-tracer.

"""
function intersectP(ray::Ray, BVH::BVHAccelerator)
    if length(BVH.nodes) == 0
        return false
    end
    thread_id = Threads.threadid()
    nodeN = 1
    todo_offset = 0
    todo = TODO_ARRAY
    for i = 1:64
        todo[i,thread_id] = 0
    end
    for i = 1:3
        dir_is_neg[i,thread_id] = ray.direction[i] < 0.0
    end
    while true
        node = BVH.nodes[nodeN]
        if intersectP(ray, node.BBox)
            if node.leaf
                if intersectP(ray, BVH.primitives[node.offset])
                    return true
                end
                if todo_offset == 0
                    break
                end
                nodeN = todo[todo_offset,thread_id]
                todo_offset -= 1
            else
                todo_offset += 1
                if dir_is_neg[node.axis,thread_id]
                    todo[todo_offset,thread_id] = nodeN + 1
                    nodeN = node.offset
                else
                    todo[todo_offset,thread_id] = node.offset
                    nodeN = nodeN + 1
                end
            end
        else
            if todo_offset == 0
                break
            end
            nodeN = todo[todo_offset,thread_id]
            todo_offset -= 1

        end
    end
    return false
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using FFTW, LinearAlgebra, Test, WaveFD

@testset "Wavelets" begin
    @testset "Minimum phase Ricker wavelet tests" begin
        dt = .001
        nt = 512
        wav = WaveletMinPhaseRicker()
        for t0 in [-51 -50 -49 0 49 50 51]
            w = get(wav, dt*t0 .+ dt * collect(0:nt-1))
            @test length(w) == nt
        end

        # check amplitude spectrum of min-phase ricker
        nt=1024
        dt=.004
        w1 = get(WaveletRicker(f=15.0), dt*collect(-div(nt,2):div(nt,2)-1))
        w2 = get(WaveletMinPhaseRicker(f=15.0), dt*collect(0:nt-1))
        W1 = abs.(rfft(w1))
        W2 = abs.(rfft(w2))
        err = norm(abs.(W2) - abs.(W1))/length(W1)
        @test isapprox(err, 0.0, atol=.01*maximum(W1))
    end

    @testset "Types, wtype=$(wtype)" for wtype in (WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletOrmsby, WaveletMinPhaseOrmsby)
        @test eltype(get(wtype(), Float64(.004)*collect(0:511))) == Float64
        @test eltype(get(wtype(), Float32(.004)*collect(0:511))) == Float32
    end

    @testset "copy methods" begin
        for W in (WaveletSine, WaveletRicker, WaveletMinPhaseRicker, WaveletDerivRicker, WaveletCausalRicker, WaveletOrmsby, WaveletMinPhaseOrmsby)
            global w = W()
            global wcopy = copy(w)
            for f in fieldnames(typeof(w))
                v = eval(Meta.parse("w.$(f)"))
                vcopy = eval(Meta.parse("wcopy.$(f)"))
                @test v ≈ vcopy
            end
        end
    end
end

nothing

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD