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

reads a 2-D grid file and fits a low-order polynomial trend to these data by
[optionally weighted] least-squares.

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

Parameters
----------

- **N** : **model** : -- Str or Number --

    Sets the number of model parameters to fit.
	[`-N`](http://gmt.soest.hawaii.edu/doc/latest/grdtrend.html#n)
- **D** : **diff** : -- Str or [] --

	Compute the difference (input data - trend). Optionaly provide a file name to save result on disk.
	[`-D`](http://gmt.soest.hawaii.edu/doc/latest/grdtrend.html#d)
- **T** : **trend** : -- Str or [] --

	Compute the trend surface. Optionaly provide a file name to save result on disk.
	[`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdtrend.html#t)
- **W** : **weights** : -- Str --

    If weight.nc exists, it will be read and used to solve a weighted least-squares problem.
	[`-W`](http://gmt.soest.hawaii.edu/doc/latest/grdtrend.html#w)
- $(GMT.opt_R)
- $(GMT.opt_V)
"""
function grdtrend(cmd0::String="", arg1=[], arg2=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("grdtrend", cmd0, arg1, arg2)	# Speedy mode

	if (isempty(cmd0) && isempty_(arg1))
		error("Must provide the grid to work with.")
	end

	d = KW(kwargs)

	cmd = add_opt("", 'N', d, [:N :model])
	if (!occursin("-N", cmd))
		error("The 'model' parameter is mandatory")
	end

	cmd, = parse_R(cmd, d)
	cmd = parse_V_params(cmd, d)
	cmd = add_opt(cmd, 'D', d, [:D :diff])
	cmd = add_opt(cmd, 'T', d, [:T :trend])

	if (occursin("-D", cmd) && occursin("-T", cmd))
		warning("Usage error, both difference and trend were required. Ignoring the trend request.")
	elseif (!occursin("-D", cmd) && !occursin("-T", cmd))
		cmd = cmd * " -T" 			# No -T -or -D provided so default to -T
	end

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end

	for sym in [:W :weights]
		if (haskey(d, sym))
			if (!isa(d[sym], GMTgrid))
				cmd = cmd * " -W" * arg2str(d[sym])
			else
				cmd, N_used = put_in_slot(cmd, d[sym], 'W', [arg1, arg2])
				if (N_used == 1)     arg1 = d[sym]
				elseif (N_used == 2) arg2 = d[sym]
				end
			end
			break
		end
	end

	if (isempty_(arg2))
		return common_grd(d, cmd, got_fname, 1, "grdtrend", arg1)		# Finish build cmd and run it
	else
		return common_grd(d, cmd, got_fname, 2, "grdtrend", arg1, arg2)
	end
end

# ---------------------------------------------------------------------------------------------------
grdtrend(arg1=[], arg2=[], cmd0::String=""; kw...) = grdtrend(cmd0, arg1, arg2; kw...)


text:
export MtlCompileOptions

const MTLCompileOptions = Ptr{MtCompileOptions}

mutable struct MtlCompileOptions
    handle::MTLCompileOptions
end

Base.convert(::Type{MTLCompileOptions}, opts::MtlCompileOptions) = opts.handle
Base.unsafe_convert(::Type{MTLCompileOptions}, opts::MtlCompileOptions) = convert(MTLCompileOptions, opts.handle)

Base.:(==)(a::MtlCompileOptions, b::MtlCompileOptions) = a.handle == b.handle
Base.hash(opts::MtlCompileOptions, h::UInt) = hash(opts.handle, h)

function MtlCompileOptions()
    handle = mtNewCompileOpts()
    obj = MtlCompileOptions(handle)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(opts::MtlCompileOptions)
    opts.handle !== C_NULL && mtCompileOptsRelease(opts)
end

Base.propertynames(o::MtlCompileOptions) = (:fastmath, :languageversion)
function Base.getproperty(o::MtlCompileOptions, f::Symbol)
    if f === :handle
        return getfield(o, :handle)
    elseif f === :fastmath
        return mtCompileOptsFastMath(o)
    elseif f === :languageversion
        return mtCompileOptsLanguageVersion(o)
    else
        error("CompileOptions does not have field $f")
    end
end

function Base.setproperty!(o::MtlCompileOptions, f::Symbol, val)
    if f === :fastmath
        return mtCompileOptsFastMathSet(o, val)
    elseif f === :languageversion
        return mtCompileOptsLanguageVersionSet(o, val)
    else
        error("CompileOptions does not have field $f")
    end
end

function Base.show(io::IO, l::MtlCompileOptions)
    print(io, "CompileOptions(…)")
end

function Base.show(io::IO, ::MIME"text/plain", l::MtlCompileOptions)
    println(io, "CompileOptions:")
    println(io, " fastmath : ", l.fastmath)
    print(io, " version  : ", l.languageversion)
end

ext: jl


text:
seed!(42)
L = 3
D = 10
N = 20
b = 5
x = [rand(D) for _ in 1:N]
y = rand(N)

@testset "AnalyticVI" begin
    i = AnalyticVI()
    @test AnalyticVI(ϵ = 0.0001f0) isa AnalyticVI{Float32,1}
    @test AnalyticSVI(10, ϵ = 0.0001f0) isa AnalyticVI{Float32,1}
    @test repr(i) == "Analytic Variational Inference"
    xview = view(x, 1:N)
    yview = view(y, 1:N)
    i = AGP.tuple_inference(i, L, D, N, N, xview, yview)
    @test AGP.getρ(i) == 1.0
    @test AGP.isStochastic(i) == false
    @test i isa AnalyticVI{Float64,L}

    i = AGP.tuple_inference(AnalyticSVI(b), L, D, N, b, xview, yview)
    @test i isa AnalyticVI{Float64,L}
    @test AGP.getρ(i) == N / b
    @test AGP.isStochastic(i) == true
end


text:
module ExactOptimalTransport

using Distances
using MathOptInterface
using Distributions
using PDMats
using QuadGK
using StatsBase: StatsBase

using LinearAlgebra
using SparseArrays

export emd, emd2
export ot_cost, ot_plan, wasserstein, squared2wasserstein

const MOI = MathOptInterface

include("distances/bures.jl")
include("utils.jl")
include("exact.jl")
include("wasserstein.jl")

end


text:
export OSProc, Context

"""
    Processor

An abstract type representing a processing device and associated memory, where
data can be stored and operated on. Subtypes should be immutable, and
instances should compare equal if they represent the same logical processing
device/memory. Subtype instances should be serializable between different
nodes. Subtype instances may contain a pointer to a "parent" `Processor` to
make it easy to transfer data to/from other types of `Processor` at runtime.
"""
abstract type Processor end

const PROCESSOR_CALLBACKS = []

"""
    execute!(proc::Processor, f, args...) -> Any

Executes the function `f` with arguments `args` on processor `proc`. This
function can be overloaded by `Processor` subtypes to allow executing function
calls differently than normal Julia.
"""
execute!

"""
    iscompatible(proc::Processor, opts, f, args...) -> Bool

Indicates whether `proc` can execute `f` over `args` given `opts`. `Processor`
subtypes should overload this function to return `true` if and only if it is
essentially guaranteed that `f(args...)` is supported. Additionally,
`iscompatible_func` and `iscompatible_arg` can be overriden to determine
compatibility of `f` and `args` individually. The default implementation
returns `false`.
"""
iscompatible(proc::Processor, opts, f, args...) =
    iscompatible_func(proc, opts, f) &&
    all(x->iscompatible_arg(proc, opts, x), args)
iscompatible_func(proc::Processor, opts, f) = false
iscompatible_arg(proc::Processor, opts, x) = false

"""
    default_enabled(proc::Processor) -> Bool

Returns whether processor `proc` is enabled by default (opt-out). `Processor` subtypes can override this function to make themselves opt-in (default returns `false`).
"""
default_enabled(proc::Processor) = false

"""
    get_processors(proc::Processor) -> Vector{T} where T<:Processor

Returns the full list of processors contained in `proc`, if any. `Processor`
subtypes should overload this function if they can contain sub-processors. The
default method will return a `Vector` containing `proc` itself.
"""
get_processors(proc::Processor) = Processor[proc]

"""
    get_parent(proc::Processor) -> Processor

Returns the parent processor for `proc`. The ultimate parent processor is an
`OSProc`. `Processor` subtypes should overload this to return their most
direct parent.
"""
get_parent

"""
    move(from_proc::Processor, to_proc::Processor, x)

Moves and/or converts `x` such that it's available and suitable for usage on
the `to_proc` processor. This function can be overloaded by `Processor`
subtypes to transport arguments and convert them to an appropriate form before
being used for exection. This is additionally called on thunk results when
moving back to `from_proc` before being serialized over the wire as needed.
The default implementation breaks a single `move` call down into a sequence of
`move` calls, and is not intended to be maximally efficient.
"""
function move(ctx, from_proc::Processor, to_proc::Processor, x)
    if from_proc == to_proc
        return x
    end
    @debug "Initiating generic move"
    # Move to remote OSProc
    @debug "(Remote) moving $parent_proc to $grandparent_proc"
    root = get_parent_osproc(from_proc)
    x, parent_proc = remotecall_fetch(move_to_osproc, root.pid, from_proc, x)

    # Move to local OSProc
    remote_proc = parent_proc
    local_proc = OSProc()
    @debug "(Network) moving $remote_proc to $local_proc"
    x = move(ctx, remote_proc, local_proc, x)

    # Move to to_proc
    parent_proc = get_parent(to_proc)
    path = Processor[to_proc, parent_proc]
    while parent_proc != local_proc && !(parent_proc isa OSProc)
        parent_proc = get_parent(parent_proc)
        push!(path, parent_proc)
    end
    last_proc = local_proc
    while !isempty(path)
        next_proc = pop!(path)
        @debug "(Local) moving $last_proc to $next_proc"
        x = move(ctx, last_proc, next_proc, x)
        last_proc = next_proc
    end
    return x
end
function get_parent_osproc(proc)
    while !(proc isa OSProc)
        proc = get_parent(proc)
    end
    proc
end
function move_to_osproc(parent_proc, x)
    ctx = Context()
    while !(parent_proc isa OSProc)
        grandparent_proc = get_parent(parent_proc)
        x = move(ctx, parent_proc, grandparent_proc, x)
        parent_proc = grandparent_proc
    end
    return x, parent_proc
end

"""
    OSProc <: Processor

Julia CPU (OS) process, identified by Distributed pid. Executes thunks when
threads or other "children" processors are not available, and/or the user has
not opted to use those processors.
"""
struct OSProc <: Processor
    pid::Int
    attrs::Dict{Symbol,Any}
    children::Vector{Processor}
    queue::Vector{Processor}
end
const OSPROC_CACHE = Dict{Int,OSProc}()
OSProc(pid::Int=myid()) = get!(OSPROC_CACHE, pid) do
    remotecall_fetch(get_osproc, pid, pid)
end
function get_osproc(pid::Int)
    proc = OSProc(pid, Dict{Symbol,Any}(), Processor[], Processor[])
    for cb in PROCESSOR_CALLBACKS
        try
            child = Base.invokelatest(cb, proc)
            child !== nothing && push!(proc.children, child)
        catch err
            @error "Error in processor callback" exception=(err,catch_backtrace())
        end
    end
    proc
end
function add_callback!(func)
    push!(Dagger.PROCESSOR_CALLBACKS, func)
    empty!(OSPROC_CACHE)
end
Base.:(==)(proc1::OSProc, proc2::OSProc) = proc1.pid == proc2.pid
iscompatible(proc::OSProc, opts, f, args...) =
    any(child->iscompatible(child, opts, f, args...), proc.children)
iscompatible_func(proc::OSProc, opts, f) =
    any(child->iscompatible_func(child, opts, f), proc.children)
iscompatible_arg(proc::OSProc, opts, args...) =
    any(child->
        all(arg->iscompatible_arg(child, opts, arg), args),
    proc.children)
get_processors(proc::OSProc) =
    vcat((get_processors(child) for child in proc.children)...,)
function choose_processor(from_proc::OSProc, options, f, args)
    if isempty(from_proc.queue)
        for child in from_proc.children
            grandchildren = get_processors(child)
            append!(from_proc.queue, grandchildren)
        end
    end
    @assert !isempty(from_proc.queue)
    for i in 1:length(from_proc.queue)
        proc = popfirst!(from_proc.queue)
        push!(from_proc.queue, proc)
        if !iscompatible(proc, options, f, args...)
            continue
        end
        if default_enabled(proc) && isempty(options.proctypes)
            return proc
        elseif any(p->proc isa p, options.proctypes)
            return proc
        end
    end
    throw(ProcessorSelectionException(options.proctypes, from_proc.queue, f, args))
end
struct ProcessorSelectionException <: Exception
    proctypes::Vector{Type}
    procsavail::Vector{Processor}
    f
    args
end
function Base.show(io::IO, pex::ProcessorSelectionException)
    println(io, "(Worker $(myid())) Exhausted all available processor types!")
    println(io, "  Proctypes: $(pex.proctypes)")
    println(io, "  Procs Available: $(pex.procsavail)")
    println(io, "  Function: $(pex.f)")
    print(io, "  Arguments: $(pex.args)")
end

move(ctx, from_proc::OSProc, to_proc::OSProc, x) = x
execute!(proc::OSProc, f, args...) = f(args...)
default_enabled(proc::OSProc) = true

"""
    ThreadProc <: Processor

Julia CPU (OS) thread, identified by Julia thread ID.
"""
struct ThreadProc <: Processor
    owner::Int
    tid::Int
end
iscompatible(proc::ThreadProc, opts, f, args...) = true
iscompatible_func(proc::ThreadProc, opts, f) = true
iscompatible_arg(proc::ThreadProc, opts, x) = true
move(ctx, from_proc::OSProc, to_proc::ThreadProc, x) = x
move(ctx, from_proc::ThreadProc, to_proc::OSProc, x) = x
@static if VERSION >= v"1.3.0-DEV.573"
    execute!(proc::ThreadProc, f, args...) = fetch(Threads.@spawn f(args...))
else
    # TODO: Use Threads.@threads?
    execute!(proc::ThreadProc, f, args...) = fetch(@async f(args...))
end
get_parent(proc::ThreadProc) = OSProc(proc.owner)
default_enabled(proc::ThreadProc) = true

# TODO: ThreadGroupProc?

"A context represents a set of processors to use for an operation."
mutable struct Context
    procs::Vector{Processor}
    log_sink::Any
    profile::Bool
    options
end

"""
    Context(;nthreads=Threads.nthreads()) -> Context
    Context(xs::Vector{OSProc}) -> Context
    Context(xs::Vector{Int}) -> Context

Create a Context, by default adding each available worker once from
every available thread. Use the `nthreads` keyword to use a different
number of threads.

It is also possible to create a Context from a vector of [`OSProc`](@ref),
or equivalently the underlying process ids can also be passed directly
as a `Vector{Int}`.
"""
function Context(xs)
    Context(xs, NoOpLog(), false, nothing) # By default don't log events
end
Context(xs::Vector{Int}) = Context(map(OSProc, xs))
function Context()
    procs = [OSProc(w) for w in workers()]
    Context(procs)
end
procs(ctx::Context) = ctx.procs

"""
    write_event(ctx::Context, event::Event)

Write a log event
"""
function write_event(ctx::Context, event::Event)
    write_event(ctx.log_sink, event)
end


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

pols = Dict("neg" =>        :neg,
            "negative" =>   :neg,
            "-" =>          :neg,
            "pos" =>        :pos,
            "positive" =>   :pos,
            "+" =>          :pos)

include("Box.jl")
include("RectangularCuboid.jl")
include("Tube.jl")
include("Cone.jl")
include("Sphere.jl")
include("Mesh_primitive.jl")
include("plot_recipes.jl")


text:
immutable AccelSteeringAngle
    a::Float64 # accel [m/s²]
    δ::Float64 # steering angle [rad]
end
Base.show(io::IO, a::AccelSteeringAngle) = @printf(io, "AccelSteeringAngle(%6.3f,%6.3f)", a.a, a.δ)
Base.length(::Type{AccelSteeringAngle}) = 2
Base.convert(::Type{AccelSteeringAngle}, v::Vector{Float64}) = AccelSteeringAngle(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::AccelSteeringAngle)
    v[1] = a.a
    v[2] = a.δ
    v
end
function propagate(veh::Entity{VehicleState, BicycleModel, Int}, action::AccelSteeringAngle, roadway::Roadway, Δt::Float64)

    L = veh.def.a + veh.def.b
    l = -veh.def.b

    a = action.a # accel [m/s²]
    δ = action.δ # steering wheel angle [rad]

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    s = v*Δt + a*Δt*Δt/2 # distance covered
    v′ = v + a*Δt

    if abs(δ) < 0.01 # just drive straight
        posG = veh.state.posG + polar(s, θ)
    else # drive in circle

        R = L/tan(δ) # turn radius

        β = s/R
        xc = x - R*sin(θ) + l*cos(θ)
        yc = y + R*cos(θ) + l*sin(θ)

        θ′ = mod(θ+β, 2π)
        x′ = xc + R*sin(θ+β) - l*cos(θ′)
        y′ = yc - R*cos(θ+β) - l*sin(θ′)

        posG = VecSE2(x′, y′, θ′)
    end

    VehicleState(posG, roadway, v′)
end


text:
@testset "Testing Custom SF and Meas submodule" begin

customSF = joinpath(datadir,"BrukerStore","20200603_172811_EasyAxisContrast_1_1","245")
customMeas = joinpath(datadir,"BrukerStore","20200603_172811_EasyAxisContrast_1_1","247")

fnCustomSF_FF = joinpath(datadir,"mdf","customSystemMatrixFF.mdf")

fSF = MPIFile(customSF)
fMeas = MPIFile(customMeas)

data = MPIFiles.loadAndProcessFFData(fSF, 100, 4)
@test size(data) == (9,817,3,1)

convertCustomSF(fnCustomSF_FF, fSF, fSF, 100, 4)
dataSet = loadDataset(MPIFile(fnCustomSF_FF))
@test size(dataSet[:measData]) == (18,817,3,1)

timeData, pos = getFFdataPerPos(fMeas, 100, 4)
@test size(timeData) == (3,1632,3,1,100)
@test size(pos) == (3,1,3)

fileNames, paramsArray = prepareAsMDFSingleMeasurement(fMeas, 100, 4)
@test length(fileNames) == 3
@test length(paramsArray) == 3
end


text:
using LinearAlgebra

"""
    Far Fiel Condition: r_emitter^2 / 2r_observation ≪ 1

Condition comes from is the Eq (4.4) in
    https://engineering.purdue.edu/wcchew/ece604f18/latex%20lecture%20notes/LectureNotes20.pdf
(A copy of the file is inside the `benchmarks/pdf` folder)

In practice, we verify if the condition is smaller than `1e-1`,
because we tested that the estimate error is around 1%.
"""
function test_farField_validity(r_observation, r_emitter)
    if abs( norm(r_emitter)^2/(2*norm(r_observation))) ≤ 1e-1
        return true
    else
        return false
    end
    nothing
end

r_atom = 10*rand(3)
scalings = 10.0.^range(1, 7, length=100)

all_error_real = []
all_error_imag = []
conditions = []
for rep=1:200
    for scaling_sensor = scalings
        r_sensor = scaling_sensor.*rand()*ones(3)

        vector_distance_SensorAtom = r_sensor - r_atom
        distance_SensorAtom = norm(vector_distance_SensorAtom)

        distance_sensor = norm(r_sensor)
        direction_sensor = r_sensor / distance_sensor

        E_exact = -0.5im * cis(distance_SensorAtom) / (im * distance_SensorAtom)
        E_farfield = -0.5im *cis(distance_sensor) * cis(-dot(direction_sensor, r_atom)) / (im * distance_sensor)

        diffE = E_exact - E_farfield
        error_real = 100*abs(real(diffE)/real(E_exact))
        error_imag = 100*abs(imag(diffE)/imag(E_exact))

        push!(all_error_real, error_real)
        push!(all_error_imag, error_imag)

        condition = abs( norm(r_atom)^2/(2*norm(r_sensor)))
        push!(conditions, condition)
    end
end

using Plots

scatter(conditions, all_error_real, label="error real",markershape=:circle,markersize=6, xflip=true)
scatter!(conditions, all_error_imag, label="error imag", markershape=:rect, markersize=5,scale=:log10)
xlabel!("Validity Condition")
ylabel!("Percentage Error")
hline!([1], c=:black, label="1%", linestyle=:dash, lw=2)
vline!([1e-1], c=:red, label="Aceptable Valid Condition", lw=2)

#=
    I state that an error aceptable is around 1% - the black dashed line.
    Therefore, the `Validity Condition` required is around 10^-1.
    Which is exact reason of the `10^-1` used in test_farField_validity() function.

    Now that I have a defined expression, I can figure out the Far Field distance:

    |r_emitter|^2/(2*|r_observation|) ≤ 1e-1
    |r_emitter|^2/(2*1e-1) ≤ |r_observation|
    |r_emitter|^2/(0.2) ≤ |r_observation|
    |r_emitter|^2/(1/5) ≤ |r_observation|

    |r_observation| ≥ 5|r_emitter|^2

    Just to be a little bit above this condition, I will define
    a function with "5.01" and not "5.0"
=#
function estimate_far_Field(r_emitter)
    return 5.01*(norm(r_emitter)^2)
end


text:
## Julia 1.1.0

λ = map

λ(x -> print("$x "),
    λ(x -> x^2 + 2x - 1, [1, 2, 3, 4]))

println()

λ(x -> print("$x "),
    λ(x -> begin
        if x < 0 && iseven(x)
            return 0
        elseif x == 0
            return 1
        else
            return x
        end
    end, [0, 1, 2, 3]))

println()

map(x -> print("$x "),
    map([0, 1, 2, 3]) do x
        if x < 0 && iseven(x)
            return 0
        elseif x == 0
            return 1
        else
            return x
        end
    end)

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
in_size = 28 * 28

std_transform = true
pca_transform = true

##

using Knet: dir ; include(dir("data","mnist.jl"))

m_nist = mnist()

samples1 = []
for i in 1:size(m_nist[1])[end]
    sample = reshape(m_nist[1][:,:,:,i], 1, in_size)
    push!(samples1, sample)
end ;
samples1_ = vcat(samples1...)

samples2 = []
for i in 1:size(m_nist[3])[end]
    sample = reshape(m_nist[3][:,:,:,i], 1, in_size)
    push!(samples2, sample)
end ;
samples2_ = vcat(samples2...)

nist = [
        reshape(samples1_, size(m_nist[1])[end], in_size),
        m_nist[2],
        reshape(samples2_, size(m_nist[3])[end], in_size),
        m_nist[4],
       ]

import StatsBase
# using StatsBase: fit, transform, reconstruct
# using StatsBase: ZScoreTransform, UnitRangeTransform

import MultivariateStats
# using MultivariateStats: fit, transform, reconstruct
# using MultivariateStats: PCA, outdim

std_transform ? (begin

    training_and_test = vcat(nist[1],nist[3])

    pca_transform ? (begin

        transform2 = MultivariateStats.fit(MultivariateStats.PCA, training_and_test', maxoutdim=in_size, pratio=1.0)
        training_and_test = MultivariateStats.transform(transform2,training_and_test')'

        in_size = MultivariateStats.outdim(transform2)

    end) : (begin

        transform1 = StatsBase.fit(StatsBase.ZScoreTransform, training_and_test,dims=1)
        training_and_test = StatsBase.transform(transform1,training_and_test)

        for e in training_and_test # TODO : find out why and remove me.
            if isnan(e)
                println("nan detected. transform1")
            end
        end

    end)

    transform3 = StatsBase.fit(StatsBase.UnitRangeTransform, training_and_test,dims=1)
    traning_and_test = StatsBase.transform(transform3,training_and_test)

    nist[1] = training_and_test[1:size(nist[1])[1],:]
    nist[3] = training_and_test[size(nist[1])[1]+1:end,:]

end) : ()

to_label(int_val) =
begin
    lbl = zeros(1, 10)
    lbl[int_val] += 1
lbl
end

fn_01_to_minus1plus1 = x->(x.*2).-1

data_train = []
data_dev  = []

data_train2 = [[] for _ in 1:10]
data_dev2 = [[] for _ in 1:10]

for i in 1:size(nist[1])[1]
    sample = reshape(nist[1][i,:], 1, in_size)
    sample = fn_01_to_minus1plus1(sample)
    binary ? sample = binarize_data(sample) : ()
    push!(data_train, sample)
    push!(data_train2[nist[2][i]], sample)
end ;

for i in 1:size(nist[3])[1]
    sample = reshape(nist[3][i,:], 1, in_size)
    sample = fn_01_to_minus1plus1(sample)
    binary ? sample = binarize_data(sample) : ()
    push!(data_dev, sample)
    push!(data_dev2[nist[4][i]], sample)
end ;

##

# function scale_features(X)
#     μ = mean(X, dims=1)
#     σ = std(X, dims=1)
#
#     X_norm = (X .- μ) ./ σ
#
#     return (X_norm, μ, σ)
# end
#
#
# # Scale the testing features using the learned parameters
# function transform_features(X, μ, σ)
#     X_norm = (X .- μ) ./ σ
#     return X_norm
# end


text:
@enum DepwarnFlag DepwarnOff=0 DepwarnOn=1 DepwarnError=2

doc"""
- `switch_depwarn!(flag :: Bool)`
- `switch_depwarn!(flag :: DepwarnFlag)`
Enable/Disable deprecation warning.
- `DepwarnOff` or `false` : switch off deprecation warning
- `DepwarnOn` or `true` : switch on deprecation warning
- `DepwarnError` : turn deprecation warning into error
"""
switch_depwarn!(flag :: Bool) = switch_depwarn!(flag ? DepwarnOn : DepwarnOff)
function switch_depwarn!(flag :: DepwarnFlag)
    old_opt = Base.JLOptions()
    params = map(fieldnames(Base.JLOptions)) do name
        name == :depwarn ? Int(flag) : getfield(old_opt, name)
    end
    new_opt = Base.JLOptions(params...)
    unsafe_store!(cglobal(:jl_options, Base.JLOptions), new_opt)
    flag
end

# one-liner
# unsafe_store!(cglobal(:jl_options, Base.JLOptions), Base.JLOptions(map(fieldnames(Base.JLOptions)) do name; name==:depwarn?0:getfield(Base.JLOptions(), name) end...))


text:
module TestCLIWrapper

using TSML
using Statistics
using Test

using TSML.CLIWrappers:tsmlrun,rawstat,aggregatedstat,aggregatedoutput,imputedstat,imputedoutput

const inputname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
const outputname = joinpath(dirname(pathof(TSML)),"../data/testdata_output.csv")

function test_cli()
  op1=tsmlrun(inputname,"")
  op2=imputedoutput(inputname,"")
  @test mean(op1[:Value] .== op2[:Value]) == 1.0
  op5=tsmlrun(inputname,"","dd/mm/yyyy HH:MM","stat")
  op6=imputedstat(inputname,"")
  @test sum(op5[1,:] .== op6[1,:]) == 20
end
@testset "CLIWrapper: check cliwrapper" begin
  test_cli()
end

end


text:
# Julia markdown to HTML for HW 3
using Weave
weave("saulnier-hw3.jmd"; doctype = "md2html", out_path = :pwd)


text:
module Letters

using Random

export randstring

const MY_NAME = "Letters"

function rand()
  Random.rand('A':'Z')
end

function randstring()
  [rand() for _ in 1:10] |> join
end

include("module_name.jl")

# include("Numbers.jl")

end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
using Libxc
include("../xc/xc_evaluate.jl")

"""
Exchange-correlation term, defined by a list of functionals and usually evaluated through libxc.
"""
struct Xc
    functionals::Vector{Functional}
    scaling_factor::Real  # to scale by an arbitrary factor (useful for exploration)
end
Xc(functionals::Vector; scaling_factor=1) = Xc(functionals, scaling_factor)
Xc(functionals::Vector{Symbol}; kwargs...) = Xc(Functional.(functionals); kwargs...)
Xc(functional::Symbol; kwargs...) = Xc([Functional(functional)]; kwargs...)
Xc(functionals::Symbol...; kwargs...) = Xc([functionals...]; kwargs...)
(xc::Xc)(basis) = TermXc(basis, xc.functionals, xc.scaling_factor)

struct TermXc <: Term
    basis::PlaneWaveBasis
    functionals::Vector{Functional}
    scaling_factor::Real
end

function ene_ops(term::TermXc, ψ, occ; ρ, kwargs...)
    basis = term.basis
    T = eltype(basis)
    model = basis.model
    @assert all(xc.family in (:lda, :gga) for xc in term.functionals)

    if isempty(term.functionals)
        ops = [NoopOperator(term.basis, kpoint) for kpoint in term.basis.kpoints]
        return (E=0, ops=ops)
    end

    # Take derivatives of the density if needed.
    max_ρ_derivs = maximum(max_required_derivative, term.functionals)
    density = DensityDerivatives(basis, max_ρ_derivs, ρ)

    potential = zeros(T, basis.fft_size)
    zk = zeros(T, basis.fft_size)  # Energy per unit particle
    E = zero(T)
    for xc in term.functionals
        # Evaluate the functional and its first derivative (potential)
        terms = evaluate(xc; input_kwargs(xc.family, density)..., zk=zk)

        # Add energy contribution
        dVol = basis.model.unit_cell_volume / prod(basis.fft_size)
        E += sum(terms.zk .* ρ.real) * dVol

        # Add potential contributions Vρ -2 ∇⋅(Vσ ∇ρ)
        potential .+= terms.vrho
        if haskey(terms, :vsigma)
            potential .+= -2 * divergence_real(α -> terms.vsigma .* density.∇ρ_real[α],
                                               density.basis)
        end
    end
    if term.scaling_factor != 1
        E *= term.scaling_factor
        potential .*= term.scaling_factor
    end

    ops = [RealSpaceMultiplication(basis, kpoint, potential) for kpoint in basis.kpoints]
    (E=E, ops=ops)
end

function compute_kernel(term::TermXc; ρ::RealFourierArray, kwargs...)
    kernel = similar(ρ.real)
    kernel .= 0
    for xc in term.functionals
        xc.family == :lda || error("compute_kernel only implemented for LDA")
        terms = evaluate(xc; rho=ρ.real, derivatives=2:2)  # only valid for LDA
        kernel .+= terms.v2rho2
    end
    Diagonal(vec(term.scaling_factor .* kernel))
end

function apply_kernel(term::TermXc, dρ; ρ::RealFourierArray, kwargs...)
    basis = term.basis
    T = eltype(basis)
    @assert all(xc.family in (:lda, :gga) for xc in term.functionals)

    # Take derivatives of the density and the perturbation if needed.
    max_ρ_derivs = maximum(max_required_derivative, term.functionals)
    density = DensityDerivatives(basis, max_ρ_derivs, ρ)
    perturbation = DensityDerivatives(basis, max_ρ_derivs, from_real(basis, dρ))

    result = similar(ρ.real)
    result .= 0
    for xc in term.functionals
        # TODO LDA actually only needs the 2nd derivatives for this ... could be optimised
        terms = evaluate(xc; input_kwargs(xc.family, density)..., derivatives=1:2)

        # Accumulate LDA and GGA terms in result
        result .+= terms.v2rho2 .* dρ
        if haskey(terms, :v2rhosigma)
            result .+= apply_kernel_term_gga(terms, density, perturbation)
        end
    end
    term.scaling_factor .* result
end

#=
The total energy is
Etot = ∫ ρ E(ρ,σ), where σ = |∇ρ|^2 and E(ρ,σ) is the energy per unit particle
libxc provides the scalars
Vρ = ∂(ρ E)/∂ρ
Vσ = ∂(ρ E)/∂σ

Consider a variation dϕi of an orbital (considered real for
simplicity), and let dEtot be the corresponding variation of the
energy. Then the potential Vxc is defined by
dEtot = 2 ∫ Vxc ϕi dϕi

dρ = 2 ϕi dϕi
dσ = 2 ∇ρ ⋅ ∇dρ = 4 ∇ρ ⋅ ∇(ϕi dϕi)
dEtot = ∫ Vρ dρ + Vσ dσ
      = 2 ∫ Vρ ϕi dϕi + 4 ∫ Vσ ∇ρ ⋅ ∇(ϕi dϕi)
      = 2 ∫ Vρ ϕi dϕi - 4 ∫ div(Vσ ∇ρ) ϕi dϕi
where we performed an integration by parts in the last equation
(boundary terms drop by periodicity).

Therefore,
Vxc = Vρ - 2 div(Vσ ∇ρ)

See eg Richard Martin, Electronic stucture, p. 158

A more algebraic way of deriving these for GGA, which helps see why
this is correct at the discrete level, is the following. Let

B the (Fourier, spherical) space of orbitals
C the (Fourier, cube) space of Fourier densities/potentials
C* the (real-space, cube) space of real-space densities/potentials

X : B -> C be the extension operator
R : C -> B the restriction (adjoint of X)
IF : C -> C* be the IFFT
F : C* -> be the FFT, adjoint of IF
K : C -> C be the gradient operator (multiplication by iK, assuming dimension 1 to
simplify notations)

Then

ρ = IF(Xϕ).^2
ρf = F ρ
∇ρf = K ρf
∇ρ = IF ∇ρf

dρ = 2 IF(Xϕ) .* IF(X dϕ)
d∇ρ = IF K F dρ
    = 2 IF K F (IF(Xϕ) .* IF(X dϕ))

Etot = sum(ρ . * E(ρ,∇ρ.^2))
dEtot = sum(Vρ .* dρ + 2 (Vσ .* ∇ρ) .* d∇ρ)
      = sum(Vρ .* dρ + (IF K F) (Vσ .* ∇ρ) .* dρ)
(IF K F is self-adjoint; this is the discrete integration by parts)

Now note that

sum(V .* dρ) = sum(V .* IF(Xϕ) .* IF(X dϕ))
             = sum(R(F(V .* IF(Xϕ))) .* dϕ)
and the result follows.
=#

function max_required_derivative(functional)
    functional.family == :lda && return 0
    functional.family == :gga && return 1
    error("Functional family $(functional.family) not known.")
end

struct DensityDerivatives
    basis
    max_derivative::Int
    ρ         # density on real-space grid
    ∇ρ_real   # density gradient on real-space grid
    σ_real    # contracted density gradient on real-space grid
end

"""
DOCME compute density in real space and its derivatives starting from ρ
"""
function DensityDerivatives(basis, max_derivative::Integer, ρ::RealFourierArray)
    model = basis.model
    @assert model.spin_polarization == :none "Only spin_polarization == :none implemented."
    function ifft(x)
        tmp = G_to_r(basis, x)
        check_real(tmp)
        real(tmp)
    end

    ρF = ρ.fourier
    σ_real = nothing
    ∇ρ_real = nothing
    if max_derivative < 0 || max_derivative > 1
        error("max_derivative not in [0, 1]")
    elseif max_derivative > 0
        ∇ρ_real = [ifft(im * [(model.recip_lattice * G)[α] for G in G_vectors(basis)] .* ρF)
                   for α in 1:3]
        # TODO The above assumes CPU arrays
        σ_real = sum(∇ρ_real[α] .* ∇ρ_real[α] for α in 1:3)
    end

    DensityDerivatives(basis, max_derivative, ρ.real, ∇ρ_real, σ_real)
end

function input_kwargs(family, density)
    family == :lda && return (rho=density.ρ, )
    family == :gga && return (rho=density.ρ, sigma=density.σ_real)
    return NamedTuple()
end

"""
Compute divergence of an operand function, which returns the cartesian x,y,z
components in real space when called with the arguments 1 to 3.
The divergence is also returned as a real-space array.
"""
function divergence_real(operand, basis)
    gradsum = sum(1:3) do α
        operand_α = r_to_G(basis, complex(operand(α)))
        del_α = im * [(basis.model.recip_lattice * G)[α] for G in G_vectors(basis)]
        del_α .* operand_α
    end
    real(G_to_r(basis, gradsum))
end

function apply_kernel_term_gga(terms, density, perturbation)
    error("apply_kernel_term_gga is not yet working")
    # Follows DOI 10.1103/PhysRevLett.107.216402
    #
    # For GGA V = Vρ - 2 ∇⋅(Vσ ∇ρ) = (∂ε/∂ρ) - 2 ∇⋅((∂ε/∂σ) ∇ρ)
    #
    # dV(r) = f(r,r') dρ(r') = (∂V/∂ρ) dρ + (∂V/∂σ) dσ
    #
    # therefore
    # dV(r) = (∂^2ε/∂ρ^2) dρ - 2 ∇⋅((∂^2ε/∂σ∂ρ) ∇ρ) dρ
    #       + (∂^2ε/∂ρ∂σ) dσ - 2 ∇⋅((∂^ε/∂σ^2) ∇ρ + (∂ε/∂σ) (∂∇ρ/∂σ)) dσ
    #
    # Note dσ = 2∇ρ ∇dρ = 2∇ρ ∇dρ, therefore
    #    - 2 ∇⋅((∂ε/∂σ) (∂∇ρ/∂σ)) dσ = - 2 ∇⋅((∂ε/∂σ) ∇dρ)
    #
    # Note that below the LDA term (∂^2ε/∂ρ^2) dρ is ignored (already dealt with)
    ρ   = density.ρ
    ∇ρ  = density.∇ρ_real
    dρ  = perturbation.ρ
    ∇dρ = perturbation.∇ρ_real
    dσ = 2sum(∇ρ[α] .* ∇dρ[α] for α in 1:3)

    (
        terms.v2rhosigma .* dσ + divergence_real(density.basis) do α
            @. begin
                -2 * terms.v2rhosigma * ∇ρ[α] * dρ
                -2 * terms.v2sigma2 * ∇ρ[α] * dσ
                -2 * terms.vsigma * ∇dρ[α]
            end
        end
    )
end


text:
reload("EMST.jl")
using EMST

function f_benchmark( d , n , reps ; nmin=64)
    ti = Vector{Float64}()
    for zi=1:reps
        tic()
        x_sol = EMST.compute_emst(rand(d,n);nmin=nmin)
        time_i = toc()
        push!(ti,time_i)
    end
    return ( mean(ti) , ti )
end

function test_optimal_nmin( d::Vector{Int64} , n::Vector{Int64} , nmins::Vector{Int64} , reps )
    # results[i,j][k] contains results for d[i] , n[j], and min=nmins[k]
    results      = Array{Vector{Float64},2}(length(d),length(n))
    results_full = Array{Vector{Vector{Float64}},2}(length(d),length(n))

    for zi in 1:length(d)
        for zj in 1:length(n)
            results_ij = Vector{Float64}(length(nmins))
            full_results_ij = Vector{Vector{Float64}}(length(nmins))
            for zk in 1:length(nmins)
                (t_m,t_all) = f_benchmark( d[zi] , n[zj] , reps ; nmin=nmins[zk])
                results_ij[zk] = t_m
                full_results_ij[zk] = t_all
            end
            results[zi,zj]      = results_ij
            results_full[zi,zj] = full_results_ij
        end
    end
    return ( results , results_full )
end

(results , results_full) = test_optimal_nmin( [2,4,20] , [8000] , [48,64,80,96,128,256,512] , 4)

(results , results_full) = test_optimal_nmin( [2,4,20] , [1000,2000,40000,8000,16000] , [48,64,80,96,128,256,512] , 4)

results[2,2]


text:
include("test_RL.jl")

function solve_tspd_RL(x::Vector{Float64}, y::Vector{Float64}; n_samples=1, device="cpu")
    supported_sizes = [11, 15, 20, 50, 100]
    n_nodes = length(x)
    if !in(n_nodes, supported_sizes)
        error("Not supported size: n = $(n_nodes). RL support n ∈ $supported_sizes only.")
    end
    # input x, y: depot is x[1], y[1]

    # in Python RL model, depot is the last element
    xx = [x[2:end]; x[1]]
    yy = [y[2:end]; y[1]]

    # table = Dict([11 => 25, 20 => 35, 50 => 85, 100 => 180])
    # decode_len = table[n_nodes]
    decode_len = round(Int, n_nodes * 1.8)

    pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)

    py"""
    def solve_tspd(x, y, n_samples, device, decode_len):
        import main
        return main.main(x, y, n_samples, device, decode_len)
    """

    obj, route_t_, route_d_ = py"solve_tspd"(xx, yy, n_samples, device, decode_len)

    # +1 for python-julia indexing difference
    # another +1 for depot numbering difference
    route_t = Int.(route_t_) .+ 2
    route_d = Int.(route_d_) .+ 2

    # remove same nodes
    unique!(route_t)
    unique!(route_d)
    # add the origin depot
    pushfirst!(route_t, 1)
    pushfirst!(route_d, 1)

    # return obj, route_t, route_d

    Ct, Cd = cost_matrices_with_dummy(x, y, 1.0, 0.5)
    result = TSPDroneResult(obj, route_t, route_d, Ct, Cd, Inf)
    return result
end


text:
function mlaxor(attrib_mod)
    b=GeoPhyInv.Fields(attrib_mod)
    cols = ["red","blue","black","green","blue2","orange","brown","tan1","blue4","grey","green4","purple4","orange4","brown4"]
    i=1
    @pdf begin
    xlim=(0,100)
    zlim=(0,100)
        for c in b
        mgrid=GeoPhyInv.get_mgrid(eval(:(GeoPhyInv.$c())),attrib_mod, range(xlim[1], stop=xlim[2], step=10), range(zlim[1], stop=zlim[2], step=10))
        println(mgrid)
         setcolor(cols[i])
         println(cols[i])
        setmode("overlap")
            for z in mgrid[2], x in mgrid[1]
            circle(Point(z, x),1,:fill)
            # fontsize(2)
            # fontface("Georgia")
            # [Luxor.text(string(z),Point(z,x))]
        if (z==mgrid[2][1] && x==mgrid[1][1])
            fontsize(5)
            fontface("Georgia")
            [Luxor.text(string(c),Point(xlim[2]+10,zlim[1]+(i-1)*5))]
            end
        end
        i+=1
    end
    end 600 600
    finish()
    preview()
end


text:
using Kinetic, Plots, Flux
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress

function split_regime!(regime, regime0, ks, ctr, nn)
    @inbounds for i in eachindex(regime)
        if i == 0
            sw = (ctr[i+1].w .- ctr[i].w) / ks.ps.dx[i]
        else
            sw = (ctr[i].w .- ctr[i-1].w) / ks.ps.dx[i]
        end
        tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
        regime[i] = Int(round(nn([ctr[i].w[1:2]; ctr[i].w[end]; sw[1:2]; sw[end]; tau])[1]))

        if regime[i] == 1 && regime0[i] == 0
            Mu, Mv, Mw, _, _1 = gauss_moments(ctr[i].prim, ks.gas.K)
            a = pdf_slope(ctr[i].prim, sw, ks.gas.K)
            swt = -ctr[i].prim[1] .* moments_conserve_slope(a, Mu, Mv, Mw, 1, 0, 0)
            A = pdf_slope(ctr[i].prim, swt, ks.gas.K)
            tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
            ctr[i].f .= chapman_enskog(ks.vs.u, ks.vs.v, ks.vs.w, ctr[i].prim, a, zero(a), zero(a), A, tau)
            #ctr[i].f .= maxwellian(ks.vs.u, ks.vs.v, ks.vs.w, ctr[i].prim)
        end
    end

    regime0 .= regime

    return nothing
end

function kngll_regime!(regime, regime0, ks, ctr, nn)
    @inbounds for i in eachindex(regime)
        if i == 0
            sw = (ctr[i+1].w .- ctr[i].w) / ks.ps.dx[i]
        else
            sw = (ctr[i].w .- ctr[i-1].w) / ks.ps.dx[i]
        end

        L = abs(ctr[i].w[1] / sw[1])
        ℓ = (1/ctr[i].prim[end])^ks.gas.ω / ctr[i].prim[1] * sqrt(ctr[i].prim[end]) * ks.gas.Kn
        regime[i] = ifelse(ℓ / L >= 0.01, 1, 0)

        if regime[i] == 1 && regime0[i] == 0
            Mu, Mv, Mw, _, _1 = gauss_moments(ctr[i].prim, ks.gas.K)
            a = pdf_slope(ctr[i].prim, sw, ks.gas.K)
            swt = -ctr[i].prim[1] .* moments_conserve_slope(a, Mu, Mv, Mw, 1, 0, 0)
            A = pdf_slope(ctr[i].prim, swt, ks.gas.K)
            tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
            ctr[i].f .= chapman_enskog(ks.vs.u, ks.vs.v, ks.vs.w, ctr[i].prim, a, zero(a), zero(a), A, tau)
        end
    end

    regime0 .= regime

    return nothing
end

function up!(KS, ctr, regime, p)
    kn_bz, nm, phi, psi, phipsi = p

    res = zeros(5)
    avg = zeros(5)

    @inbounds for i = 1:ks.ps.nx
        if false#regime[i] == 0
            KitBase.step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                ks.gas.γ,
                ks.ps.dx[i],
                res,
                avg,
            )
        else
            KitBase.step!(
                face[i].fw,
                face[i].ff,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[i+1].fw,
                face[i+1].ff,
                KS.gas.γ,
                kn_bz,
                nm,
                phi,
                psi,
                phipsi,
                KS.ps.dx[i],
                dt,
                res,
                avg,
                :fsm,
            )
        end
    end
end

cd(@__DIR__)
#@load "../nn_scalar.jld2" nn
@load "/home2/vavrines/Coding/Flowmachine/src/1d/sampler/nn_rif.jld2" nn

set = Setup(case = "sod", space = "1d1f3v", collision = "boltzmann", maxTime = 0.1, cfl = 0.3)
ps = PSpace1D(0.0, 1.0, 100, 1)
vs = VSpace3D(-8.0, 8.0, 80, -8.0, 8.0, 28, -8.0, 8.0, 28)
gas = Gas(Kn = 1e-3, Pr = 2/3, K = 0.0)
ib = IB1F(ib_sod(set, ps, vs, gas)...)
ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks, ks.ps, :dynamic_array; structarray = false)

begin
    kn_bzm = hs_boltz_kn(ks.gas.μᵣ, 1.0)
    phi, psi, phipsi = kernel_mode(
        5,
        ks.vs.u1,
        ks.vs.v1,
        ks.vs.w1,
        ks.vs.du[1, 1, 1],
        ks.vs.dv[1, 1, 1],
        ks.vs.dw[1, 1, 1],
        ks.vs.nu,
        ks.vs.nv,
        ks.vs.nw,
        1.0,
    )
    p = (kn_bzm, 5, phi, psi, phipsi)
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ctr[1].w)
regime = ones(Int, axes(ks.ps.x))
regime0 = deepcopy(regime)

@time for iter = 1:nt
    println("iteration: $iter of $nt")

    split_regime!(regime, regime0, ks, ctr, nn)
    #kngll_regime!(regime, regime0, ks, ctr, nn)

    @inbounds for i in eachindex(face)
        if false#1 ∉ (regime[i], regime[i-1])
            flux_gks!(
                face[i].fw,
                ctr[i-1].w .+ ctr[i-1].sw .* ks.ps.dx[i-1] / 2,
                ctr[i].w .- ctr[i].sw .* ks.ps.dx[i] / 2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1] / 2,
                ks.ps.dx[i] / 2,
                1.0,
                ctr[i-1].sw,
                ctr[i].sw,
            )
        else
            flux_kfvs!(
                face[i].fw,
                face[i].ff,
                ctr[i-1].f .+ ctr[i-1].sf .* ks.ps.dx[i-1] / 2,
                ctr[i].f .- ctr[i].sf .* ks.ps.dx[i] / 2,
                ks.vs.u,
                ks.vs.v,
                ks.vs.w,
                ks.vs.weights,
                dt,
                ctr[i-1].sf,
                ctr[i].sf,
            )
        end
    end

    #KitBase.update!(ks, ctr, face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))
    up!(ks, ctr, regime, p)

    t += dt
end

plot(ks, ctr)
plot!(ks.ps.x[1:ks.ps.nx], regime[1:ks.ps.nx])

"""
Benchmark computation cost:

Kn = 1e-4
dvm: 1211.636247 seconds (86.95 M allocations: 4.535 TiB, 5.42% gc time, 0.36% compilation time)
adaptive: 377.118386 seconds (50.25 M allocations: 1.383 TiB, 4.01% gc time, 1.12% compilation time)

Kn = 1e-3
dvm: 1172.045585 seconds (63.84 M allocations: 4.534 TiB, 6.01% gc time)
adaptive: 382.871913 seconds (27.70 M allocations: 1.426 TiB, 6.65% gc time)

Kn = 1e-2
dvm: 1177.973543 seconds (63.90 M allocations: 4.534 TiB, 6.02% gc time, 0.00% compilation time)
adaptive: 741.156948 seconds (43.94 M allocations: 2.819 TiB, 4.68% gc time)

"""

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
min_creation(file) = unix2datetime(min(mtime(file), ctime(file)))

function tryparsedatetime(x)
    try
        dt = DateTime(x)
        return dt
    catch
        return nothing
    end
end

function getdatetime(file, minimum_dt)
    @label askdatetime
    dialog = Dialog(string(min_creation(joinpath(coffeesource, file))))
    _dt = requestᵐ("Specify a creation date & time for: $file", dialog)
    dt = tryparsedatetime(_dt)
    if dt ≡ nothing
        @warn "the format of the date or time is wrong, try something like:" now()
        @goto askdatetime
    end
    if dt < minimum_dt
        @warn "the date & time should come after the end of the previous video segment" minimum_dt
        @goto askdatetime
    end
    dt
end

function getorder(_files)
    length(_files) ≤ 1 && return _files
    perms = permutations(_files)
    options = [join(files, ", ") for files in perms]
    menu = RadioMenu(options)
    i = request("Which chronological order of the video files is correct?", menu)
    first(Iterators.drop(perms, i - 1))
end

const videocomment_dialog = Dialog()
const video_type_menu = RadioMenu(["segmented: multiple segments with no temporal gaps", "disjointed: video files with temppral gaps between them"], 1)

function _newvideo(files)
    @label newvideol
    @info "Registering video file/s:" files
    sunregistered = getorder(files)
    durations = [Millisecond(round(Int, 1000VideoIO.get_duration(joinpath(coffeesource, file_name)) )) for file_name in sunregistered]
    comment = requestᵐ("Any comments about this video as a whole?", videocomment_dialog)
    nfiles = length(files)
    if nfiles == 1
        date_time = getdatetime(sunregistered[], Date(0))
        try
            return WholeVideo(VideoFile(files[], date_time, durations[]), comment)
        catch ex
            @warn "something was wrong" ex
            @goto newvideol
        end
    else
        i = requestᵐ("Is this video segmented or disjointed?", video_type_menu)
        date_times = Vector{DateTime}(undef, nfiles)
        if i == 1
            date_times[1] = getdatetime(sunregistered[1], Date(0))
            for i in 2:nfiles
                date_times[i] = date_times[i - 1] + durations[i - 1]
            end
            try
                return FragmentedVideo(VideoFile.(sunregistered, date_times, durations), comment)
            catch ex
                @warn "something was wrong" ex
                @goto newvideol
            end
        else
            last_dt = DateTime(0)
            for (i, file_name) in enumerate(sunregistered)
                date_times[i] = getdatetime(file_name, last_dt)
                last_dt = date_times[i] + durations[i]
            end
            try
                return DisjointVideo(VideoFile.(sunregistered, date_times, durations), comment)
            catch ex
                @warn "something was wrong" ex
                @goto newvideol
            end
        end
    end
end

function selectfiles(files)
    length(files) == 1 && return files
    menu = MultiSelectMenu(files, Set(collect(1:length(files))))
    i = request("Select the file, or multiple files in case of a segmented video, that constitute/s a single video:", menu)
    files[collect(i)]
end

goodvideo(file) = first(file) ≠ '.' && occursin(r"mts|mp4|avi|mpg|mov|mkv"i, last(splitext(file))) && isfile(joinpath(coffeesource, file))

function newvideo(existing)
    files = String[file for file in readdir(coffeesource) if file ∉ existing && goodvideo(file)]
    if isempty(files)
        @warn "found no new unregistered video files…"
        return nothing
    end
    files = selectfiles(files)
    if isempty(files)
        @warn "no files were selected for registration…"
        return nothing
    end
    _newvideo(files)
end

function newtl(tl::WholeVideo, _, start)
    file = tl.file
    WholeVideo(tl, (; start = start))
end
function newtl(tl::T, j, start) where T <: AbstractTimeLine
    tl.files[j] = VideoFile(tl.files[j], (; start = start))
    tl
end

function edit_video()
    i = request("Which video file do you want to edit the registration of?", videofile_menu)
    vf = videofiles[i]
    video, ii, j = findvideoindex(vf)
    @info """the current date & time of "$(vf.name)" is:""" start(vf)
    st = getdatetime(vf.name, Date(0))
    new_video = newtl(video, ii, st)
    if video == new_video
        @info "nothing changed"
    else
        videofiles[i] = files(new_video)[ii]
        videos[j] = new_video
        @info "changes are saved to file"
        JLSO.save(joinpath(sourcefolder, "videos.jlso"), videos)
    end
end


text:
#MAIN ENTRY POINT
#final columns for the decompressed file

#merges crsp and compustat
# note refresh merge will be automatically run if any of the other two are run
function loadccm(;
  datapath::String=DATA_PATH,
  ccmname::String = CCM_NAME,
  incsvstream::Function = IN_CSV_STREAM,
  csvextension::String = CSV_EXTENSION)

  local ccm::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  ccm = incsvstream("$datapath\\$ccmname.$csvextension") |> CSV.read |> DataFrame

  return ccm

end

function mergecompccm!(comp::DataFrame;
    testoutput::Bool = TEST_OUTPUT)

  #first
  local ccm::DataFrame = loadccm()
  ccm = preprocessccm!(ccm)

  testoutput && CSV.write("output\\ccm.csv", ccm)

  select!(ccm, Not([:cusip]))
  #println("Comp rows before merge: $(size(comp, 1))")
  comp = join(comp, ccm, on=[:gvkey], kind=:inner)
  comp.keep = trues(size(comp,1))

  comp.keep = (r::DataFrameRow->
    r.adate ∈ r.linkeffdate:Day(1):r.linkenddate).(eachrow(comp))

  comp = comp[comp.keep,:]
  select!(comp, Not(:keep))
  println("Comp rows after merge, post-filter: $(size(comp, 1))")

  comp = reconcilecompccmintervals!(comp)

  return comp
end

#formats the ccm table in a reasonable way
function preprocessccm!(ccm::DataFrame;
  retainedcolsccm::Vector{Symbol} = RETAINED_COLS_CCM,
  lastdatestring::String = LAST_DATE_STRING,
  testoutput::Bool = TEST_OUTPUT)

  #names ot lower case
  cleannames::Vector{Symbol} = (s::Symbol->Symbol(lowercase(string(s)))).(names(ccm))
  rename!(ccm, ((oldn, newn)->oldn=>newn).(names(ccm), cleannames))

  #fix the end date
  ccm.linkenddt .= (s::MString->
    coalesce(s,"")=="E" ? lastdatestring : s).(ccm.linkenddt)

  #deleterows!(ccm, (!).(completecases(ccm[!,[:gvkey, :lpermno]])))

  #below didn't work for some reason
  #filter!(r::DataFrameRow->((!ismissing(r.gvkey)) && (!ismissing(r.lpermno))), ccm)
  ccm = ccm[(r::DataFrameRow->(
    (!ismissing(r.gvkey)) && (!ismissing(r.lpermno)))).(eachrow(ccm)),:]
  ccm.linkprim = (s::String->length(s)==1 ? s[1] : error("invalid linkprim in ccm")).(ccm.linkprim)
  filter!(r::DataFrameRow->r.linkprim =='P' || r.linkprim =='C', ccm)

  #parse the dates
  ccm.linkeffdate = (s->parseccm(Date, s)).(ccm.linkdt)
  ccm.linkenddate = (s->parseccm(Date, s)).(ccm.linkenddt)

  #dedup if possible
  ccm.tokeep = trues(size(ccm,1))
  ccm.dateranges = ((linkeffdate::Date,linkenddate::Date)->
    linkeffdate:Day(1):linkenddate).(ccm.linkeffdate,ccm.linkenddate)

  gccm::GroupedDataFrame = groupby(ccm, [:gvkey, :lpermno])
  @mpar for i ∈ 1:length(gccm)
    sccm::SubDataFrame = gccm[i]
    Nsccm::Int = size(sccm,1)

    #construct the union of all of the date sets
    if Nsccm >  1
      local periodrange::StepRange{Date,Day} = minimum(
        sccm.linkeffdate):Day(1):maximum(sccm.linkenddate)
      local dateunionlength::Int = length(unique([sccm.dateranges...;]))

      if dateunionlength == length(periodrange)
        sccm.linkeffdate[1] = minimum(periodrange)
        sccm.linkenddate[1] = maximum(periodrange)
        sccm.tokeep[2:end] .= false
      end

    end
  end

  println("ccm size pre-dedup $(size(ccm))")
  ccm = ccm[ccm.tokeep,:]
  println("ccm size post-dedup $(size(ccm))")

  if testoutput
    testoutput && CSV.write("output\\ccm_preproc.csv", ccm)
  end
  #rename!(ccm, :indfmt=>:indfmtold)
  #ccm.indfmt = (s::MString->parsecomp(Symbol, s)).(ccm.indfmtold)

  #only need some of the fields
  ccm = select!(ccm, retainedcolsccm)

  return ccm
end

#this function reconciles the link intervals with the announcement intervals
#see
function reconcilecompccmintervals!(comp::DataFrame)
  #subtract one day since the effective interval is inclusive of the
  #first date while the begin-end date is exclusve of the first date
  comp.begindate .= (max).(comp.begindate, comp.linkeffdate .- Day(1))
  comp.enddate .= (min).(comp.enddate, comp.linkenddate)

  #this gets rid of one day records
  #(this assumes the previous effective record is the current record)
  #not many of these ~140
  #println("records before interval reconciliation: $(size(comp,1))")
  #println(comp[comp.begindate .≥ comp.enddate,
  #  [:gvkey, :adateq, :begindate, :enddate, :linkeffdate, :linkenddate]])
  comp = comp[comp.begindate .< comp.enddate, :]
  println("records after interval reconciliation: $(size(comp,1))")

  return comp
end

function parseccm(::Type{T}, s::String;
    #=parsedmissings = CCM_PARSED_MISSINGS=#) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  #(!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

#the date case
parseccm(::Type{Date}, s::String;
    ccmdateformat::DateFormat = CCM_DATE_FORMAT) = Dates.Date(s, ccmdateformat)

#helper methods and special cases
parseccm(::Type{<:Any}, v::Missing) = missing
parseccm(::Type{Date}, i::Int) = parseccm(Date, "$i")
parseccm(::Type{Symbol}, s::String) = Symbol(s)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
struct IndicatorBridge{T, B<:MOIBC.SetMapBridge{T}, A, IS_MOI_OR_CS, S} <: MOIBC.AbstractBridge
    con_idx::CI
end

function MOI.supports_constraint(
    ::Type{<:IndicatorBridge{T, B}},
    ::Type{F},
    ::Type{<:Union{MOI.IndicatorSet{A,S}, CS.IndicatorSet{A,S}}}
) where {T, B, F<:MOI.VectorAffineFunction, A, S}
    is_supported = MOI.supports_constraint(B, MOIU.scalar_type(F), S)
    !is_supported && return false
    S <: BoolSet && return true

    concrete_B = MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    return supports_concreteB(concrete_B)
end

function MOIBC.concrete_bridge_type(
    ::Type{<:IndicatorBridge{T,B}},
    G::Type{<:MOI.VectorAffineFunction},
    ::Type{IS},
) where {T,B,A,S,IS<:Union{MOI.IndicatorSet{A,S}, CS.IndicatorSet{A,S}}}
    if S <: BoolSet
        concrete_B = B
    else
        concrete_B = MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    end
    IS_MOI_OR_CS = IS <: MOI.IndicatorSet ? MOI.IndicatorSet : CS.IndicatorSet
    return IndicatorBridge{T,concrete_B,A,IS_MOI_OR_CS,S}
end

function MOIB.added_constraint_types(
    ::Type{<:IndicatorBridge{T,B,A,IS_MOI_OR_CS,S}}
) where {T,B,A,IS_MOI_OR_CS,S}
    if S <: BoolSet
        added_constraints = added_constraint_types(B, S)
    else
        added_constraints = MOIB.added_constraint_types(B)
    end
    return [(MOI.VectorAffineFunction{T}, IS_MOI_OR_CS{A,added_constraints[1][2]})]
end

function MOIB.added_constrained_variable_types(::Type{<:IndicatorBridge{T,B}}) where {T,B}
    return MOIB.added_constrained_variable_types(B)
end

function MOIBC.bridge_constraint(::Type{<:IndicatorBridge{T, B, A, IS_MOI_OR_CS, S}}, model, func, set) where {T, B, A, IS_MOI_OR_CS, S}
    f = MOIU.eachscalar(func)
    if S <: BoolSet
        new_func = MOIU.operate(vcat, T, f[1], map_function(B, f[2:end], set.set))
        new_inner_set = map_set(B, set.set)
    else
        new_func = MOIU.operate(vcat, T, f[1], MOIBC.map_function(B, f[2:end]))
        new_inner_set = MOIBC.map_set(B, set.set)
    end
    new_set = IS_MOI_OR_CS{A}(new_inner_set)
    return IndicatorBridge{T,B,A,IS_MOI_OR_CS,S}(MOI.add_constraint(model, new_func, new_set))
end


text:
using Test
using BlackScholesGreeks.GreekAnalytical

@testset "GreekAnalytical.jl" begin
    # Delta Test Case 1
    test_delta = blackScholesDelta(true, 105.0, 100.0, 0.5, 0.1, 0.0, 0.36)
    @test test_delta ≈ 0.5946 atol=1e-4
    # Delta Test Case 2
    test_delta = blackScholesDelta(false, 105.0, 100.0, 0.5, 0.1, 0.0, 0.36)
    @test test_delta ≈ -0.3566 atol=1e-4
    # Vanna Test Case 1
    test_vanna = blackScholesVanna(false, 90.0, 80.0, 0.25, 0.05, 0.05, 0.2)
    @test test_vanna ≈ -1.0008 atol=1e-4
    # Charm Test Case 1
    test_charm = blackScholesCharm(false, 105.0, 90.0, 0.25, 0.14, 0.0, 0.24)
    @test test_charm ≈ 0.3700 atol=1e-4
    # Gamma Test Case 1
    test_gamma = blackScholesGamma(true, 55.0, 60.0, 0.75, 0.1, 0.1, 0.3)
    @test test_gamma ≈ 0.0278 atol=1e-4
    # Zomma Test Case 1
    test_zomma = blackScholesZomma(false, 100.0, 80.0, 0.25, 0.05, 0.0, 0.26)
    @test test_zomma ≈ 0.0463 atol=1e-4
    # Speed Test Case 1
    test_speed = blackScholesSpeed(false, 50.0, 48.0, 0.0833, 0.06, 0.01, 0.2)
    @test test_speed ≈ -0.0291 atol=1e-4
    # Vega Test Case 1
    test_vega = blackScholesVega(true, 55.0, 60.0, 0.75, 0.105, 0.0695, 0.3)
    @test test_vega ≈ 18.5027 atol=1e-4
    # Volga Test Case 1
    test_volga = blackScholesVolga(false, 90.0, 130.0, 0.75, 0.05, 0.0, 0.28)
    @test test_volga ≈ 92.3444 atol=1e-4
    # Theta Test Case 1
    test_theta = blackScholesTheta(false, 430.0, 405.0, 0.0833, 0.07, 0.02, 0.2)
    @test test_theta ≈ -31.1924 atol=1e-4
    # Driftless Theta Test Case 1
    test_driftless_theta = blackScholesDriftlessTheta(false, 430.0, 405.0, 0.0833, 0.07, 0.02, 0.2)
    @test test_driftless_theta ≈ -32.6214 atol=1e-4
    # Rho Test Case 1
    test_rho = blackScholesRho(true, 72.0, 75.0, 1.0, 0.09, 0.09, 0.19)
    @test test_rho ≈ 38.7325 atol=1e-4
end

max_stars_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_issues_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_forks_repo_name: Stochastic-Adventure/BlackScholesGreeks


text:
"""
  `imadjustintensity(img [, (minval,maxval)]) -> Image`

   Map intensities over the interval `(minval,maxval)` to the interval
   `[0,1]`. This is equivalent to `map(ScaleMinMax(eltype(img), minval,
   maxval), img)`.  (minval,maxval) defaults to `extrema(img)`.
"""
imadjustintensity{T}(img::AbstractArray{T}, range) = map(ScaleMinMax(T, range...), img)
imadjustintensity{T}(img::AbstractArray{T}) = map(ScaleAutoMinMax(T), img)

_imstretch{T}(img::AbstractArray{T}, m::Number, slope::Number) = shareproperties(img, map(i -> 1 / (1 + (m / (i + eps(T))) ^ slope), img))

"""
`imgs = imstretch(img, m, slope)` enhances or reduces (for
slope > 1 or < 1, respectively) the contrast near saturation (0 and 1). This is
essentially a symmetric gamma-correction. For a pixel of brightness `p`, the new
intensity is `1/(1+(m/(p+eps))^slope)`.

This assumes the input `img` has intensities between 0 and 1.
"""
imstretch(img::AbstractArray, m::Number, slope::Number) = _imstretch(float(img), m, slope)

"""
```
complement_img = imcomplement(img)
```

Returns the complement of an image.
"""
imcomplement{T}(img::AbstractArray{T}) = map(complement, img)

imcomplement(img::AbstractImage) = copyproperties(img, imcomplement(data(img)))
complement(x) = one(x)-x
complement(x::TransparentColor) = typeof(x)(complement(color(x)), alpha(x))

imhist{T<:Colorant}(img::AbstractArray{T}, nbins::Integer = 400) = imhist(convert(Array{Gray}, data(img)), nbins)

function imhist{T<:Union{Gray,Number}}(img::AbstractArray{T}, nbins::Integer = 400)
    minval = minfinite(img)
    maxval = maxfinite(img)
    imhist(img, nbins, minval, maxval)
end

"""
```
edges, count = imhist(img, nbins)
edges, count = imhist(img, nbins, minval, maxval)
```

Generates a histogram for the image over nbins spread between `(minval, maxval]`.
If `minval` and `maxval` are not given, then the minimum and
maximum values present in the image are taken.

`edges` is a vector that specifies how the range is divided;
`count[i+1]` is the number of values `x` that satisfy `edges[i] <= x < edges[i+1]`.
`count[1]` is the number satisfying `x < edges[1]`, and
`count[end]` is the number satisfying `x >= edges[end]`. Consequently,
`length(count) == length(edges)+1`.
"""
function imhist(img::AbstractArray, nbins::Integer, minval::Union{Gray,Real}, maxval::Union{Gray,Real})
    edges = StatsBase.histrange([Float64(minval), Float64(maxval)], nbins, :left)
    imhist(img, edges)
end

function imhist(img::AbstractArray, edges::Range)
    histogram = zeros(Int, length(edges) + 1)
    o = Base.Order.Forward
    G = graytype(eltype(img))
    for v in img
        val = convert(G, v)
        if val >= edges[end]
            histogram[end] += 1
            continue
        end
        index = searchsortedlast(edges, val, o)
        histogram[index + 1] += 1
    end
    edges, histogram
end

function _histeq_pixel_rescale{T<:Union{Gray,Number}}(pixel::T, cdf, minval, maxval)
    n = length(cdf)
    bin_pixel = clamp(ceil(Int, (pixel - minval) * length(cdf) / (maxval - minval)), 1, n)
    rescaled_pixel = minval + ((cdf[bin_pixel] - cdf[1]) * (maxval - minval) / (cdf[end] - cdf[1]))
    convert(T, rescaled_pixel)
end
function _histeq_pixel_rescale{C<:Color}(pixel::C, cdf, minval, maxval)
    yiq = convert(YIQ, pixel)
    y = _histeq_pixel_rescale(yiq.y, cdf, minval, maxval)
    convert(C, YIQ(y, yiq.i, yiq.q))
end
function _histeq_pixel_rescale{C<:TransparentColor}(pixel::C, cdf, minval, maxval)
    base_colorant_type(C)(_histeq_pixel_rescale(color(pixel), cdf, minval, maxval), alpha(pixel))
end

"""
```
hist_equalised_img = histeq(img, nbins)
hist_equalised_img = histeq(img, nbins, minval, maxval)
```

Returns a histogram equalised image with a granularity of approximately `nbins`
number of bins.

The `histeq` function can handle a variety of input types. The returned image depends
on the input type. If the input is an `Image` then the resulting image is of the same type
and has the same properties.

For coloured images, the input is converted to YIQ type and the Y channel is equalised. This
is the combined with the I and Q channels and the resulting image converted to the same type
as the input.

If minval and maxval are specified then intensities are equalized to the range
(minval, maxval). The default values are 0 and 1.

"""
function histeq(img::AbstractArray, nbins::Integer, minval::Union{Number,Gray}, maxval::Union{Number,Gray})
    bins, histogram = imhist(img, nbins, minval, maxval)
    cdf = cumsum(histogram[2:end-1])
    img_shape = size(img)
    minval == maxval && return map(identity, img)
    # Would like to use `map` here, but see https://github.com/timholy/Images.jl/pull/523#issuecomment-235236460
    hist_equalised_img = similar(img)
    for I in eachindex(img)
        hist_equalised_img[I] = _histeq_pixel_rescale(img[I], cdf, minval, maxval)
    end
    hist_equalised_img
end

function histeq(img::AbstractArray, nbins::Integer)
    T = graytype(eltype(img))
    histeq(img, nbins, zero(T), one(T))
end

function histeq(img::AbstractImage, nbins::Integer, minval::Union{Number,Gray}, maxval::Union{Number,Gray})
    newimg = histeq(data(img), nbins, minval, maxval)
    shareproperties(img, newimg)
end

histeq(img::AbstractImage, nbins::Integer) = shareproperties(img, histeq(data(img), nbins))

adjust_gamma(img::AbstractImage, gamma::Number) = shareproperties(img, adjust_gamma(data(img), gamma))

_gamma_pixel_rescale{T<:Union{Gray, Number}}(pixel::T, gamma::Number) = pixel ^ gamma

function _gamma_pixel_rescale{C<:Color}(pixel::C, gamma::Number)
    yiq = convert(YIQ, pixel)
    y = _gamma_pixel_rescale(yiq.y, gamma)
    convert(C, YIQ(y, yiq.i, yiq.q))
end

function _gamma_pixel_rescale{C<:TransparentColor}(pixel::C, gamma::Number)
    base_colorant_type(C)(_gamma_pixel_rescale(color(pixel), gamma), alpha(pixel))
end

function _gamma_pixel_rescale(original_val::Number, gamma::Number, minval::Number, maxval::Number)
    Float32(minval + (maxval - minval) * ((original_val - minval) / (maxval - minval)) ^ gamma)
end

"""
```
gamma_corrected_img = adjust_gamma(img, gamma)
```

Returns a gamma corrected image.

The `adjust_gamma` function can handle a variety of input types. The returned image depends
on the input type. If the input is an `Image` then the resulting image is of the same type
and has the same properties.

For coloured images, the input is converted to YIQ type and the Y channel is gamma corrected.
This is the combined with the I and Q channels and the resulting image converted to the same
type as the input.

"""
function adjust_gamma{T<:FixedPointNumbers.UFixed}(img::AbstractArray{Gray{T}}, gamma::Number)
    raw_type = FixedPointNumbers.rawtype(T)
    gamma_inv = 1.0 / gamma
    table = zeros(T, typemax(raw_type) + 1)
    for i in zero(raw_type):typemax(raw_type)
        table[i + 1] = T((i / typemax(raw_type)) ^ gamma_inv)
    end
    gamma_corrected_img = similar(img)
    for I in eachindex(img)
        gamma_corrected_img[I] = Gray(table[convert(base_colorant_type(typeof(img[I])){T}, img[I]).val.i + 1])
    end
    gamma_corrected_img
end

adjust_gamma{T<:Number}(img::AbstractArray{T}, gamma::Number) = _adjust_gamma(img, gamma, Float64)
adjust_gamma{T<:Colorant}(img::AbstractArray{T}, gamma::Number) = _adjust_gamma(img, gamma, T)

function _adjust_gamma(img::AbstractArray, gamma::Number, C::Type)
    gamma_corrected_img = zeros(C, size(img))
    for I in eachindex(img)
        gamma_corrected_img[I] = _gamma_pixel_rescale(img[I], gamma)
    end
    gamma_corrected_img
end

function adjust_gamma{T<:Number}(img::AbstractArray{T}, gamma::Number, minval::Number, maxval::Number)
    gamma_corrected_img = zeros(Float64, size(img))
    for I in eachindex(img)
        gamma_corrected_img[I] = _gamma_pixel_rescale(img[I], gamma, minval, maxval)
    end
    gamma_corrected_img
end

"""
```
hist_matched_img = histmatch(img, oimg, nbins)
```

Returns a grayscale histogram matched image with a granularity of `nbins` number of bins. `img` is the image to be
matched and `oimg` is the image having the desired histogram to be matched to.

"""
histmatch(img::AbstractImage, oimg::AbstractArray, nbins::Integer = 400) = shareproperties(img, histmatch(data(img), oimg, nbins))

_hist_match_pixel{T<:Union{Gray, Number}}(pixel::T, bins, lookup_table) = T(bins[lookup_table[searchsortedlast(bins, pixel)]])

function _hist_match_pixel{T<:Color}(pixel::T, bins, lookup_table)
    yiq = convert(YIQ, pixel)
    y = _hist_match_pixel(yiq.y, bins, lookup_table)
    convert(T, YIQ(y, yiq.i, yiq.q))
end

_hist_match_pixel{T<:TransparentColor}(pixel::T, bins, lookup_table) = base_colorant_type(T)(_hist_match_pixel(color(pixel), bins, lookup_table), alpha(pixel))

function histmatch{T<:Colorant}(img::AbstractArray{T}, oimg::AbstractArray, nbins::Integer = 400)
    el_gray = graytype(eltype(img))
    oedges, ohist = imhist(oimg, nbins, zero(el_gray), one(el_gray))
    _histmatch(img, oedges, ohist)
end

function _histmatch(img::AbstractArray, oedges::Range, ohist::AbstractArray{Int})
    bins, histogram = imhist(img, oedges)
    ohist[1] = zero(eltype(ohist))
    ohist[end] = zero(eltype(ohist))
    histogram[1] = zero(eltype(histogram))
    histogram[end] = zero(eltype(histogram))
    cdf = cumsum(histogram)
    norm_cdf = cdf / cdf[end]
    ocdf = cumsum(ohist)
    norm_ocdf = ocdf / ocdf[end]
    lookup_table = zeros(Int, length(norm_cdf))
    @compat for I in eachindex(cdf)
        lookup_table[I] = indmin(abs.(norm_ocdf .- norm_cdf[I]))
    end
    hist_matched_img = similar(img)
    for I in eachindex(img)
        hist_matched_img[I] = _hist_match_pixel(img[I], bins, lookup_table)
    end
    hist_matched_img
end

"""
```
hist_equalised_img = clahe(img, nbins, xblocks = 8, yblocks = 8, clip = 3)

```

Performs Contrast Limited Adaptive Histogram Equalisation (CLAHE) on the input image. It differs from ordinary histogram
equalization in the respect that the adaptive method computes several histograms, each corresponding to a distinct section
of the image, and uses them to redistribute the lightness values of the image. It is therefore suitable for improving the
local contrast and enhancing the definitions of edges in each region of an image.

In the straightforward form, CLAHE is done by calculation a histogram of a window around each pixel and using the transformation
function of the equalised histogram to rescale the pixel. Since this is computationally expensive, we use interpolation which gives
a significant rise in efficiency without compromising the result. The image is divided into a grid and equalised histograms are
calculated for each block. Then, each pixel is interpolated using the closest histograms.

The `xblocks` and `yblocks` specify the number of blocks to divide the input image into in each direction. `nbins` specifies
the granularity of histogram calculation of each local region. `clip` specifies the value at which the histogram is clipped.
The excess in the histogram bins with value exceeding `clip` is redistributed among the other bins.

"""
function clahe{C}(img::AbstractArray{C, 2}, nbins::Integer = 100; xblocks::Integer = 8, yblocks::Integer = 8, clip::Number = 3)
    h, w = size(img)
    y_padded = ceil(Int, h / (2 * yblocks)) * 2 * yblocks
    x_padded = ceil(Int, w / (2 * xblocks)) * 2 * xblocks

    img_padded = imresize(img, (y_padded, x_padded))

    hist_equalised_img = _clahe(img_padded, nbins, xblocks, yblocks, clip)
    imresize(hist_equalised_img, (h, w))
end

function clahe(img::AbstractImage, nbins::Integer = 100; xblocks::Integer = 8, yblocks::Integer = 8, clip::Number = 3)
    shareproperties(clahe(data(img), nbins, xblocks = xblocks, yblocks = yblocks, clip = clip), img)
end

function _clahe{C}(img::AbstractArray{C, 2}, nbins::Integer = 100, xblocks::Integer = 8, yblocks::Integer = 8, clip::Number = 3)
    h, w = size(img)
    xb = 0:1:xblocks - 1
    yb = 0:1:yblocks - 1
    blockw = Int(w / xblocks)
    blockh = Int(h / yblocks)
    temp_cdf = Array{Float64, 1}[]
    T = graytype(eltype(img))
    edges = StatsBase.histrange([Float64(zero(T)), Float64(one(T))], nbins, :left)

    for i in xb
        for j in yb
            temp_block = img[j * blockh + 1 : (j + 1) * blockh, i * blockw + 1 : (i + 1) * blockw]
            _, histogram = imhist(temp_block, edges)
            clipped_hist = cliphist(histogram[2:end - 1], clip)
            cdf = cumsum(clipped_hist)
            n_cdf = cdf / cdf[end]
            push!(temp_cdf, n_cdf)
        end
    end

    norm_cdf = reshape(temp_cdf, yblocks, xblocks)
    res_img = zeros(C, size(img))

    #Interpolations
    xb = 1:2:xblocks * 2 - 2
    yb = 1:2:yblocks * 2 - 2

    for j in yb
        for i in xb
            p_block = img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)]
            bnum_y = floor(Int, (j + 1) / 2)
            bnum_x = floor(Int, (i + 1) / 2)
            top_left = norm_cdf[bnum_y, bnum_x]
            top_right = norm_cdf[bnum_y, bnum_x + 1]
            bot_left = norm_cdf[bnum_y + 1, bnum_x]
            bot_right = norm_cdf[bnum_y + 1, bnum_x + 1]
            temp_block = zeros(C, size(p_block))
            for l in 1:blockw
                for m in 1:blockh
                    temp_block[m, l] = _clahe_pixel_rescale(p_block[m, l], top_left, top_right, bot_left, bot_right, edges, l, m, blockw, blockh)
                end
            end
            res_img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)] = temp_block
        end
    end

    #Corners

    block_tl = img[1 : Int(blockh / 2), 1 : Int(blockw / 2)]
    corner_tl = map(i -> _clahe_pixel_rescale(i, norm_cdf[1, 1], edges), block_tl)
    res_img[1 : Int(blockh / 2), 1 : Int(blockw / 2)] = corner_tl

    block_tr = img[1 : Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)]
    corner_tr = map(i -> _clahe_pixel_rescale(i, norm_cdf[1, xblocks], edges), block_tr)
    res_img[1 : Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)] = corner_tr

    block_bl = img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), 1 : Int(blockw / 2)]
    corner_bl = map(i -> _clahe_pixel_rescale(i, norm_cdf[yblocks, 1], edges), block_bl)
    res_img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), 1 : Int(blockw / 2)] = corner_bl

    block_br = img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)]
    corner_br = map(i -> _clahe_pixel_rescale(i, norm_cdf[yblocks, xblocks], edges), block_br)
    res_img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)] = corner_br

    #Horizontal Borders

    for j in [0, yblocks * 2 - 1]
        for i in xb
            p_block = img[j * Int(blockh / 2) + 1 : (j + 1) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)]
            bnum_x = floor(Int, (i + 1) / 2)
            left = norm_cdf[ceil(Int, (j + 1) / 2), bnum_x]
            right = norm_cdf[ceil(Int, (j + 1) / 2), bnum_x + 1]
            temp_block = zeros(C, size(p_block))
            for l in 1:blockw
                for m in 1:Int(blockh / 2)
                    temp_block[m, l] = _clahe_pixel_rescale(p_block[m, l], left, right, edges, l, blockw)
                end
            end
            res_img[j * Int(blockh / 2) + 1 : (j + 1) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)] = temp_block
        end
    end

    #Vertical Borders

    for i in [0, xblocks * 2 - 1]
        for j in yb
            p_block = img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 1) * Int(blockw / 2)]
            bnum_y = floor(Int, (j + 1) / 2)
            top = norm_cdf[bnum_y, ceil(Int, (i + 1) / 2)]
            bot = norm_cdf[bnum_y + 1, ceil(Int, (i + 1) / 2)]
            temp_block = zeros(C, size(p_block))
            for l in 1:Int(blockw / 2)
                for m in 1:blockh
                    temp_block[m, l] = _clahe_pixel_rescale(p_block[m, l], top, bot, edges, m, blockh)
                end
            end
            res_img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 1) * Int(blockw / 2)] = temp_block
        end
    end
    res_img
end

_clahe_pixel_rescale{T<:Union{Gray, Number}}(pixel::T, cdf, edges) = cdf[searchsortedlast(edges, pixel, Base.Order.Forward)]

function _clahe_pixel_rescale{T<:Union{Gray, Number}}(pixel::T, first, second, edges, pos, length)
    id = searchsortedlast(edges, pixel, Base.Order.Forward)
    f = first[id]
    s = second[id]
    T(((length - pos) * f + (pos - 1) * s) / (length - 1))
end

function _clahe_pixel_rescale{T<:Union{Gray, Number}}(pixel::T, top_left, top_right, bot_left, bot_right, edges, i, j, w, h)
    id = searchsortedlast(edges, pixel, Base.Order.Forward)
    tl = top_left[id]
    tr = top_right[id]
    bl = bot_left[id]
    br = bot_right[id]
    r1 = ((w - i) * tl + (i - 1) * tr) / (w - 1)
    r2 = ((w - i) * bl + (i - 1) * br) / (w - 1)
    T(((h - j) * r1 + (j - 1) * r2) / (h - 1))
end

function _clahe_pixel_rescale{C<:Color}(pixel::C, args...)
    yiq = convert(YIQ, pixel)
    y = _clahe_pixel_rescale(yiq.y, args...)
    convert(C, YIQ(y, yiq.i, yiq.q))
end

function _clahe_pixel_rescale{C<:TransparentColor}(pixel::C, args...)
    base_colorant_type(C)(_clahe_pixel_rescale(color(pixel), args...), alpha(pixel))
end

"""
```
clipped_hist = cliphist(hist, clip)
```

Clips the histogram above a certain value `clip`. The excess left in the bins
exceeding `clip` is redistributed among the remaining bins.
"""
function cliphist{T}(hist::AbstractArray{T, 1}, clip::Number)
    hist_length = length(hist)
    excess = sum(map(i -> i > clip ? i - clip : zero(i - clip), hist))
    increase = excess / hist_length
    clipped_hist = zeros(Float64, hist_length)
    removed_sum = zero(Float64)

    for i in 1:hist_length
        if hist[i] > clip - increase
            clipped_hist[i] = Float64(clip)
            if hist[i] < clip
                removed_sum += hist[i] - (clip - increase)
            end
        else
            clipped_hist[i] = hist[i] + Float64(increase)
            removed_sum += increase
        end
    end

    leftover = excess - removed_sum

    while true
        oleftover = leftover
        for i in 1:hist_length
            leftover <= 0 && break
            step = ceil(Int, max(1, hist_length / leftover))
            for h in i:step:hist_length
                leftover <= 0 && break
                diff = clip - clipped_hist[h]
                if diff > 1
                    clipped_hist[h] += 1
                    leftover -= 1
                elseif diff > 0
                    clipped_hist[h] = clip
                    leftover -= diff
                end
            end
        end
        (leftover <= 0 || leftover >= oleftover) && break
    end
    clipped_hist
end


text:
"""
    LeafModels(interception, energy, photosynthesis, stomatal_conductance, status)
    LeafModels(;interception = missing, energy = missing, photosynthesis = missing,
        stomatal_conductance = missing,status...)

LeafModels component, which is a subtype of `AbstractComponentModel` implementing a component with
a photosynthetic activity. It could be a leaf, or a leaflet, or whatever kind of component
that is photosynthetic. The name `LeafModels` was chosen not because it is generic, but because it
is short, simple and self-explanatory.

# Arguments

- `interception <: Union{Missing,AbstractInterceptionModel}`: An interception model.
- `energy <: Union{Missing,AbstractEnergyModel}`: An energy model, e.g. [`Monteith`](@ref).
- `photosynthesis <: Union{Missing,AbstractAModel}`: A photosynthesis model, e.g. [`Fvcb`](@ref)
- `stomatal_conductance <: Union{Missing,AbstractGsModel}`: A stomatal conductance model, e.g. [`Medlyn`](@ref) or
[`ConstantGs`](@ref)
- `status <: MutableNamedTuple`: a mutable named tuple to track the status (*i.e.* the variables) of
the leaf. Values are set to `0.0` if not provided as VarArgs (see examples)

# Details

The status field depends on the input models. You can get the variables needed by a model using
[`variables`](@ref) on the instantiation of a model. Generally the variables are:

## Light interception model  (see [`AbstractInterceptionModel`](@ref))

- `Rₛ` (W m-2): net shortwave radiation (PAR + NIR). Often computed from a light interception model
- `PPFD` (μmol m-2 s-1): absorbed Photosynthetic Photon Flux Density
- `skyFraction` (0-2): view factor between the object and the sky for both faces.

## Energy balance model (see [`energy_balance`](@ref))

- `Tₗ` (°C): temperature of the object. Often computed from an energy balance model.
- `Rₗₗ` (W m-2): net longwave radiation for the object (TIR)
- `Cₛ` (ppm): stomatal CO₂ concentration
- `ψₗ` (kPa): leaf water potential
- `H` (W m-2): sensible heat flux
- `λE` (W m-2): latent heat flux
- `Dₗ` (kPa): vapour pressure difference between the surface and the saturation vapour

## Photosynthesis model (see [`photosynthesis`](@ref))

- `A` (μmol m-2 s-1): carbon assimilation
- `Gbₕ` (m s-1): boundary conductance for heat (free + forced convection)
- `Cᵢ` (ppm): intercellular CO₂ concentration

# Examples

```julia
# A leaf with a width of 0.03 m, that uses the Monteith and Unsworth (2013) model for energy
# balance, The Farquhar et al. (1980) model for photosynthesis, and a constant stomatal
# conductance for CO₂ of 0.0011 with no residual conductance. The status of
# the leaf is not set yet, all are initialised at `0.0`:
LeafModels(energy = Monteith(),
     photosynthesis = Fvcb(),
     stomatal_conductance = ConstantGs(0.0, 0.0011))

# If we need to initialise some variables at different values, we can call the leaf as:

LeafModels(photosynthesis = Fvcb(),Cᵢ = 380.0)

# Or again:
LeafModels(photosynthesis = Fvcb(), energy = Monteith(), Cᵢ = 380.0, Tₗ = 20.0)
```
"""
struct LeafModels{I <: Union{Missing,AbstractInterceptionModel},E <: Union{Missing,AbstractEnergyModel},A <: Union{Missing,AbstractAModel},Gs <: Union{Missing,AbstractGsModel},S <: MutableNamedTuple} <: AbstractComponentModel
    interception::I
    energy::E
    photosynthesis::A
    stomatal_conductance::Gs
    status::S
end

function LeafModels(;interception = missing, energy = missing,
                photosynthesis = missing, stomatal_conductance = missing,status...)
    status = init_variables_manual(interception, energy, photosynthesis,
        stomatal_conductance;status...)
    LeafModels(interception, energy, photosynthesis, stomatal_conductance, status)
end

"""
    Base.copy(l::T)

Copy a [`LeafModels`](@ref)
"""
function Base.copy(l::T) where T <: LeafModels
    LeafModels(l.interception, l.energy, l.photosynthesis, l.stomatal_conductance, MutableNamedTuple(; zip(keys(l.status), values(l.status))...))
end

"""
    Component(interception, energy, status)
    Component(;interception = missing, energy = missing, status...)

Generic component, which is a subtype of `AbstractComponentModel` implementing a component with
an interception model and an energy balance model. It can be anything such as a trunk, a
solar panel or else.

# Arguments

- `interception <: Union{Missing,AbstractInterceptionModel}`: An interception model.
- `energy <: Union{Missing,AbstractEnergyModel}`: An energy model.
- `status <: MutableNamedTuple`: a mutable named tuple to track the status (*i.e.* the variables) of
the component. Values are set to `0.0` if not provided as VarArgs (see examples)

# Examples

```julia
# An internode in a plant:
Component(energy = Monteith())
```
"""
struct Component{I <: Union{Missing,AbstractInterceptionModel},E <: Union{Missing,AbstractEnergyModel},S <: MutableNamedTuple} <: AbstractComponentModel
    interception::I
    energy::E
    status::S
end

function Component(;interception = missing, energy = missing,status...)
    status = init_variables_manual(interception, energy;status...)
    Component(interception, energy, status)
end


text:
module ImportanceSampling

using Distributions
using LinearAlgebra: norm

function pairwise_dist(Z)
  K = size(Z, 2)
  return [norm(Z[:, k] - Z[:, h], 1) for k in 2:K for h in 1:(k - 1)]
end

function rep_fn(Z, phi)
  d = pairwise_dist(Z)
  if phi == 0
    return 1.0
  else
    log_rep_term = phi * sum(log1p.(-exp.(-d)))
    rep_term = exp(log_rep_term)
    return rep_term
  end
end

"""
Approximate the expectation of a scalar function of Z ~ rFAM E[f(Z)|phi] using
importance sampling.
"""
function exp_f_rfam(; N::Integer, J::Integer, K::Integer,
                    repulsive_fn::Function, f::Function)
  # NOTE: In the importance sampling, the density of the regular
  # IBP component is cancelled, leaving the repulsive component.
  # Therefore, we only need binary random matrices to approximate
  # mean functions of Z ~ rfam.
  Zs = [rand(Bool, J, K) for _ in 1:N]

  weights = repulsive_fn.(Zs)
  numer = sum(f.(Zs) .* weights)
  denom = sum(weights)

  if denom == 0
    println("WARNING: weights are all 0!")
    return 0.0
  else
    return numer / denom
  end
end

min_pairwise_dist_ge_u(Z, u) = minimum(pairwise_dist(Z)) .>= u

end  # ImportanceSampling


text:
using DrWatson, Test
using BSON, DataFrames, FileIO, JLD2

@testset "Collect Results ($ending)" for ending ∈ ["bson", "jld2"]

###############################################################################
#                        Setup Folder structure                               #
###############################################################################
# %%
cd(@__DIR__)
isdir("testdir") && rm("testdir", recursive=true)
mkdir("testdir")
initialize_project("testdir"; git = false)
quickactivate("testdir")

###############################################################################
#                           Create Dummy Data                                 #
###############################################################################
mkdir(datadir("results"))
cd(datadir("results"))

d = Dict("a" => 1, "b" => "2", "c" => rand(10))
DrWatson.wsave(savename(d)*"."*ending, d)

d = Dict("a" => 3, "b" => "4", "c" => rand(10), "d" => Float64)
DrWatson.wsave(savename(d)*"."*ending, d)

d = Dict("a" => 3, "c" => rand(10), "d" => Float64)
DrWatson.wsave(savename(d)*"."*ending, d)

mkdir("subfolder")
cd("subfolder")

d = Dict("a" => 4., "b" => "twenty" , "d" => Int)
DrWatson.wsave(savename(d)*"."*ending, d)

###############################################################################
#                           Collect Data Into DataFrame                       #
###############################################################################
using Statistics
special_list = [ :lv_mean => data -> mean(data["c"]),
                :lv_var  => data -> var(data["c"])]

black_list = ["c"]

folder = datadir("results")

defaultname = joinpath(dirname(folder), "results_$(basename(folder))."*ending)
isfile(defaultname) && rm(defaultname)
cres = collect_results!(defaultname, folder;
    subfolders = true, special_list=special_list, black_list = black_list)

@test size(cres) == (4, 6)
for n in ("a", "b", "lv_mean")
    @test n ∈ String.(names(cres))
end
@test "c" ∉ names(cres)
@test all(startswith.(cres[!,"path"], projectdir()))

relpathname = joinpath(dirname(folder), "results_relpath_$(basename(folder))."*ending)
cres_relpath = collect_results!(relpathname, folder;
    subfolders = true, special_list=special_list, black_list = black_list,
    rpath = projectdir())
@info all(startswith.(cres[!,"path"], "data"))

###############################################################################
#                           Add another file in a sub sub folder              #
###############################################################################

@test isfile(defaultname)

mkdir("subsubfolder")
cd("subsubfolder")
d = Dict("b" => 35., "d" => Number, "c" => rand(5), "e" => "new_column")
DrWatson.wsave(savename(d)*".bson", d)
DrWatson.wsave(savename(d)*".jld2", d)

cres2 = collect_results!(defaultname, folder;
    subfolders = true, special_list=special_list, black_list = black_list)

@test size(cres2) == (6, 7)
@test all(names(cres) .∈ Ref(names(cres2)))

###############################################################################
#               Test additional syntax for special list                       #
###############################################################################

special_list2 = [ :lv_mean => data -> mean(data["c"]),
                data -> :lv_var  =>  var(data["c"]),
                data -> [:lv_mean2 => mean(data["c"]),
                        :lv_var2  =>  var(data["c"])]]
black_list = ["c"]

folder = datadir("results")
defaultname2 = joinpath(dirname(folder), "results_betterspeciallist."*ending)
isfile(defaultname2) && rm(defaultname2)
cres10 = collect_results!(defaultname2, folder;
    subfolders = true, special_list=special_list2, black_list = black_list)

@test size(cres10) == (6, 9)
for n in ("a", "b", "lv_mean", "lv_var", "lv_mean2", "lv_var2")
    @test n ∈ String.(names(cres10))
end
@test "c" ∉ names(cres10)
###############################################################################
#                           Load and analyze  DataFrame                       #
###############################################################################

df = load(defaultname)["df"]
@test size(df) == size(cres2)
@test sort(names(df)) == sort(names(cres2))

###############################################################################
#                            test empty whitelist                             #
###############################################################################

rm(defaultname)
cres_empty = collect_results!(defaultname, folder;
    subfolders = true, special_list=special_list, white_list=[])

@test dropmissing(cres2[!,[:lv_mean, :lv_var, :path]]) == dropmissing(cres_empty)

###############################################################################
#                           test out-of-place form                            #
###############################################################################
cd(@__DIR__)

cres3 = collect_results(folder;
subfolders = true, special_list=special_list, black_list = black_list)

@test sort(names(cres3)) == sort(names(cres2))
@test size(cres3) == size(cres2)

###############################################################################
#                              Quickactivate macro                            #
###############################################################################

cd(@__DIR__)
isdir("testdir") && rm("testdir", recursive=true)
initialize_project("testdir"; git = false)
open(joinpath("testdir","testinclude.jl"),"w") do f
    write(f,"@quickactivate\n")
end
include(joinpath("testdir", "testinclude.jl"))
@test Base.active_project() == abspath(joinpath("testdir","Project.toml"))

###############################################################################
#                                 Delete Folders                              #
###############################################################################

cd(@__DIR__)
rm("testdir", recursive=true)

end


text:
# Translation of code from Yuliy Sannikov, payoff_policy_growth, with the modification that instead of S, we plug in S²
function upwind_parabolic_pde(X, R, μ, Σ², G, V, dt_div_1pdt)
    N = length(X)
    dX = diff(X)

    # Perform upwind scheme w/centered difference on diffusion term
    Σ²0 = zeros(N)
    Σ²0[2:N-1] .= Σ²[2:N-1] ./ (dX[1:N-2] + dX[2:N-1]) # approx Σ² / (2 * dx): this term is the Σ²/2 coefficient
    DU = -(max.(μ[1:N-1], 0.) + Σ²0[1:N-1]) ./ dX .* dt_div_1pdt # up diagonal, μ divided by dX, Σ²0 is Σ² / (2 * dx^2), th
    DD = -(max.(-μ[2:N], 0.) + Σ²0[2:N]) ./ dX .* dt_div_1pdt # down diagonal, note should be negative b/c FD scheme makes DD negative, multiplied by negative drift ⇒ positive, then subtracted ⇒ negative

    # observe: Σ² and μ are zero at endpoints, hence Σ²0 zero at endpts too ->
    # boundary conditions for our PDE

    D0 = (1 - dt_div_1pdt) .* ones(N) + dt_div_1pdt .* R # diagonal
    D0[1:N-1] = D0[1:N-1] - DU
    D0[2:N] = D0[2:N] - DD # subtract twice b/c centered diff
    A = spdiagm(0 => D0, 1 => DU, -1 => DD) # + spdiagm(DU,1,N,N) + spdiagm(DD[1:N-1],-1,N,N)
    F = A \ (G .* dt_div_1pdt + V .* (1 - dt_div_1pdt)) # solve linear system

#=    # Equivalent to this code, which constructs the first and second finite difference matrices separately
    DU1 = zeros(N)
    DU2 = zeros(N)
    DD1 = zeros(N)
    DD2 = zeros(N)
    D01 = zeros(N)
    D02 = zeros(N)    DU1[2:N] = max.(μ[1:N - 1], 0.) ./ dX
    DU2[2:N] = Σ²0[1:N - 1] ./ dX
    DD1[1:N - 1] = max.(-μ[2:N], 0.) ./ dX
    DD2[1:N - 1] = Σ²0[2:N] ./ dX
    D01[1:N - 1] .= -DU1[2:N]
    D01[2:N] .-= DD1[1:N - 1]
    D02[1:N - 1] .= -DU2[2:N]
    D02[2:N] .-= DD2[1:N - 1]
    L1 = spdiagm(0 => D01, -1 => DD1[1:N - 1], 1 => DU1[2:N])
    L2 = spdiagm(0 => D02, -1 => DD2[1:N - 1], 1 => DU2[2:N])
    Acheck = spdiagm(0 => (1 - dt_div_1pdt) .* ones(N) + dt_div_1pdt .* R) - dt_div_1pdt .* (L1 + L2)=#

    return F
end


text:
using CLArrays, GPUArrays
using FileIO, Interpolations, Colors, ColorVectorSpace, FixedPointNumbers

function poincare_inner{N}(rv, result, c, π, ::Val{N}, n)
    # find next spiking neuron
    ϕ₁, ϕ₂, ϕ₃ = rv[1], rv[2], rv[3]
    πh = π / 2f0
    π2 = π * 2f0
    for unused = 1:N
        if ϕ₁ > ϕ₂
            if ϕ₁ > ϕ₃
                # first neuron is spiking
                dt = πh - ϕ₁
                # evolve phases till next spike time
                ϕ₁ = -πh
                ϕ₂ = atan(tan(ϕ₂ + dt) - c)
                ϕ₃ += dt
                # save state of neuron 2 and 3
                x = Cuint(max(round(((ϕ₂ + πh) / π) * (Float32(n) - 1f0)) + 1f0, 1f0))
                y = Cuint(max(round(((ϕ₃ + πh) / π) * (Float32(n) - 1f0)) + 1f0, 1f0))
                i1d = GPUArrays.gpu_sub2ind((n, n), (x, y))
                @inbounds if i1d <= Cuint(n * n) && i1d > Cuint(0)
                    accum = result[i1d]
                    # this is unsafe, since it could read + write from different threads, but good enough for the stochastic kind of process we're doing
                    result[i1d] = accum + 1f0
                end
                continue
            end
        else
            if ϕ₂ > ϕ₃
                # second neuron is spiking
                dt = πh - ϕ₂
                # evolve phases till next spike time
                ϕ₁ += dt
                ϕ₂ = -πh
                ϕ₃ = atan(tan(ϕ₃ + dt) - c)
                continue
            end
        end
        # third neuron is spikinga
        dt = πh - ϕ₃
        # evolve phases till next spike time
        ϕ₁ += dt
        ϕ₂ = atan(tan(ϕ₂ + dt) - c)
        ϕ₃ = -πh
    end
    return
end

function poincare_inner(n, seeds::GPUArray, result, c, π, val::Val{N}) where N
    foreach(poincare_inner, seeds, result, c, Float32(pi), val, n)
end

c = 1f0; divisor = 2^10
srand(2)
N = 10^10
ND = Cuint(2048)
AT = CLArray
result = AT(zeros(Float32, ND, ND))
_n = div(N, divisor)
jl_seeds = [ntuple(i-> rand(Float32), Val{3}) for x in 1:divisor]
seeds = AT(jl_seeds)
poincare_inner(ND, seeds, Base.RefValue(result), c, Float32(pi), Val{_n}())

cmap = interpolate(([
    RGB(0.0, 0.0, 0),
    RGB(0.2, 0.2, 0.9),
    RGB(0.2, 0.6, 0.9),
    RGB(0.7, 0.7, 0.98),
    RGB(0.8, 0.8, 0.9),
    RGB(0.82, 0.8, 1.0)
]), BSpline(Linear()), OnCell())

cn = length(cmap)
resultcpu = Array(result)

img_color = map(resultcpu) do val
    val = log(1.0 + log(1.0 + val))
    val = clamp(1.0 - val, 0f0, 1f0);
    idx = (val * (cn - 1)) + 1.0
    RGB{N0f8}(cmap[idx])
end
#save as an image
save(joinpath(@__DIR__, "poincare.png"), img_color)