text:
if !@isdefined(TestIGraph)
    struct TestIGraph <: AGraph  end
    struct TestIDiGraph <: ADiGraph end
    struct TestIEdge <: AEdge end
end

@testset "$TEST $G" begin

g = TestIGraph()
h = TestIDiGraph()
e = TestIEdge()

@test_throws ErrorException src(e)
@test_throws ErrorException dst(e)
@test_throws ErrorException reverse(e)

@test_throws ErrorException nv(g)
@test_throws ErrorException nv(h)
@test_throws ErrorException neighbors(g, 1)
@test edge(g, 1, 2) == Edge{Int}(1, 2)
@test edge(h, 1, 2) == Edge{Int}(1, 2)
@test_throws MethodError idx(Edge(1, 2))
@test_throws ErrorException neighbors(g, 1)
@test_throws ErrorException in_neighbors(h, 1)
@test_throws ErrorException out_neighbors(h, 1)
@test_throws ErrorException graphtype(h)
@test_throws ErrorException digraphtype(g)
@test edgetype(g) == Edge{Int}
@test_throws ErrorException  pop_vertex!(g)
@test vertextype(g) == Int
@test_throws ErrorException add_edge!(g, 1, 2)
@test_throws ErrorException ne(g)
@test_throws ErrorException rem_edge!(g, 1, 2)
@test_throws ErrorException add_vertex!(g)
@test_throws ErrorException swap_vertices!(g, 1, 2)

@test graphtype(Graph) == Graph{Int}
@test graphtype(DiGraph) == Graph{Int}
@test digraphtype(Graph) == DiGraph{Int}
@test digraphtype(DiGraph) == DiGraph{Int}

end # testset


text:
using Documenter, StanVariational

makedocs(
    modules = [StanVariational],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "StanVariational.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/goedman/StanVariational.jl.git",
)

ext: jl


text:
module CollectiveVariables

# Signal to Revise.jl that this file should be tracked as a package
__revise_mode__ = :eval

# Dependecies
using ChainRulesCore
using CUDA
using Parameters
using StaticArrays
using Zygote

# Imports
using LinearAlgebra: ×, ⋅, norm
using Zygote: _pullback, sensitivity, tailmemaybe

# Exports
export DihedralAngle, DynamicCV, value_and_gradient

# Types
abstract type CVCache <: Function end

struct NoCache <: CVCache
end

struct ValGradCache{V, G, E} <: CVCache
    val::V
    grad::G
    extgrad::E

    function ValGradCache(cv, rs::AbstractArray{T}) where {T}
        val, V = withtype(value_storage(rs))
        grad, G = withtype(grad_storage(cv, rs))
        extgrad, E = withtype(reinterpret(reshape, SVector{3, T}, zero(rs)))
        return new{V, G, E}(val, grad, extgrad)
    end
end

abstract type CollectiveVariable <: Function end

abstract type GroupCV      <: CollectiveVariable end
abstract type FourPointCV  <: CollectiveVariable end
abstract type ThreePointCV <: CollectiveVariable end

struct DynamicCV{CV, T, C <: CVCache} <: CollectiveVariable
    ξ::CV
    indices::T
    cache::C

    function DynamicCV(cv::CV, rs, indices::T) where {CV, T}
        cache, C = withtype(ValGradCache(cv, rs))
        return new{CV, T, C}(cv, indices, cache)
    end

    function DynamicCV(cv::CV, rs, indices::T, cache::NoCache) where {CV, T}
        return new{CV, T, NoCache}(cv, indices, cache)
    end
end

struct DihedralAngle <: FourPointCV
    indices::SVector{4, Int}
end

DihedralAngle(i, j, k, l) = DihedralAngle(SVector(i, j, k, l))

# Methods
cv_function(::DihedralAngle) = dihedral_angle
cv_function(cv::DynamicCV) = cv_function(cv.ξ)

get_cache(cv::DynamicCV) = cv.cache

get_indices(cv::CollectiveVariable) = cv.indices

get_tags(cv::DihedralAngle) = cv.indices
get_tags(cv::DynamicCV) = get_tags(cv.ξ)

"""    dihedral_angle(p₁, p₂, p₃, p₄)

Computes the dihedral (or torsion) angle defined by four points in space (around the line
defined by the two central points).
"""
@inline dihedral_angle(p₁, p₂, p₃, p₄) = dihedral_angle(p₂ - p₁, p₃ - p₂, p₄ - p₃)
#
@inline function dihedral_angle(a, b, c)
    p = a × b
    q = b × c
    return atan((p × q) ⋅ b, (p ⋅ q) * norm(b))
end

# Gradient evaluation
get_value(cache::CVCache) = cv.val
get_gradient(cache::CVCache) = cv.grad
get_extended_gradient(cache::CVCache) = cv.extgrad

value_storage(rs::AbstractArray{T}) where {T} = Ref(zero(T))
value_storage(rs::AnyCuArray{T}) where {T} = CUDA.fill(zero(T))

grad_storage(cv::FourPointCV, rs::AbstractArray{T}) where {T} = nothing
#
function grad_storage(cv::FourPointCV, rs::AnyCuArray{T}) where {T}
    return CUDA.fill(zeros(SVector{3, T}), 4)
end

@inline function withstaticgrad(f::F, ps) where {F}
    ξ, pb = _pullback(Zygote.Context(), f, ps...)
    ∂ξ = tailmemaybe(pb(sensitivity(ξ)))
    return ξ, SVector(∂ξ)
end

function (cv::DynamicCV)(rs::AbstractMatrix{T}) where {T <: Real}
    I = get_indices(cv)
    tags = get_tags(cv)
    ps = reinterpret(reshape, SVector{3, T}, rs)
    f = ps -> withstaticgrad(cv_function(cv), ps)
    ξ, ∂ξ = evaluate!(f, get_cache(cv), ps, I, tags)
    return ξ, ∂ξ
end

function evaluate!(f::F, cache::ValGradCache, ps, I::AbstractVector, tags) where {F}
    ξ = get_value(cache)
    ∂p = get_extended_gradient(cache)
    fill!(∂p, zero(eltype(∂p)))

    inds = SVector(view(I, tags) .+ 1)
    ξ[], ∂ξ = f(view(ps, inds))
    ∂p[inds] .= ∂ξ

    return ξ[], vec(reinterpret(reshape, eltype(ξ), ∂p))
end

@inline function evaluate!(f::F, cache::ValGradCache, ps, I::AnyCuVector, tags) where {F}
    ξ = get_value(cache)
    ∂ξ = get_gradient(cache)
    ∂p = get_extended_gradient(cache)
    ∂x = reinterpret(reshape, SVector{4, eltype(ps)}, ∂ξ)
    fill!(∂p, zero(eltype(∂p)))

    function kernel(ξ, ∂ξ, ∂x, ∂p, ps, I)
        i = threadIdx().x
        @inbounds if i ≤ length(∂ξ)
            ∂ξ[i] = ps[I[tags[i]] + 1]
        end
        sync_threads()
        @inbounds if i == length(∂ξ) + 1
            ξ[], ∂x[] = f(∂x[])
        end
        sync_threads()
        @inbounds if i ≤ length(∂ξ)
            ∂p[I[tags[i]] + 1] = ∂ξ[i]
        end
        return nothing
    end

    @cuda threads=(length(∂ξ) + 1) kernel(ξ, ∂ξ, ∂x, ∂p, ps, I)

    return CUDA.@allowscalar(ξ[]), vec(reinterpret(reshape, eltype(ξ), ∂p))
end

function value_and_gradient(cv::CollectiveVariable, rs::AbstractMatrix{T}) where {T <: Real}
    ids = indices(cv)
    ps = reinterpret(reshape, SVector{3, T}, view(parent(rs), 1:3, ids))
    # TODO: Preallocate the memory for ∂ξ outside the function
    ∂ξ = reinterpret(reshape, SVector{3, T}, zero(rs))
    ξ, pb = _pullback(cv_function(cv), ps...)
    ∂p = tailmemaybe(pb(sensitivity(ξ)))
    if ∂p !== nothing
        ∂ξ[ids] .= ∂p
    end
    return ξ, vec(reinterpret(reshape, T, ∂ξ))
end

# Some (type-pirating ☠) bug fixes to ChainRulesCore

function ChainRulesCore.ProjectTo(x::AbstractArray{T}) where {T <: AbstractFloat}
    return ProjectTo{AbstractArray}(; element=ProjectTo(zero(T)), axes=axes(x))
end

function (project::ProjectTo{AbstractArray})(dx::AbstractArray{S}) where {S <: Number}
    T = ChainRulesCore.project_type(project.element)
    return S <: T ? dx : map(project.element, dx)
end

#function ChainRulesCore.rrule(
#    ::typeof(reinterpret), ::typeof(reshape), ::Type{R}, A::AbstractMatrix{T}
#) where {N, T, R <: SVector{N, T}}
#    y = reinterpret(reshape, R, A)
#    @inline function pullback(V̄)
#        ∂A = reinterpret(reshape, R, zero(A))
#        for (i, v) in pairs(V̄)
#            if v isa AbstractVector
#                ∂A[i] = v
#            end
#        end
#        #∂A = mapreduce(v -> v isa R ? v : zero(R), vcat, Ā; init = similar(A, 0))
#        return (NO_FIELDS, NO_FIELDS, DoesNotExist(), reinterpret(reshape, T, ∂A))
#    end
#    return y, pullback
#end

# Other utils

withtype(x::T) where {T} = (x, T)

end  # module CollectiveVariables


text:
using stoosc
rc("font", family="arial")

function params_prc_create_model!(params)

    println("Simulating reference ...")

    # Get input parameters
    parameter_names = params["parameter_names"]
    parameter_values = params["parameter_values"]
    trajectories = params["trajectories"]
    T = 1.0  # the model should have free running period 1
    n_pulses = params["n_pulses"]
    show_plots = params["show_plots"]
    model_name = params["model_name"]
    problem_type = params["problem_type"]
    i_variable = params["i_variable"]

    tend = 45.0  # final simulation time
    toffset = 25.0  # offset to remove the effect of initial conditions

    # Create a model
    model = load_model(model_name, problem_type)
    model = set_timespan(model, (0.0, tend))
    model = set_solver(model; saveat=0.01, maxiters=1e10)
    model = set_parameter_value(model, parameter_names, parameter_values)
    outfun = sol -> hcat(sol[i_variable, :])
    model = set_output(model, outfun)

    # Generate LD cycle for the whole tspan
    events = create_events([(:LD, tend/T + 1, T/2, T/2)])

    # Keep only 10 LD cycles after timepoint toffset
    idx = events[:, 1] .> toffset
    idx_first_event = findall(idx)[1]  # first event after toffset
    time_first_events = events[idx_first_event, 1]
    idx_last_event = idx_first_event + 9
    events = events[1:idx_last_event, :]  # remove tailing events to create DD

    # Add events to the model and run the simulation
    model = set_input(model, events)
    solution = simulate_population(model, trajectories, save_trajectories=false,
        show_progress=true)

    if show_plots
        _, ax = subplots(figsize=(6, 3))
        plot_solution(solution, ax=ax)
    end

    # Select solution from the first event after time toffset
    solution_selection = select_time(solution, mintime=time_first_events)
    events_selection = solution_selection.events

    # Find peaks in the selected solution
    t = solution_selection.time
    x = solution_selection.mean[:, 1]
    pr = findpeaks(x, t, npeaks=15, sortstr="descend", sortref="prominence")
    pr = sort(pr, ref="location")
    pks = peakheights(pr)
    locs = peaklocations(pr)

    # Calculate entrainment phase, free-running period, and pulse starts
    C0 = mean(locs[1:10] .- events_selection[:, 1])  # entrainment phase
    FRP = locs[13] - locs[12]  # free-running period
    # reference_location = locs[15]  # reference peak
    pulse_starts = range(locs[12] - C0, locs[13] - C0, length=n_pulses)
    reference_time = pulse_starts[end] + T

    idx = t .> reference_time
    xarr = [x[idx]]
    tarr = [t[idx]]

    # Plot the calculated properties for visual verification
    if show_plots
        _, ax = subplots(figsize=(6, 3))
        ax.plot(t, x, color="black")
        ax.plot(locs[1:10], pks[1:10], "o", color="blue")
        ax.plot(locs[12:13], pks[12:13], "o", color="red")
        for pulse_start in pulse_starts
            ax.vlines(pulse_start, minimum(x), maximum(x), color="blue")
        end
        ax.plot(tarr[1], xarr[1], color="green")
        ax.vlines(reference_time, minimum(x), maximum(x), color="green")
        plot_events(events_selection, ax=ax)
    end

    # Save important variables
    params["model"] = model
    params["events"] = events
    params["time_first_events"] = time_first_events
    params["entrainment_phase"] = C0
    params["free_running_period"] = FRP
    params["reference_time"] = reference_time
    params["pulse_starts"] = time_first_events .+ pulse_starts
    params["xarr"] = xarr
    params["tarr"] = tarr

end

function params_prc_estimate_pulse_responses!(params)

    println("Simulating pulse responses ...")

    T = 1.0  # params["ode_free_running_period"]
    pulse_legth = T / 2
    pulse_starts = params["pulse_starts"]
    events = params["events"]
    model = params["model"]
    time_first_events = params["time_first_events"]
    xarr = params["xarr"]
    tarr = params["tarr"]
    reference_time = params["reference_time"]
    trajectories = params["trajectories"]
    n_pulses = params["n_pulses"]
    show_plots = params["show_plots"]

    if show_plots
        _, axarr = subplots(n_pulses, figsize=(6, 4))
    end

    progressmeter = ProgressMeter.Progress(n_pulses; barlen=20)
    for i in 1:n_pulses

        pulse_start = pulse_starts[i]
        PRC_events = vcat(events, [pulse_start pulse_start + pulse_legth])
        PRC_model = set_input(model, PRC_events)

        success = false
        counter = 0
        solution = stoosc.PopulationSolution()
        while !success
            try
                solution = simulate_population(PRC_model, trajectories,
                    save_trajectories=false, show_progress=false)
                success = true
            catch err
                counter += 1
                if counter == 5
                    rethrow(err)
                else
                    @warn "[$counter] Simulation threw an error: $err"
                end
            end

        end

        # Select solution from the first event after time toffset
        solution_selection = select_time(solution, mintime=time_first_events)
        events_selection = solution_selection.events

        # Find peaks in the selected solution
        t = solution_selection.time
        x = solution_selection.mean[:, 1]

        idx = t .> reference_time
        push!(xarr, x[idx])
        push!(tarr, t[idx])

        ProgressMeter.next!(progressmeter)

        # Plot the calculated properties for visual verification
        if show_plots
            axarr[i].plot(t, x, color="black")
            plot_events(events_selection, ax=axarr[i])
            axarr[i].plot(tarr[end], xarr[end], color="green")
            axarr[i].axis("off")
        end

    end

end

function params_prc_estimate_prc!(params)

    show_plots = params["show_plots"]
    tarr = params["tarr"]
    xarr = params["xarr"]
    free_running_period = 1.0
    pulse_starts = (params["pulse_starts"] .- params["pulse_starts"][1]) ./ free_running_period
    PRC = Float64[]

    n = length(tarr)
    tref = tarr[1]
    xref = xarr[1]
    Ts = mean(diff(tref))
    for i = 2:n

        R = crosscor(xarr[i], xref, ((1-length(tref)):(length(tref)-1)))
        lags = ((1-length(tref)):(length(tref)-1)).*Ts
        pr = findpeaks(R, lags)
        locs = peaklocations(pr)
        pks = peakheights(pr)
        idx = argmin(abs.(locs))
        pk = pks[idx]
        loc = locs[idx]
        push!(PRC, loc ./ free_running_period)

    end

    params["PRC"] = DataFrame([pulse_starts, PRC],
        ["time", "phase_shift"])

    if show_plots
        fig, ax = subplots()
        df = params["PRC"]
        ax.scatter(df[:, "time"], df[:, "phase_shift"])
    end

end

function params_prc_save_prc(params)

    PRC = params["PRC"]
    filename = params["filename"]
    save_data(PRC, filename)

end

params = Dict{String, Any}(
    "parameter_names" => ["d", "B", "τ", "I", "σ"],  # parameter names
    "parameter_values" => [-0.1, 1, 6.2, 3.5, 0.01],  # parameter values
    "trajectories" => 1000,  # number of oscillators in the population
    "n_pulses" => 30,  # number of pulses
    "model_name" => "model02",  # model01 or model02
    "problem_type" => "sde",  # ode or sde
    "i_variable" => 1,  # variable to estimate PRC
    "filename" => "test_prc.csv",  # save the PRC here
    "show_plots" => true  # visual control
)

params_prc_create_model!(params)
params_prc_estimate_pulse_responses!(params)
params_prc_estimate_prc!(params)
params_prc_save_prc(params)

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
# Autogenerated wrapper script for ranger_jll for x86_64-apple-darwin
export libranger

JLLWrappers.@generate_wrapper_header("ranger")
JLLWrappers.@declare_library_product(libranger, "@rpath/libranger.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libranger,
        "lib/libranger.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Gato

export Tablero, nvo_juego, tiro_valido, turno
import Base.show

# El tipo Tablero guarda la configuración de un juego

type Tablero
	juego_A::Array{Int}
	juego_B::Array{Int}
	representacion::Array{ASCIIString}
	turno::Bool
	cuenta_turnos::Int
end

function show(io::IO, t::Tablero)
	rep  = t.representacion
	for reng in 1:3
		print(io, rep[reng,1], "\t", rep[reng,2], "\t", rep[reng,3], "\n")
	end
end

# Las funciones

# Crea un juego vacío

function nvo_juego()
	turno = true
	cuenta_turnos = 0
	A = zeros(Int, 3, 3)
	B = zeros(Int, 3, 3)

	repr = Array(ASCIIString, (3, 3))
	fill!(repr, ".")

	Tablero(A, B, repr, turno, cuenta_turnos)
end

# Prueba si el arreglo que se le da como argumento es una configuración ganadora

function prueba_ganador(M::Array{Int})
	v = sum(M, 1)
	h = sum(M, 2)
	diagonal = diag(M)
	otra_diag = [M[1, 3], M[2, 2], M[3, 1]]

	if diagonal == [1, 1, 1] || otra_diag == [1, 1, 1]
		return true
	else
		for i in 1:3
			if v[i] == 3 || h[i] == 3
				return true
			end
		end
	end

	false
end

# Verifica que la posción que se escoge es una casilla vacía

function tiro_valido(t::Tablero, pos::Int)
	if t.juego_A[pos] == 0 && t.juego_B[pos] == 0
		true
	else
		false
	end
end

# Hace efectiva la tirada de un jugador

function tirada(t::Tablero, pos::Int)
	if t.turno
		t.juego_A[pos] = 1
		t.representacion[pos] = "X"
	else
		t.juego_B[pos] = 1
		t.representacion[pos] = "O"
	end
end

# Un turno del juego

function turno(t::Tablero, pos::Int)
	tirada(t, pos)

	if t.turno
		termina = prueba_ganador(t.juego_A)
		ganador = "X"
	else
		termina = prueba_ganador(t.juego_B)
		ganador = "O"
	end

	if termina
		return true, ganador
	else
		# Si no hay un ganador, es el turno del siguiente jugador
		t.turno = !t.turno
		return false, ""
	end
end

end # Termina módulo Gato

max_stars_repo_name: rodrigolece/juego-gato

max_issues_repo_name: rodrigolece/juego-gato

max_forks_repo_name: rodrigolece/juego-gato


text:
######### Stan batch program example  ###########

using StanSample, Test

dyes ="
data {
  int BATCHES;
  int SAMPLES;
  real y[BATCHES, SAMPLES];
  // vector[SAMPLES] y[BATCHES];
}

parameters {
  real<lower=0> tau_between;
  real<lower=0> tau_within;
  real theta;
  real mu[BATCHES];
}

transformed parameters {
  real sigma_between;
  real sigma_within;
  sigma_between <- 1/sqrt(tau_between);
  sigma_within <- 1/sqrt(tau_within);
}

model {
  theta ~ normal(0.0, 1E5);
  tau_between ~ gamma(.001, .001);
  tau_within ~ gamma(.001, .001);

  mu ~ normal(theta, sigma_between);
  for (n in 1:BATCHES)
    y[n] ~ normal(mu[n], sigma_within);
}

generated quantities {
  real sigmasq_between;
  real sigmasq_within;

  sigmasq_between <- 1 / tau_between;
  sigmasq_within <- 1 / tau_within;
}
"

dyesdata = Dict("BATCHES" => 6,
    "SAMPLES" => 5,
    "y" => reshape([
      [1545, 1540, 1595, 1445, 1595];
      [1520, 1440, 1555, 1550, 1440];
      [1630, 1455, 1440, 1490, 1605];
      [1595, 1515, 1450, 1520, 1560];
      [1510, 1465, 1635, 1480, 1580];
      [1495, 1560, 1545, 1625, 1445]
    ], 6, 5)
  )

sm = SampleModel("dyes", dyes);

rc = stan_sample(sm, data=dyesdata)

if success(rc)
  samples = read_samples(sm)

  # Fetch cmdstan summary data frame.
  df = read_summary(sm)
  @test df[df.parameters .== :theta, :mean][1] ≈ 1527.5 atol=50.0
  @test df[df.parameters .== Symbol("mu[1]"), :mean][1] ≈ 1512.5 rtol=0.1
  @test df[df.parameters .== Symbol("mu[2]"), :mean][1] ≈ 1528.1 rtol=0.1
  @test df[df.parameters .== Symbol("mu[3]"), :mean][1] ≈ 1553.8 rtol=0.1
  @test df[df.parameters .== Symbol("mu[4]"), :mean][1] ≈ 1507.2 rtol=0.1
  @test df[df.parameters .== Symbol("mu[5]"), :mean][1] ≈ 1578.7 rtol=0.1
  @test df[df.parameters .== Symbol("mu[6]"), :mean][1] ≈ 1487.3 rtol=0.1

end


text:
using ApproxFun, SingularIntegralEquations , Plots
gr()

k = 50.
ω = 2π
d = (1,-1)
d = d[1]/hypot(d[1],d[2]),d[2]/hypot(d[1],d[2])
ui = (x,y) -> exp(im*k*(d[1]*x+d[2]*y))

# The Helmholtz Green's function, split into singular and nonsingular pieces.
g1 = (x,y) -> -besselj0(k*abs(y-x))/2
g2 = (x,y) -> x == y ? -(log(k/2)+γ)/2/π + im/4 : im/4*hankelh1(0,k*abs(y-x)) - g1(x,y).*logabs(y-x)/π

Γ = Curve(Fun(x->x+im*x^3))

#Γ = Segment(0.,1.0+im)

sp = Space(Γ)
cwsp = CauchyWeight(sp⊗sp,0)
uiΓ,⨍ = Fun(t->ui(real(t),imag(t)),sp),DefiniteLineIntegral(Γ)

@time G = GreensFun(g1,cwsp;method=:Cholesky) + GreensFun(g2,sp⊗sp;method=:Cholesky)

@time ∂u∂n = ⨍[G]\uiΓ

N=100
x = linspace(-3,3,N);y = linspace(-2,2,N)

us = (x,y) -> -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

@time uvals = ui.(x,y') .+ us.(x,y')

contourf(x,y,real(uvals'))
    plot!(Γ)

println("The length of ∂u∂n is: ",ncoefficients(∂u∂n))
us = (x,y) -> -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

t=0.1;us(t,t^3)-ui(t,t^3)
t=0.2;us(t,t^3)-ui(t,t^3)

t=0.1;us(t,t)+ui(t,t)
t=0.2;us(t,t)+ui(t,t)

Γ = Curve(Fun(x->x+im*x^3)) ∪ Curve(Fun(x->x-1+im*(x^4+4)))

plot(Γ)

sp = Space(Γ)
cwsp = CauchyWeight(sp⊗sp,0)

uiΓ,⨍ = Fun(t->ui(real(t),imag(t)),sp),DefiniteLineIntegral(Γ)

@time G = GreensFun(g1,cwsp;method=:Cholesky) + GreensFun(g2,sp⊗sp;method=:Cholesky)

@time ∂u∂n = ⨍[G]\uiΓ
println("The length of ∂u∂n is: ",ncoefficients(∂u∂n))
us = (x,y) -> -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

N=300
x = linspace(-4,3,N);y = linspace(-2,6,N)

us = (x,y) -> -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

@time uvals = ui.(x,y') .+ us.(x,y')

writecsv("/Users/solver/Desktop/uvals.csv",uvals)

contourf(x,y,real(uvals');legend=false,size=(1000,600))
    plot!(Γ;color=:black)
png("scat.png")

3

@time us(0.1,0.2)

writecsv("dudncoeffs.csv",∂u∂n.coefficients)

# on 0.5
dsp = domainspace(DefiniteLineIntegral(Γ))
cfs = eval.(parse.(vec(readcsv("/Users/solver/Desktop/dudncoeffs.csv"))))

∂u∂n = Fun(dsp,chop!(cfs,1E-5))

us = (x,y) -> -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

us(0.1,0.2)

N=300
x = linspace(-4,3,N);y = linspace(-2,6,N)
xx,yy = x.+0.*y',0.*x.+y'

@time uvals = ui.(xx,yy) .+ us.(xx,yy)

using Plots
pyplot()

plot(Γ;color="black")

umax = maximum(abs,uvals)
contourf(x,y,real(uvals)')#,cmap="cubehelix")
    plot!(Γ;color="black")

xlabel("\$x\$");ylabel("\$y\$");colorbar(ax=gca(),shrink=0.685)#515#2/3)
savefig("ScatteringNeumannPDEplot.png";dpi=600,bbox_inches="tight")

function makegif(x,y,u,L;plotfunction=Main.PyPlot.contourf,seconds=1,cmap="seismic",vert=1)
    tm=string(time_ns())
    dr = pwd()*"/"*tm*"mov"
    mkdir(dr)

    umax = maxabs(u)
    fps = 24
    MLen = seconds*fps

    tt = linspace(-1.,1.,100)

    for k=1:MLen
        t = 2π/ω*(k-1)/fps
        Main.PyPlot.clf()
        Main.PyPlot.axes(aspect="equal")
        Main.PyPlot.plot(tt,tt.^3;color=:black)
        Main.PyPlot.plot(tt-1,tt.^4.+4;color=:black)
        plotfunction(x,y,real(u*exp(-im*ω*t))';vmin=-umax*vert,vmax=umax*vert,cmap="seismic")
        xlabel!("\$x\$");ylabel!("\$y\$")
        Main.PyPlot.savefig(dr * "/" * lpad(k,max(4,ceil(Int,log10(MLen))),0) * ".png";dpi=150,bbox_inches="tight")
    end
    # Requires: brew install imagemagick
    run(`convert -delay 6 -loop 0 $dr/*.png $dr/scattering.gif`)
    run(`open $dr/scattering.gif`)
end

makegif(x,y,uvals,100;seconds=1,cmap="seismic",vert=0.5)

fps = 24
seconds=1
MLen = seconds*fps

k=MLen
u=uvals;
x=xx;y=yy;
k=24
    t = 2π/ω*(k-1)/fps
    real(u*exp(-im*ω*t))
    Main.PyPlot.contourf(x,y,real(u*exp(-im*ω*t))')
    plot!(Γ;color=:black,legend=false)
xlabel!("\$x\$");ylabel!("\$y\$")

tm=string(time_ns())
dr = pwd()*"/"*tm*"mov"
mkdir(dr)

Main.PyPlot.savefig(dr * "/" * lpad(k,max(4,ceil(Int,log10(MLen))),0) * ".png";dpi=150,bbox_inches="tight")

max_stars_repo_path: examples/HelmHoltz on a curve.jl

max_issues_repo_path: examples/HelmHoltz on a curve.jl

max_forks_repo_path: examples/HelmHoltz on a curve.jl


text:
"""
This `OCIM2_48L` module is used to load the OCIM2-48L matrices and grid for use in AIBECS.

!!! tip
    To load the default OCIM2_48L matrix and grid, do
    ```
    julia> grd, T = OCIM2_48L.load()
    ```
    But you can also load the other matrices by specifying which version you want, e.g.,
    ```
    julia> grd, T = OCIM2_48L.load(version="OCIM2_48L_KiHIGH_noHe")
    ```
    See *DeVries and Holzer* (2019) for more details

!!! note
    The files, that are downloaded from a public and persistant URL in FigShare,
    were created with the code available at https://github.com/briochemc/OceanCirculations.
"""
module OCIM2_48L

using SparseArrays          # for sparse matrix
using DataDeps              # for automated download/storage of data
using Downloads
using MAT                   # for reading OCIM2_48L transport operator matrix from MAT file
using NCDatasets            # for reading OCIM2_48L grid from NetCDF file
using Unitful               # for units
using Unitful: s, yr, m, km, °
using Reexport
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Downloads.download(remotepath, localpath)
    return localpath
end

# OCIM2_48L URL
url() = "https://files.figshare.com/28468077/OCIM2_48L_base.tar.gz"

# OCIM2_48L Hashes
sha() = "7c6e7981df69122957c9f2475be7ffe958819ff44d02be6d004886b6dfc3fac7"

OCIM2versionerror(version) = error("""`$version` is not a valid OCIM2 version name.

                                   Valid versions are `CTL_He`, `CTL_noHe`, `KiHIGH_He`, `KiHIGH_noHe`, `KiLOW_He`, `KiLOW_noHe`, `KvHIGH_He`, `KvHIGH_KiHIGH_noHe`, `KvHIGH_KiLOW_He`, `KvHIGH_KiLOW_noHe`, and `KvHIGH_noHe`.

                                   See *DeVries and Holzer* (2019) for more details on OCIM2 configurations.""")

# Create registry entry for OCIM in JLD2 format
function register_OCIM2_48L()
    register(
        DataDep(
            "AIBECS-OCIM2_48L",
            """
            References:
            - $(citations())
            """,
            url(),
            sha(),
            fetch_method = fallback_download,
            post_fetch_method = unpack
        )
    )
    return nothing
end

"""
    load

Returns the grid and the transport matrix of the OCIM2_48L model.

See *DeVries and Holzer* (2019) and *Holzer et al.* (2021) for more details.
"""
function load()
    register_OCIM2_48L()
    files_path = @datadep_str "AIBECS-OCIM2_48L"
    @info """You are about to use the OCIM2_48L model.
          If you use it for research, please cite:

          $(citations())

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "DeVries_Holzer_2019" and "Holzer_etal_2021" keys.)
          """

    # Convert convergence T in yr⁻¹ from original MAT file to divergence T in s⁻¹
    T = -ustrip.(s^-1, matread(joinpath(files_path, "OCIM2_48L_base_transport.mat"))["TR"] * yr^-1);

    # Create grid object from NetCDF file variables
    # TODO: Add Hyrothermal He fluxes as for the OCIM2 load function
    grid = Dataset(joinpath(files_path, "OCIM2_48L_base_data.nc"), "r") do ds
        wet3D = ds["ocnmask"][:] .== 1
        nlat, nlon, ndepth = size(wet3D)
        lat_3D = ds["tlat"][:]°
        lon_3D = ds["tlon"][:]°
        depth_3D = ds["tz"][:]m
        lat = unique(lat_3D)
        lon = unique(lon_3D)
        depth = unique(depth_3D)
        ulat = unique(ds["ulon"][:])° # ulat↔ulon in OCIM2-48L original files
        ulon = unique(ds["ulat"][:])° # ulat↔ulon in OCIM2-48L original files
        depth_top_3D = ds["wz"][:]m
        depth_top = unique(depth_top_3D)
        δlat = 2(ulat - lat)
        δlon = 2(ulon - lon)
        δdepth = 2(depth - depth_top)
        R = 6371.0km
        δy = R * δlat ./ 360°
        δy_3D = repeat(reshape(δy, (nlat,1,1)), outer=(1,nlon,ndepth))
        A_3D = ds["area"][:]m^2
        δx_3D = A_3D ./ δy_3D
        volume_3D = ds["vol"]m^3
        δz_3D = volume_3D ./ A_3D
        A_2D = A_3D[:,:,1]
        nboxes = count(wet3D)
        OceanRectilinearGrid(
            lat,
            lon,
            depth,
            δlat,
            δlon,
            δdepth,
            lat_3D,
            lon_3D,
            depth_3D,
            δy,
            δx_3D,
            δy_3D,
            δz_3D,
            volume_3D,
            depth_top,
            depth_top_3D,
            A_2D,
            wet3D,
            nlon,
            nlat,
            ndepth,
            nboxes
        )
    end

    return grid, T

end

citations() = """
              - Holzer, M., DeVries, T. & de Lavergne, C. Diffusion controls the ventilation of a Pacific Shadow Zone above abyssal overturning. Nat Commun 12, 4348 (2021). https://doi.org/10.1038/s41467-021-24648-x
              - DeVries, T., & Holzer, M. (2019). Radiocarbon and helium isotope constraints on deep ocean ventilation and mantle-3He sources. Journal of Geophysical Research: Oceans, 124, 3036–3057. https://doi.org/10.1029/2018JC014716
              """

end # end module

export OCIM2_48L


text:
module compare_table_to_schema

export compare

using ..handle_validvalues
using ..types

include("compare/common.jl")
include("compare/inmemory_table.jl")
include("compare/ondisk_table.jl")

using .common
using .inmemory_table
using .ondisk_table

"""
Compares a table to a TableSchema and returns:
- A copy of the input table, transformed as much as possible according to the schema.
- A table of the ways in which the input table doesn't comply with the schema.
- A table of the ways in which the output table doesn't comply with the schema.

There are currently 2 methods for comparing a table to a schema:

1. `compare(table, tableschema)` compares an in-memory table.

2. `compare(tableschema, input_data_file::String; output_data_file="", input_issues_file="", output_issues_file="")` compares a table stored on disk in `input_data_file`.

   This method is designed for tables that are too big for RAM.
   It examines one row at a time.
   The 3 tables of results (see above) are stored on disk. By default they are stored in the same directory as the input table.
"""
compare(tableschema::TableSchema, table; sorted_by_primarykey::Bool=false) = inmemory_table.compare(tableschema, table, sorted_by_primarykey)

function compare(tableschema::TableSchema, input_data_file::String;
                 output_data_file::String="", input_issues_file::String="", output_issues_file::String="", sorted_by_primarykey::Bool=false)
    !isfile(input_data_file) && error("The input data file does not exist.")
    fname, ext = splitext(input_data_file)
    output_data_file   = output_data_file   == "" ? "$(fname)_transformed.tsv"   : output_data_file
    input_issues_file  = input_issues_file  == "" ? "$(fname)_input_issues.tsv"  : input_issues_file
    output_issues_file = output_issues_file == "" ? "$(fname)_output_issues.tsv" : output_issues_file
    outdir  = dirname(output_data_file)  # outdir = "" means output_data_file is in the pwd()
    outdir != "" && !isdir(outdir) && error("The directory containing the specified output file does not exist.")
    ondisk_table.compare(tableschema, input_data_file, output_data_file, input_issues_file, output_issues_file, sorted_by_primarykey)
end

end


text:
module MPI

using Libdl, Serialization

macro mpichk(expr)
    @assert expr isa Expr && expr.head == :call && expr.args[1] == :ccall
    if Sys.iswindows()
        insert!(expr.args, 3, :stdcall)
    end
    esc(expr)
end

function serialize(x)
    s = IOBuffer()
    Serialization.serialize(s, x)
    take!(s)
end

function deserialize(x)
    s = IOBuffer(x)
    Serialization.deserialize(s)
end

primitive type SentinelPtr
    Sys.WORD_SIZE
end
Base.cconvert(::Type{Ptr{T}}, sptr::SentinelPtr) where {T} = reinterpret(Ptr{T}, sptr)

include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
include("handle.jl")
include("info.jl")
include("comm.jl")
include("environment.jl")
include("datatypes.jl")
include("operators.jl")
include("pointtopoint.jl")
include("collective.jl")
include("topology.jl")
include("window.jl")

include("cman.jl")
include("deprecated.jl")

function __init__()
    @static if Sys.isunix()
        # need to open libmpi with RTLD_GLOBAL flag for Linux, before
        # any ccall cannot use RTLD_DEEPBIND; this leads to segfaults
        # at least on Ubuntu 15.10
        Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end

    if filesize(dlpath(libmpi)) != libmpi_size
        error("MPI library has changed, re-run Pkg.build(\"MPI\")")
    end
end

end


text:
struct Data
    mag::AbstractArray
    phase::AbstractArray
    header
    TEs::AbstractArray
end

mutable struct Options
    σ::AbstractArray
    unwrapping::Symbol
    mode::Symbol
    level::Number
    combination_type::Union{Symbol, Pair{Symbol,Tuple{Symbol,Symbol}}}
    sensitivity::Union{AbstractArray, Nothing}
    writedir::Union{AbstractString, Nothing}
    writesteps::Bool
end

function Options(;σ = [4, 4, 2], unwrapping = :laplacian, mode = :tanh, level = 5, combination_type = :SNR, sensitivity = nothing, writedir = nothing, writesteps = false)
    Options(σ, unwrapping, mode, level, combination_type, sensitivity, writedir, writesteps)
end

Base.string(t::Pair{Symbol,Tuple{Symbol,Symbol}}) = "$(t[1])_$(t[2][1])_$(t[2][2])"
Base.string(s::Nothing) = "nothing"
function saveconfiguration(options)
    open(joinpath(options.writedir, "settings_swi.txt"), "w") do io
        for fname in fieldnames(typeof(options))
            val = getfield(options, fname)
            if !(typeof(val) <: AbstractArray)
                println(io, "$fname: " * string(val))
            end
        end
    end
end


text:
# Control whether expensive assertions are enabled or not
if get(ENV, "DFTK_PROFILE", "develop") == "performance"
    assert_expensive_enabled() = false
else
    assert_expensive_enabled() = true
end

macro assert_expensive(expr)
    esc(:(if $(@__MODULE__).assert_expensive_enabled()
          @assert($expr) end))
end

macro assert_expensive(expr, text)
    esc(:(if $(@__MODULE__).assert_expensive_enabled()
          @assert($expr, $text) end))
end


text:
using BridgeSPDE
using Test

@testset "BridgeSPDE.jl" begin
    # Write your own tests here.
end


text:
"""
    Diagnostics

Accumulate mean fields and covariance statistics on the computational grid.
"""
module Diagnostics

export DiagnosticsGroup,
    setup_atmos_default_diagnostics,
    setup_atmos_core_diagnostics,
    setup_atmos_default_perturbations,
    setup_atmos_refstate_perturbations,
    setup_atmos_turbulence_stats,
    setup_atmos_mass_energy_loss,
    setup_atmos_spectra_diagnostics,
    setup_dump_state_diagnostics,
    setup_dump_aux_diagnostics

using CUDA
using Dates
using FileIO
using JLD2
using KernelAbstractions
using MPI
using OrderedCollections
using Printf
using StaticArrays
import KernelAbstractions: CPU

using ..ConfigTypes
using ..DGMethods
using ..BalanceLaws
using ..Mesh.Interpolation
using ..MPIStateArrays
using ..VariableTemplates
using ..Writers
import ..GenericCallbacks
using ..TicToc
using ..Spectra

using CLIMAParameters
using CLIMAParameters.Planet: planet_radius

Base.@kwdef mutable struct Diagnostic_Settings
    mpicomm::MPI.Comm = MPI.COMM_WORLD
    param_set::Union{Nothing, AbstractParameterSet} = nothing
    dg::Union{Nothing, DGModel} = nothing
    Q::Union{Nothing, MPIStateArray} = nothing
    starttime::Union{Nothing, String} = nothing
    output_dir::Union{Nothing, String} = nothing
end
const Settings = Diagnostic_Settings()

"""
    init(mpicomm, param_set, dg, Q, starttime, output_dir)

Initialize the diagnostics collection module -- save the parameters into
`Settings`.
"""
function init(
    mpicomm::MPI.Comm,
    param_set::AbstractParameterSet,
    dg::DGModel,
    Q::MPIStateArray,
    starttime::String,
    output_dir::String,
)
    Settings.mpicomm = mpicomm
    Settings.param_set = param_set
    Settings.dg = dg
    Settings.Q = Q
    Settings.starttime = starttime
    Settings.output_dir = output_dir
end

include("variables.jl")
include("helpers.jl")
include("atmos_common.jl")
include("thermo.jl")
include("groups.jl")

"""
    __init()__

Module initialization function. Currently, only fills in all currently
defined diagnostic variables.
"""
function __init__()
    setup_variables()
end

end # module Diagnostics

max_forks_repo_name: jm-c/CLIMA


text:
const Counter = let
    init() = 0
    op(n, _) = n + 1
    value(n) = n
    (; init, op, value)
end

struct FrequencyAnalysis end

to_nothings(v) = fill(nothing, axes(v))

function (f::FrequencyAnalysis)(input::ProcessedLayer)
    positional = vcat(input.positional, Any[map(to_nothings, first(input.positional))])
    labels = set(input.labels, length(positional) => "count")
    augmented_input = ProcessedLayer(input; positional, labels)
    return groupreduce(Counter, augmented_input)
end

"""
    frequency()

Compute a frequency table of the arguments.
"""
frequency() = transformation(FrequencyAnalysis())


text:
export MountainCarEnv, ContinuousMountainCarEnv

struct MountainCarEnvParams{T}
    min_pos::T
    max_pos::T
    max_speed::T
    goal_pos::T
    goal_velocity::T
    power::T
    gravity::T
    max_steps::Int
end

Base.show(io::IO, params::MountainCarEnvParams) = print(
    io,
    join(["$p=$(getfield(params, p))" for p in fieldnames(MountainCarEnvParams)], ","),
)

function MountainCarEnvParams(;
    T = Float64,
    min_pos = -1.2,
    max_pos = 0.6,
    max_speed = 0.07,
    goal_pos = 0.5,
    goal_velocity = 0.0,
    power = 0.001,
    gravity = 0.0025,
    max_steps = 200,
)
    MountainCarEnvParams{T}(
        min_pos,
        max_pos,
        max_speed,
        goal_pos,
        goal_velocity,
        power,
        gravity,
        max_steps,
    )
end

mutable struct MountainCarEnv{A,T,ACT,R<:AbstractRNG} <: AbstractEnv
    params::MountainCarEnvParams{T}
    action_space::A
    observation_space::Space{Vector{ClosedInterval{T}}}
    state::Vector{T}
    action::ACT
    done::Bool
    t::Int
    rng::R
end

"""
    MountainCarEnv(;kwargs...)

# Keyword arguments
- `T = Float64`
- `continuous = false`
- `rng = Random.GLOBAL_RNG`
- `min_pos = -1.2`
- `max_pos = 0.6`
- `max_speed = 0.07`
- `goal_pos = 0.5`
- `max_steps = 200`
- `goal_velocity = 0.0`
- `power = 0.001`
- `gravity = 0.0025`
"""
function MountainCarEnv(;
    T = Float64,
    continuous = false,
    rng = Random.GLOBAL_RNG,
    kwargs...,
)
    if continuous
        params = MountainCarEnvParams(; goal_pos = 0.45, power = 0.0015, T = T, kwargs...)
    else
        params = MountainCarEnvParams(; T = T, kwargs...)
    end
    action_space = continuous ? ClosedInterval{T}(-1.0, 1.0) : Base.OneTo(3)
    env = MountainCarEnv(
        params,
        action_space,
        Space([params.min_pos..params.max_pos, -params.max_speed..params.max_speed]),
        zeros(T, 2),
        rand(action_space),
        false,
        0,
        rng,
    )
    reset!(env)
    env
end

ContinuousMountainCarEnv(; kwargs...) = MountainCarEnv(; continuous = true, kwargs...)

Random.seed!(env::MountainCarEnv, seed) = Random.seed!(env.rng, seed)

RLBase.action_space(env::MountainCarEnv) = env.action_space
RLBase.state_space(env::MountainCarEnv) = env.observation_space
RLBase.reward(env::MountainCarEnv{A,T}) where {A,T} = env.done ? zero(T) : -one(T)
RLBase.is_terminated(env::MountainCarEnv) = env.done
RLBase.state(env::MountainCarEnv) = env.state

function RLBase.reset!(env::MountainCarEnv{A,T}) where {A,T}
    env.state[1] = 0.2 * rand(env.rng, T) - 0.6
    env.state[2] = 0.0
    env.done = false
    env.t = 0
    nothing
end

function (env::MountainCarEnv{<:ClosedInterval})(a::AbstractFloat)
    @assert a in env.action_space
    env.action = a
    _step!(env, a)
end

function (env::MountainCarEnv{<:Base.OneTo{Int}})(a::Int)
    @assert a in env.action_space
    env.action = a
    _step!(env, a - 2)
end

function _step!(env::MountainCarEnv, force)
    env.t += 1
    x, v = env.state
    v += force * env.params.power + cos(3 * x) * (-env.params.gravity)
    v = clamp(v, -env.params.max_speed, env.params.max_speed)
    x += v
    x = clamp(x, env.params.min_pos, env.params.max_pos)
    if x == env.params.min_pos && v < 0
        v = 0
    end
    env.done =
        x >= env.params.goal_pos && v >= env.params.goal_velocity ||
        env.t >= env.params.max_steps
    env.state[1] = x
    env.state[2] = v
    nothing
end


text:
module Dolfin

using PyCall

dolfin = PyCall.pyimport_conda("dolfin", "dolfin")
operator = PyCall.pyimport_conda("operator", "operator")

dolfin_methods =
["ALE", "AdaptiveLinearVariationalSolver", "AdaptiveNonlinearVariationalSolver",
 "And", "Argument", "Arguments", "Assembler", "AssemblerBase", "AutoSubDomain",
 "BDF1", "BackwardEuler", "BasisFunction", "BlockMatrix", "BlockVector",
 "BoolParameter", "BoostGraphOrdering", "BoundaryMesh", "BoundingBoxTree",
 "BoundingBoxTree3D", "BoxMesh", "BrokenElement", "ButcherMultiStageScheme",
 "CN2", "CRITICAL", "Cell", "CellFunction", "CellFunctionBool",
"CellFunctionDouble", "CellFunctionInt", "CellFunctionSizet", "CellNormal",
"CellSize", "CellType", "CellVolume", "Circumradius", "Coefficient",
"Coefficients", "CompiledSubDomain", "Constant", "CoordinateMatrix",
"CrankNicolson", "DBG", "DEBUG", "DOLFIN_EPS", "DOLFIN_EPS_LARGE",
"DOLFIN_LINELENGTH", "DOLFIN_PI", "DOLFIN_SQRT_EPS", "DOLFIN_TERM_WIDTH",
"DOLFIN_VERSION_GIT", "DOLFIN_VERSION_MAJOR", "DOLFIN_VERSION_MICRO",
"DOLFIN_VERSION_MINOR", "DOLFIN_VERSION_RELEASE", "DOLFIN_VERSION_STRING",
"DefaultFactory", "DirichletBC", "Dn", "DofMap", "Domain", "DomainBoundary",
"DoubleArray", "DoubleParameter", "Dx", "DynamicMeshEditor", "ERK", "ERK1",
"ERK4", "ERROR", "ESDIRK3", "ESDIRK4", "Edge", "EdgeFunction",
"EdgeFunctionBool", "EdgeFunctionDouble", "EdgeFunctionInt",
"EdgeFunctionSizet", "EigenFactory", "EigenKrylovSolver",
"EigenLUSolver", "EigenMatrix", "EigenVector", "EnrichedElement",
"Equation", "ErrorControl", "Event", "ExplicitEuler",
"ExplicitMidPoint", "Expression", "Extrapolation", "Face",
"FaceFunction", "FaceFunctionBool", "FaceFunctionDouble",
"FaceFunctionInt", "FaceFunctionSizet", "Facet", "FacetArea",
"FacetCell", "FacetElement", "FacetFunction", "FacetFunctionBool",
"FacetFunctionDouble", "FacetFunctionInt", "FacetFunctionSizet",
"FacetNormal", "File", "FiniteElement", "FiniteElementBase", "Form",
"ForwardEuler", "Function", "FunctionAXPY", "FunctionAssigner",
"FunctionSpace", "FunctionSpaceBase", "FunctionSpaceFromCpp", "GRL1",
"GRL2", "GenericAdaptiveVariationalSolver", "GenericBoundingBoxTree",
"GenericDofMap", "GenericFile", "GenericFunction", "GenericLUSolver",
"GenericLinearAlgebraFactory", "GenericLinearOperator",
"GenericLinearSolver", "GenericMatrix", "GenericPreconditioner",
"GenericSparsityPattern", "GenericTensor", "GenericVector",
"GlobalParameters", "Graph", "GraphBuilder", "H1", "H2", "HCurl",
"HDF5Attribute", "HDF5File", "HDiv", "HierarchicalDirichletBC",
"HierarchicalErrorControl", "HierarchicalForm",
"HierarchicalFunction", "HierarchicalFunctionSpace",
"HierarchicalLinearVariationalProblem", "HierarchicalMesh",
"HierarchicalMeshFunctionBool", "HierarchicalMeshFunctionDouble",
"HierarchicalMeshFunctionInt", "HierarchicalMeshFunctionSizet",
"HierarchicalMeshFunctionUInt",
"HierarchicalNonlinearVariationalProblem", "INFO", "Identity",
"ImplicitEuler", "Index", "IndexSet", "IntArray",
"IntParameter", "Integral", "InteriorElement", "IntervalMesh",
"KrylovSolver", "L2", "LUSolver", "Lagrange",
"LagrangeInterpolator", "Legendre", "LinearAlgebraObject",
"LinearOperator", "LinearSolver", "LinearVariationalProblem",
"LinearVariationalSolver", "LocalAssembler", "LocalMeshData",
"LocalSolver", "MPI", "MPICH_IGNORE_CXX_SEEK", "MPIInfo",
"MPI_Comm", "Matrix", "Max", "MaxCellEdgeLength",
"MaxFacetEdgeLength", "Measure", "Mesh", "MeshColoring",
"MeshConnectivity", "MeshCoordinates", "MeshData",
"MeshDisplacement", "MeshDomains", "MeshEditor",
"MeshEntity", "MeshFunction", "MeshFunctionBool",
"MeshFunctionDouble", "MeshFunctionInt",
"MeshFunctionSizet", "MeshGeometry", "MeshHierarchy",
"MeshPartitioning", "MeshPointIntersection",
"MeshQuality", "MeshRenumbering", "MeshTopology",
"MeshTransformation", "MeshValueCollection",
"MeshValueCollectionBool", "MeshValueCollectionDouble",
"MeshValueCollectionInt", "MeshValueCollectionSizet", "Min",
"MinCellEdgeLength", "MinFacetEdgeLength", "MixedElement", "MixedFunctionSpace",
"MultiMesh", "MultiMeshAssembler", "MultiMeshDirichletBC", "MultiMeshDofMap",
"MultiMeshForm", "MultiMeshFunction", "MultiMeshFunctionSpace",
"MultiMeshSubSpace", "MultiStageScheme", "NewtonSolver", "NonlinearProblem",
"NonlinearVariationalProblem", "NonlinearVariationalSolver", "Not",
"OpenMpAssembler", "OptimisationProblem", "Or", "OuterProductCell",
"OuterProductElement", "OuterProductVectorElement", "PETScBaseMatrix",
"PETScFactory", "PETScKrylovSolver", "PETScLUSolver", "PETScLinearOperator",
"PETScMatrix", "PETScObject", "PETScOptions", "PETScPreconditioner",
"PETScSNESSolver", "PETScUserPreconditioner", "PETScVector", "PROGRESS",
"ParameterValue", "Parameters", "PeriodicBoundaryComputation",
"PermutationSymbol", "Point", "PointIntegralSolver", "PointSource",
"ProductCell", "ProductDomain", "Progress", "RK4", "RKSolver", "RL1", "RL2",
"RectangleMesh", "RestrictedElement", "Restriction", "RushLarsenScheme",
"SCOTCH", "SLEPcEigenSolver", "STLFactory", "STLMatrix", "Scalar",
"SparsityPattern", "SparsityPatternBuilder", "SpatialCoordinate",
"SpecialFacetFunction", "StringParameter", "SubDomain", "SubMesh", "SubSpace",
"SubSystemsManager", "SubsetIterator", "SystemAssembler", "TRACE", "Table",
"TableEntry", "TensorConstant", "TensorElement", "TensorFunctionSpace",
"TensorLayout", "TensorProductElement", "TestFunction", "TestFunctions",
"Time", "TimeSeries", "Timer", "TimingClear_clear", "TimingClear_keep",
"TimingType_system", "TimingType_user", "TimingType_wall", "TraceElement",
"TrialFunction", "TrialFunctions", "UFLException", "UIntArray", "UnitCubeMesh",
"UnitHexMesh", "UnitIntervalMesh", "UnitQuadMesh", "UnitSquareMesh",
"UnitTetrahedronMesh", "UnitTriangleMesh", "VTKPlotter", "Variable", "Vector",
"VectorConstant", "VectorElement", "VectorFunctionSpace", "VectorSpaceBasis",
"Vertex", "VertexFunction", "VertexFunctionBool", "VertexFunctionDouble",
"VertexFunctionInt", "VertexFunctionSizet", "WARNING", "XDMFFile",
"acos", "action",
"adapt", "adapt_markers", "add_logfile", "adjoint", "as_backend_type",
"as_cell", "as_domain", "as_matrix", "as_tensor", "as_ufl", "as_vector",
"asin", "assemble", "assemble_system", "assign", "atan", "atan_2", "avg",
 "bessel_I", "bessel_J", "bessel_K", "bessel_Y", "between",
"cell_avg", "cells", "cofac", "common", "compile_expressions",
"compile_extension_module", "compilemodules", "compute_edge_map",
"compute_vertex_map", "conditional", "cos", "cosh", "cpp",
"create_ufc_function_spaces", "cross", "curl", "dC", "dI", "dO", "dP", "dS",
"dS_h", "dS_v", "dc", "debug", "deprecation", "derivative", "det", "dev",
"diag", "diag_vector", "diff", "div", "dof_to_vertex_map", "dolfin",
"dolfin_error", "dolfin_version", "dorfler_mark", "dot", "down_cast",
"ds", "ds_b", "ds_t", "ds_tb", "ds_v", "dump_timings_to_xml", "dx", "e",
"edges", "elem_div", "elem_mult", "elem_op", "elem_pow", "energy_norm",
"entities", "eq", "erf", "error", "errornorm", "exp", "exterior_derivative",
"f", "faces", "facet", "facet_avg", "facets", "fem", "fem_solve",
"functional", "functions", "ge", "generate_error_control",
"generate_error_control_forms", "get_global_parameters", "get_handler",
"get_log_level", "get_logger", "get_tensor_type", "git_commit_hash", "grad",
"gt", "has_cholmod", "has_hdf5", "has_krylov_solver_method",
"has_krylov_solver_preconditioner", "has_linear_algebra_backend",
"has_lu_solver_method", "has_mpi", "has_openmp", "has_parmetis", "has_petsc",
"has_petsc4py", "has_scotch", "has_slepc", "has_slepc4py", "has_tao",
"has_type", "has_umfpack", "has_zlib", "hexahedron", "homogenize", "i",
"importhandler", "in_nullspace", "increase_order", "indices", "info",
"info_blue", "info_green", "info_red", "info_stream", "info_underline",
"init", "init_parent_edge_indices", "inner", "integral_types", "interactive",
"interpolate", "intersect", "interval", "inv", "io", "ipow", "j", "jit", "jump",
"k", "krylov_solver_methods", "krylov_solver_preconditioners", "l", "la",
"la_index_dtype", "la_solve", "le", "lhs", "linear_algebra_backends",
"linear_solver_methods", "list_krylov_solver_methods",
"list_krylov_solver_preconditioners", "list_linear_algebra_backends",
"list_linear_solver_methods", "list_lu_solver_methods", "list_timings", "ln",
"log", "lt", "lt_coordinate", "lu_solver_methods", "mark", "max_value",
"memory_usage", "mesh", "min_value", "monitor_memory_usage", "mpi_comm_self",
"mpi_comm_world", "multistage", "nabla_div", "nabla_grad", "ne", "near",
"norm", "normalize", "not_working_in_parallel", "old_init", "outer", "p",
"parameters", "perp", "plot", "product", "project", "q", "quadrilateral",
"r", "rand", "rank", "refine", "register_element", "register_integral_type",
"relabel", "replace", "replace_integral_domains", "residual", "rhs", "rot",
"s", "seed", "sensitivity_rhs", "set_log_active", "set_log_level",
"set_output_stream", "shape", "show_elements", "sign", "sin", "sinh",
"sizeof_la_index", "skew", "solve", "split", "sqr", "sqrt",
"supported_elements", "supported_elements_for_plotting", "sym", "system",
"tan", "tanh", "tear", "tetrahedron", "tic", "time", "timing", "timings",
"toc", "tr", "transpose", "triangle", "unit_matrices", "unit_matrix",
"unit_vector", "unit_vectors", "variable", "vertex", "vertex_to_dof_map",
"vertices", "warning", "zero"]

function __init__()
  for meth in dolfin_methods
    ex = string(meth, "=", "dolfin[\"", meth, "\"]")
    prog = parse(ex)
    eval(prog)
    eval(Expr(:export, parse(meth)))
  end
end

import Base:+, -,*,/, ==, <<
import PyCall:isequal
+(f::PyObject, g::PyObject) = operator[:add](f,g)
+(f::PyObject, g) = operator[:add](f,Constant(g))

-(f::PyObject, g::PyObject) = operator[:sub](f,g)
-(f::PyObject, g) = operator[:sub](f,Constant(g))

*(f::PyObject, g::PyObject) = operator[:mul](f,g)
*(f::PyObject, g) = operator[:mul](f,Constant(g))

/(f::PyObject, g::PyObject) = operator[:div](f,g)
/(f::PyObject, g) = operator[:div](f,Constant(g))

==(f::PyObject, g::PyObject) = operator[:eq](f,g)
-(f::PyObject) = operator[:neg](f)
<<(f::PyObject, g::PyObject) = operator[:lshift](f, g)

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 191b2740-c20f-11eb-1bc3-4554934741e5
using DataFrames, JLD2, FileIO, Plots, Images, Colors, Swalbe, CUDA, HTTP

# ╔═╡ b2ec63db-847f-4986-b980-4d45482c7020
md"# Illustrative simulation
"

# ╔═╡ 5dc45dab-d666-43d4-b993-4985c5fe0501
md"## Dewetting of the **HI-ERN** logo

In this notebook we try to dewet the **HI-ERN** logo.
To do so we initialize a flat film and put small perturbations on top.
Under the film is a substrate on which we have printed the logo.

"

# ╔═╡ f72a2bc4-1ed8-448c-99a5-7d8fe85fd35d
html"""

"""

# ╔═╡ 80930386-a39f-41d6-9b33-c743dbde3b95
md"After some tweaking to make it as binary as possible we can add that image as an array to our simulation.
There was a default red in inkscape, I just stick to it.
This allows us to change the contact angle inside the letters and make them more wettable, such lower contact angle."

# ╔═╡ 3b5552d9-692d-4851-b353-630e5c5449dd
html"""

"""

# ╔═╡ d77cdc31-10ce-4384-adf4-90c57d8919e7
# Work in progress scrap images from web
# img = load(HTTP.get("https://jugit.fz-juelich.de/s.zitz/timedependent_wettability/-/raw/master/Figures/logo_red.png?inline=false").body)

# ╔═╡ 7cc7da7e-c1da-4ba3-9d4d-af0d9d5bfa73
# The dewetting function
function dewet_hiern(sys::Swalbe.SysConst; ϵ=1e-3, h₀=1.0, device="GPU", T=Float64)
	println("Starting logo dewetting")
	# Memory allocation
    fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py = Swalbe.Sys(sys, device, false, T)
    if device == "CPU"
        for i in 1:sys.Lx, j in 1:sys.Ly
            height[i,j] = h₀ + ϵ * randn()
        end
    elseif device == "GPU"
        h = zeros(size(height))
        for i in 1:sys.Lx, j in 1:sys.Ly
            h[i,j] = h₀ + ϵ * randn()
        end
		logo = load("E:\\JuliaStuff\\Github_Gists\\807b9a7b2226e65643288df9a8cc1f46\\logo_red.png")
        # Lower contact angle inside the letters
		theta = 2/9 .- 1/18 .* red.(reverse(rot180(logo), dims=2))
		# Forward it to the GPU
		θ = CUDA.adapt(CuArray, theta)
        height = CUDA.adapt(CuArray, h)
    end
    Swalbe.equilibrium!(fout, height, velx, vely, vsq)
    ftemp .= fout
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(height)
            if verbos
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(pressure, height, dgrad, sys.γ, theta, sys.n, sys.m, sys.hmin, sys.hcrit)
        Swalbe.∇f!(h∇px, h∇py, pressure, dgrad, height)
        Swalbe.slippage!(slipx, slipy, height, velx, vely, sys.δ, sys.μ)
        # Forces are the pressure gradient and the slippage due to substrate liquid boundary conditions
        Fx .= h∇px .+ slipx
        Fy .= h∇py .+ slipy
        # New equilibrium
        Swalbe.equilibrium!(feq, height, velx, vely, vsq)
        Swalbe.BGKandStream!(fout, feq, ftemp, -Fx, -Fy)
        # New moments
        Swalbe.moments!(height, velx, vely, fout)
        # Measurements, in this case only snapshots of simulational arrays
    end
    return height
    CUDA.reclaim()
end

# ╔═╡ a902a01d-04f5-4f2e-9414-decdc8f3bfd9
begin
	sys = Swalbe.SysConst(Lx=383, Ly=885, n=3, m=2, hmin=0.7)
	h = dewet_hiern(sys)
end

# ╔═╡ 6c154892-8375-4160-badf-fdfb16ce947e

# ╔═╡ Cell order:
# ╠═b2ec63db-847f-4986-b980-4d45482c7020
# ╠═191b2740-c20f-11eb-1bc3-4554934741e5
# ╟─5dc45dab-d666-43d4-b993-4985c5fe0501
# ╟─f72a2bc4-1ed8-448c-99a5-7d8fe85fd35d
# ╟─80930386-a39f-41d6-9b33-c743dbde3b95
# ╟─3b5552d9-692d-4851-b353-630e5c5449dd
# ╠═d77cdc31-10ce-4384-adf4-90c57d8919e7
# ╠═7cc7da7e-c1da-4ba3-9d4d-af0d9d5bfa73
# ╠═a902a01d-04f5-4f2e-9414-decdc8f3bfd9
# ╠═6c154892-8375-4160-badf-fdfb16ce947e

max_stars_repo_name: Zitzeronion/NotebooksnStuff

max_issues_repo_name: Zitzeronion/NotebooksnStuff

max_forks_repo_name: Zitzeronion/NotebooksnStuff


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ROCmOpenCLRuntime_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ROCmOpenCLRuntime")
JLLWrappers.@generate_main_file("ROCmOpenCLRuntime", UUID("10ae2a08-2eea-53f8-8c20-eec175020e9f"))
end  # module ROCmOpenCLRuntime_jll


text:
include("basic.jl")
include("operators.jl")


text:
"""
Module này thực hiện tính toán với hàng đợi

# Cách dùng

- `m = TênMôHình(tham_số...)`: Tạo một mô hình, danh sách mô hình xem ở dưới
- `λe(m)`: Tính cường độ đến hiệu quả ``\\lambda_e``.
- `pn(m)`: Tính xác suất ``P\\left[N(t) = n\\right]``.
- `L(m)`: Lượng khách trung bình trong hệ thống
- `Lq(m)`: Lượng khách trung bình trong hàng đợi
- `W(m)`: Thời gian trung bình trong hệ thống
- `Wq(m)`: Thời gian trung bình trong hàng đợi

## Ví dụ
~~~julia
julia> include("$(@__FILE__)")
Main.MMCKModels

julia> using .MMCKModels

julia> m = MM1(6, 8)
MM1(6, 8, 0.75)

julia> L(m)
3.0

julia> Lq(m)
2.25

julia> W(m)
0.5
~~~

Các loại mô hình hàng đợi:
- Hiện có: `MM1`, `MM1K`
- Dự kiến: `MMC`

Viết `?TênMôHình` để xem cách dùng.
"""
module MMCKModels

export MM1, MM1K, MMC,
	L, Lq, W, Wq, λe, pn

abstract type AbstractMMCK end;

"""
	MM1(μ, λ)

Tạo mô hình M/M/1
"""
struct MM1 <: AbstractMMCK
	λ; μ; ρ
	function MM1(λ::Real, μ::Real)
		if λ/μ >= 1
			@error "ρ phải < 1"
		end
		new(λ, μ, λ/μ)
	end
	function MM1(λ, μ)
		new(λ, μ, λ/μ)
	end
end

pn(m::MM1, n::Int) = (1 - m.ρ) * m.ρ^n
L(m::MM1) = m.ρ / (1 - m.ρ)
Lq(m::MM1) = m.ρ ^ 2 / (1 - m.ρ)
W(m::MM1) = L(m) / m.λ
Wq(m::MM1) = Lq(m) / m.λ

"""
	MM1K(λ, μ, k)

Tạo mô hình `M/M/1/K`.
"""
struct MM1K <: AbstractMMCK
	λ; μ; k; ρ
	MM1K(λ, μ, k) = new(λ, μ, k, λ/μ)
end

function pn(m::MM1K, n::Int)
	ρ = m.ρ
	k = m.k
	if ρ == 1
		1 / (k + 1)
	else
		ρ^n * (1 - ρ) / (1 - ρ^(k+1))
	end
end

λe(m::MM1K) = m.λ * (1 - pn(m, m.k))
L(m::MM1K) = let ρ = m.ρ, k = m.k
	if ρ == 1
		k / 2
	else
		l = k + 1
		ρ * (1 + k*ρ^l - l*ρ^k) / (1 - ρ) / (1 - ρ^l)
	end
end
Lq(m::MM1K) = let ρ = m.ρ, k = m.k
	if ρ == 1
		k * (k - 1) / 2 / (k + 1)
	else
		L(m) - ρ * (1 - ρ^k) / (1 - ρ^(k+1))
	end
end
W(m::MM1K) = L(m) / λe(m)
Wq(m::MM1K) = W(m) - 1/μ

"""
	MMC(μ, λ, c)

Tạo mô hình M/M/c
"""
struct MMC <: AbstractMMCK
	μ; λ; c; r; ρ
	function MMC(μ, λ, c)
		new(μ, λ, c, λ/μ, λ/μ/c)
	end
end

function pn(m::MMC, n::Integer)
	μ = m.μ
	λ = m.λ
	c = m.c
	r = m.r
	ρ = m.ρ
	if n == 0
		1 / ((c * r^c) / factorial(c) / (c - r) + sum(r^n/factorial(i) for i = 0:c-1))
	elseif n < c
		pn(m, 0) * r^n / factorial(n)
	else
		pn(m, 0) * r^n / factorial(c) / c^(n-c)
	end
end

end


text:
"""
    $(SIGNATURES)

An in-memory DistributedFactorGraph based on LightGraphs.jl with parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
mutable struct LightDFG{T <: AbstractParams, V <: AbstractDFGVariable, F <:AbstractDFGFactor} <: AbstractDFG{T}
    g::FactorGraph{Int, V, F}
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
    blobStores::Dict{Symbol, AbstractBlobStore}
end

"""
    $(SIGNATURES)

Create an in-memory LightDFG with the following parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
function LightDFG{T,V,F}(g::FactorGraph{Int,V,F}=FactorGraph{Int,V,F}();
                           description::String="LightGraphs.jl implementation",
                           userId::String="DefaultUser",
                           robotId::String="DefaultRobot",
                           sessionId::String="Session_$(string(uuid4())[1:6])",
                           userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           solverParams::T=T(),
                           blobstores=Dict{Symbol, AbstractBlobStore}()) where {T <: AbstractParams, V <:AbstractDFGVariable, F<:AbstractDFGFactor}
   # Validate the userId, robotId, and sessionId
   !isValidLabel(userId) && error("'$userId' is not a valid User ID")
   !isValidLabel(robotId) && error("'$robotId' is not a valid Robot ID")
   !isValidLabel(sessionId) && error("'$sessionId' is not a valid Session ID")
   return LightDFG{T,V,F}(g, description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], solverParams, blobstores)
end

# LightDFG{T}(; kwargs...) where T <: AbstractParams = LightDFG{T,DFGVariable,DFGFactor}(;kwargs...)

"""
    $(SIGNATURES)

Create an in-memory LightDFG with the following parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
function LightDFG{T}(g::FactorGraph{Int,DFGVariable,DFGFactor}=FactorGraph{Int,DFGVariable,DFGFactor}();
                     kwargs...) where T <: AbstractParams
    return LightDFG{T,DFGVariable,DFGFactor}(g; kwargs...)
end

function LightDFG(g::FactorGraph{Int,DFGVariable,DFGFactor}=FactorGraph{Int,DFGVariable,DFGFactor}();
                  solverParams::T=NoSolverParams(), kwargs...) where T
    return LightDFG{T,DFGVariable,DFGFactor}(g; solverParams=solverParams, kwargs...)
end

LightDFG(description::String,
         userId::String,
         robotId::String,
         sessionId::String,
         userData::Dict{Symbol, String},
         robotData::Dict{Symbol, String},
         sessionData::Dict{Symbol, String},
         solverParams::AbstractParams,
         blobstores=Dict{Symbol, AbstractBlobStore}()) =
         LightDFG(FactorGraph{Int,DFGVariable,DFGFactor}(), description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], solverParams, blobstores)

LightDFG{T,V,F}(description::String,
                userId::String,
                robotId::String,
                sessionId::String,
                userData::Dict{Symbol, String},
                robotData::Dict{Symbol, String},
                sessionData::Dict{Symbol, String},
                solverParams::T,
                blobstores=Dict{Symbol, AbstractBlobStore}()) where {T <: AbstractParams, V <:AbstractDFGVariable, F<:AbstractDFGFactor} =
                LightDFG(FactorGraph{Int,V,F}(), description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], solverParams, blobstores)


text:
using PSUtil
using Documenter

makedocs(;
    modules=[PSUtil],
    sitename="PSUtil.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://johnnychen94.github.io/PSUtil.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/johnnychen94/PSUtil.jl",
)


text:
using RoadRunnerJulia
using Test

rr = RoadRunnerJulia.createRRInstance()

@testset "RoadRunnerJulia.jl" begin
    ant_str = """
        const Xo, X1
         Xo -> S1; k1*Xo - k2*S1
         S1 -> S2; k3*S2
         S2 -> X1; k4*S2

         Xo = 1;   X1 = 0
         S1 = 0;   S2 = 0
         k1 = 0.1; k2 = 0.56
         k3 = 1.2; k4 = 0.9
    """
    rr = loada(ant_str)
    ssValues = RoadRunnerJulia.computeSteadyStateValues(rr)
    sol = [0.178571429, -4.38678616e-27]
    diff = ssValues - sol
    @test abs(diff[1]) < 1e-7 && abs(diff[2]) < 1e-7
    # Write your own tests here.
end


text:
module WebStitcher

using WebIO, FileIO, Images
using Colors
using CSSUtil

# include("projective_transform.jl")
include("Photostitching.jl")

const redraw = js"""
function redraw(ctx, cross, x, y){

    var rect = ctx.canvas.getBoundingClientRect();

    ctx.clearRect(0, 0, rect.width, rect.height); // Clears the canvas
    ctx.drawImage(ctx.img, 0, 0, ctx.imw, ctx.imh);

    ctx.lineWidth = 1;
    ctx.strokeStyle = \"D9E9FF\";
    var pad = 5;
    if(cross){
        ctx.beginPath();

        ctx.moveTo(x, y + pad);
        ctx.lineTo(x, ctx.imh);

        ctx.moveTo(x, y - pad);
        ctx.lineTo(x, 0);

        ctx.moveTo(x + pad, y);
        ctx.lineTo(ctx.imw, y);

        ctx.moveTo(x - pad, y);
        ctx.lineTo(0, y);

        ctx.stroke();

    }
    ctx.lineWidth = 2;
    ctx.strokeStyle = \"#f00\";
    ctx.lineJoin = "round";

    for(var i = 0; i < ctx.xx.length; i++){
        ctx.beginPath();
        ctx.arc(ctx.xx[i], ctx.yy[i], 5, 5, 180)
        ctx.stroke()
    }
}
"""

function to_pix_coords(scaling, xy...)
    isempty(xy[1]) && return Matrix{Float64}(0, 0)
    N = minimum(length.(xy))
    Float64[xy[j][i]/scaling for i = 1:N, j = 1:2]
end

const mouse_x = WebIO.@js function (canvas, e)
    @var rect = canvas.getBoundingClientRect();
    @var scaleX = canvas.width / rect.width #relationship bitmap vs. element for X
    return (e.clientX - rect.left) * scaleX
end

const mouse_y = WebIO.@js function (canvas, e)
    @var rect = canvas.getBoundingClientRect();
    @var scaleY = canvas.height / rect.height
    return (e.clientY - rect.top) * scaleY
end

function pointpicker(image, id, active, isdone, w; width = 200, num_points = 5)

    imh, imw = size(image)
    scaling = (width / imw)
    newh = round(Int, scaling * imh)

    image_base64 = sprint() do io
        show(Base64EncodePipe(io), MIME"image/png"(), image)
    end;

    obs = Observable(w, "obs$id", string("data:image/png;base64,", image_base64))

    x_obs = Observable(w, "x_obs$id", [], sync = true)
    y_obs = Observable(w, "y_obs$id", [], sync = true)

    is_active = Observable(w, "is_active$id", active, sync = true)

    ondependencies(w, WebIO.@js function ()

        @var el = this.dom.querySelector($("#$id"))
        @var ctx = el.getContext("2d")

        @var image = @new Image();
        image.onload = function()
            ctx.drawImage(image, 0, 0, $width, $newh);
        end
        image.src = $obs[]

        ctx.img = image;
        ctx.imw = $width;
        ctx.imh = $newh;
        ctx.xx = []
        ctx.yy = []
        ctx.is_active = $is_active[]

        window.redraw = $redraw

    end)

    add_click = WebIO.@js function add_click(e, context)

        @var el = context.dom.querySelector($("#$id"))
        @var ctx = el.getContext("2d")
        @var canvas = ctx.canvas;

        @var x = $(mouse_x)(canvas, e)
        @var y = $(mouse_y)(canvas, e)

        if ctx.is_active && !($isdone[])
            ctx.xx.push(x)
            ctx.yy.push(y)
            $x_obs[] = ctx.xx.concat([])
            $y_obs[] = ctx.yy.concat([])
            window.redraw(ctx, false, 0.0, 0.0)
            $is_active[] = false
            if !$active && ctx.xx.length >= $num_points # last image
                $isdone[] = true
            end
        end
    end

    on_mouse = WebIO.@js function on_mouse(e, context)
        @var el = context.dom.querySelector($("#$id"))
        @var ctx = el.getContext("2d")
        @var canvas = ctx.canvas;

        @var x = $(mouse_x)(canvas, e)
        @var y = $(mouse_y)(canvas, e)
        if ctx.is_active && !($isdone[])
            window.redraw(ctx, true, x, y)
        end
    end

    onjs(is_active, WebIO.@js function (val)
        @var div = this.dom.querySelector($("#selection"*id))
        @var el = this.dom.querySelector($("#$id"))
        @var ctx = el.getContext("2d")
        ctx.is_active = val

        if val
            div.style.backgroundColor = "#b0c780"
            if $active # the first image
                div.textContent = "Select a point"
            else
                div.textContent = "Select corresponing point to picture 1"
            end
        else
            div.style.backgroundColor = "#fff"
            div.textContent = ""
            window.redraw(ctx, false, 0.0, 0.0) # clear cross
        end

    end)
    onjs(isdone, WebIO.@js function (val)
        @var div = this.dom.querySelector($("#selection"*id))
        if val
            div.textContent = "DONE!"
            div.style.backgroundColor = "#999"
        end
    end)
    app = dom"canvas"(
        id = id,
        attributes = Dict(
            :width => "$(width)px", :height => "$(newh)px",
        ),
        events = Dict(
            :mouseup => add_click,
            :mousemove => on_mouse
        )
    )

    widget = dom"div"(
        id = "outline"*id,
        app, style = Dict(
            :outline => "2px solid #000",
            :width => "$(width)px", :height => "$(newh)px",
        )
    )

    selected = dom"div"(
        active ? "Select a Point" : "",

        id = "selection"*id,
        style = Dict(
            :fontSize => "200%",
            :textAlign => "center",
            :color => "#fff",
            :padding =>"15px",
            :backgroundColor => active ? "#b0c780" : "#fff",
            :outline => "2px solid #000",
            :width => "$(width)px", :height => "50px"
        )
    )
    points = map((x, y)-> to_pix_coords(scaling, x, y), x_obs, y_obs)
    vbox(selected, widget), points, is_active
end

function getcorrespondences(img1, img2, n; width = 700)
    w = Widget()
    isdone = Observable(w, "isdone", false);
    w1, points1, is_active1 = pointpicker(img1, "test1", true, isdone, w, width = width, num_points = n);
    w2, points2, is_active2 = pointpicker(img2, "test2", false, isdone, w, width = width, num_points = n);
    WebIO.onjs(is_active1, WebIO.@js function (val)
        if !val
            $is_active2[] = true
        end
    end)
    WebIO.onjs(is_active2, WebIO.@js function (val)
        if !val
            $is_active1[] = true
        end
    end)
    (points1, points2), w(vbox(w1, w2))
end

export getcorrespondences, ImageStitcher, stitchImages

end # module


text:
module PUnits
import Base: promote_rule, convert, show, sqrt, +, *, -, /, ^, .*, ./, .^, ==
# UnitSymbol was added to support prefixes
type UnitSymbol
    sym::String # unit symbol, ie "m"
    pre::Int16      # integer representing power of 10 of prefix, ie 3 represents kilo
end
UnitSymbol(sym::AbstractString, pre::Int) = UnitSymbol(convert(String, sym), Int16(pre))
Base.hash(x::UnitSymbol) = hash("$(x.sym),$(x.pre)") # make UnitSymbol act nice with Dict
Base.isequal(x::UnitSymbol, y::UnitSymbol) = x.sym==y.sym && x.pre==y.pre
PrefixSystem = Dict{Int16, String}()
reset_prefix_system() = [pop!(PrefixSystem, k) for (k,v) in PrefixSystem]
type Prefix
    pre::Int16
end
function Prefix(sym::AbstractString, pre::Int)
    PrefixSystem[Int16(pre)] = convert(String, sym)
    return Prefix(Int16(pre))
end
# Unit keeps track of unit symbols and powers, it has no idea how a symbol relates to
# any other symbol
type Unit
        d::Dict{UnitSymbol, Float64} # dict of exponents for unitsymbols
end
function Unit(x::AbstractString, prefix::Int=0)
    z = Unit()
    z.d[UnitSymbol(x,prefix)] = 1.0
    return z
end
Unit() = Unit(Dict{UnitSymbol, Float64}())
const Unitless = Unit()

function remove_prefix(x::Unit)
    z = Unit()
    total_power = 0
    for (unitsymbol, exponent) in x.d
        z.d[UnitSymbol(unitsymbol.sym,0)] = exponent + defaultget(z.d, UnitSymbol(unitsymbol.sym,0))
        total_power += unitsymbol.pre*exponent
    end
    filter!((k,v)->v!=0,z.d)
    return 10^total_power, z # prefactor, unit without prefix
end
defaultget(x::Dict, key::UnitSymbol) = haskey(x, key) ? x[key] : Float64(0)
function *(x::Unit, y::Unit)
    z = Unit()
    combined_units = union(keys(x.d), keys(y.d))
    for u in combined_units
        # round to Float32 precision so we dont end up with m^1e-12
        # going all the way to Float32 precision is probably overkill
        new_exponent = Float64(Float32(defaultget(x.d,u)+defaultget(y.d,u))) # round to Float32 precision
        if new_exponent != 0.0
                z.d[u] = new_exponent
        end
    end
    return z
end
function ^(x::Unit,y::Integer) # there must be a shorter way to define ^
    z = Unit()
    for (k,v) in x.d
            z.d[k] = v*y
    end
    return z
end
function ^(x::Unit,y::Rational)
    z = Unit()
    for (k,v) in x.dNao
            z.d[k] = v*y
    end
    return z
end
function ^(x::Unit,y::Number)
    z = Unit()
    for (k,v) in x.d
            z.d[k] = v*y
    end
    return z
end

-(x::Unit) = x^-1.0
/(x::Unit, y::Unit) = *(x,y^-1.0)
==(x::Unit, y::Unit) = x.d==y.d
+(x::Unit, y::Unit) = error("cannot add Units")
-(x::Unit, y::Unit) = error("cannot subtract Units")

function *(x::Prefix, y::Unit)
    length(y.d) != 1 ? error("Prefix*Unit only allowed for single units, length(Unit.d) = $(y.d)") : 0
    for (unitsymbol, exponent) in y.d
        exponent != 1 ? error("Prefix*Unit only allowed for single units with exponent==1, y=$(y), exponent=$exponent") : 0
        z = Unit()
        z.d[UnitSymbol(unitsymbol.sym, unitsymbol.pre+x.pre)] = Float64(1)
        return z
    end
end

# printing Text

superscript(c::Char) = "⁰¹²³⁴⁵⁶⁷⁸⁹·⁻"[[1,4,6,8,10,13,16,19,22,25,28,30][search("0123456789.-", c)]]
superscript(x::Number) = map(superscript,repr(x))
prettyround(x::Float64) = x%1 == 0 ? Int64(round(x)) : x
function show(io::IO,x::Unit)
    if isempty(x.d)
        print(io, "Unitless")
    end
    for (unitsymbol, exponent) in x.d
        if exponent == 1
            show(io, unitsymbol)
            print(io, " ")
        elseif exponent != 0
            show(io, unitsymbol)
            print(io, superscript(prettyround(exponent)))
        end
    end
end
function show(io::IO, unitsymbol::UnitSymbol)
    if haskey(PrefixSystem, unitsymbol.pre)
        print(io, PrefixSystem[unitsymbol.pre]*unitsymbol.sym)
    else
        print(io,"(10"*superscript(unitsymbol.pre)*unitsymbol.sym*")")
    end
end
function show(io::IO, prefix::Prefix)
    print(io, "Prefix $(prefix.pre) => "*get(PrefixSystem, prefix.pre, "10"*superscript(prefix.pre)))
end

export Unit, Unitless, Prefix, remove_prefix

end


text:
function get_class_name(s::AbstractString)
  replace(title(s), " ", "")
end


text:
# Transactions
# A transaction is the primary type through which the database is accessed. A
# transaction can be a single request, or it can be held open through many
# requests as a means of batching jobs together.

using Compat

export transaction, rollback, commit

struct Transaction
  conn::Connection
  commit::AbstractString
  location::AbstractString
  statements::Vector{Statement}
end

# TODO: Provide a version that accepts statements.
function transaction(conn::Connection)
  url = connurl(conn, "transaction")
  headers = connheaders(conn)
  body = Dict("statements" => [ ])

  resp = Requests.post(url; headers=headers, json=body)
  if resp.status != 201
    error("Failed to connect to database ($(resp.status)): $(conn)\n$(resp)")
  end
  respdata = Requests.json(resp)
  respheaders = Requests.headers(resp)

  Transaction(conn, respdata["commit"], respheaders["Location"], Statement[])
end

@compat function (txn::Transaction)(cypher::AbstractString, params::Pair...;
    submit::Bool=false)
  append!(txn.statements, [Statement(cypher, Dict(params))])
  if submit
    url = txn.location
    headers = connheaders(txn.conn)
    body = Dict("statements" => txn.statements)

    resp = Requests.post(url; headers=headers, json=body)

    if resp.status != 200
      error("Failed to submit transaction ($(resp.status)): $(txn)\n$(resp)")
    end
    respdata = Requests.json(resp)

    empty!(txn.statements)
    result = Result(respdata["results"], respdata["errors"])

    result
  end
end

function commit(txn::Transaction)
  url = txn.commit
  headers = connheaders(txn.conn)
  body = Dict("statements" => txn.statements)

  resp = Requests.post(url; headers=headers, json=body)

  if resp.status != 200
    error("Failed to commit transaction ($(resp.status)): $(txn)\n$(resp)")
  end
  respdata = Requests.json(resp)

  Result(respdata["results"], respdata["errors"])
end

function rollback(txn::Transaction)
  url = txn.location
  headers = connheaders(txn.conn)

  resp = Requests.delete(url; headers=headers)
  if resp.status != 200
    error("Failed to rollback transaction ($(resp.status)): $(txn)\n$(resp)")
  end
end


text:
using TexasHoldem
import TexasHoldem
TH = TexasHoldem
struct BotRandom <: AbstractAI
    prob::Float64
end
function TH.player_option!(game::Game, player::Player{BotRandom}, ::CheckRaiseFold)
    prob = TH.life_form(player).prob
    if rand() < prob; check!(game, player)
    else
        amt = Int(round(rand()*bank_roll(player), digits=0))
        amt = TH.bound_raise(game.table, player, amt) # to properly bound raise amount
        raise_to!(game, player, amt)
    end
end
function TH.player_option!(game::Game, player::Player{BotRandom}, ::CallRaiseFold)
    prob = TH.life_form(player).prob
    if rand() < 1-prob
        if rand() < prob; call!(game, player)
        else # re-raise
            amt = Int(round(rand()*bank_roll(player), digits=0))
            amt = TH.bound_raise(game.table, player, amt) # to properly bound raise amount
            raise_to!(game, player, amt)
        end
    else
        fold!(game, player)
    end
end
function TH.player_option!(game::Game, player::Player{BotRandom}, ::CallAllInFold)
    prob = TH.life_form(player).prob
    if rand() < 1-prob
        if rand() < prob; call!(game, player)
        else; raise_all_in!(game, player) # re-raise
        end
    else
        fold!(game, player)
    end
end
function TH.player_option!(game::Game, player::Player{BotRandom}, ::CallFold)
    prob = TH.life_form(player).prob
    if rand() < 1-prob; call!(game, player)
    else; fold!(game, player)
    end
end


text:
using SuiteSparseMatrixCollection
using MatrixMarket
using SuiteSparseGraphBLAS
using BenchmarkTools
using SparseArrays
include("tc.jl")
include("pr.jl")
graphs = [
    "karate",
    "com-Youtube",
    "as-Skitter",
    "com-LiveJournal",
    "com-Orkut",
    "com-Friendster",
]

ssmc = ssmc_db()
matrices = filter(row -> row.name ∈ graphs, ssmc)
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true
for name ∈ graphs
    path = fetch_ssmc(matrices[matrices.name .== name, :])[1]
    G = GBMatrix(convert(SparseMatrixCSC{Float64}, MatrixMarket.mmread(joinpath(path, "$name.mtx"))))
    SuiteSparseGraphBLAS.gbset(G, SuiteSparseGraphBLAS.FORMAT, SuiteSparseGraphBLAS.BYROW)
    GC.gc()
    println("$name | $(size(G)) | $(nnz(G)) edges")
    for centrality in [PR, TC1, TC3]
        println("Benchmarking $(string(centrality)) on $(name)")
        result = @benchmark $centrality($G) samples=3 seconds=600
        show(stdout,MIME"text/plain"(),result)
    end
end


text:
# stat(2)

function fstat(fd)
    s = Ref(C.stat64())
#    @cerr C.fstat64(fd, s)
     C.__fxstat64(0, fd, s)
    return s[]
end

stat_mode(stat) = stat.st_mode & C.S_IFMT

# End of file: stat.jl


text:
using BusinessDays
using Test
import Dates

# Issue #18
Dict(Any[])

# Issue #30
# list holidays for all available calendars
for c in [:BRSettlement, :BrazilExchange, :USNYSE, :USGovernmentBond, :USSettlement, :CanadaTSX, :CanadaSettlement, :EuroZone, :UKSettlement, :AustraliaASX]
    listholidays(c, Dates.Date(1900,1,1), Dates.Date(2100,1,1))
end

# Types
bhc = BusinessDays.Brazil()
ushc = BusinessDays.USSettlement()
ukhc = BusinessDays.UKSettlement()
usnysehc = BusinessDays.USNYSE()
usgovbondhc = BusinessDays.USGovernmentBond()
targethc = BusinessDays.TARGET()
hc_composite_BR_USA = CompositeHolidayCalendar([BusinessDays.Brazil(), BusinessDays.USSettlement()])
all_calendars_vec = [bhc, ushc, ukhc, hc_composite_BR_USA, usnysehc, usgovbondhc]

# two different instances of the same HolidayCalendar subtype should be equal
@test bhc == BusinessDays.Brazil()
@test ushc == BusinessDays.USSettlement()
@test bhc != ushc
@test string(bhc) == "BusinessDays.BRSettlement"

# Check typing system
@test isa(bhc,HolidayCalendar)
@test isa(ushc,HolidayCalendar)
@test isa(ukhc,HolidayCalendar)

include("easter_dates.jl")
include("easter-min-max.jl")

###############
# findweekday
###############
# weekday values:
# const Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7
# see query.jl on Dates module
# See also dayofweek(dt) function.
# function findweekday(weekday_target :: Int, yy :: Int, mm:: Int, occurrence :: Int, ascending :: Bool )
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 1, true) == Dates.Date(2015, 07, 06)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 2, true) == Dates.Date(2015, 07, 13)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 3, true) == Dates.Date(2015, 07, 20)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 4, true) == Dates.Date(2015, 07, 27)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 5, true) == Dates.Date(2015, 08, 03)

@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 1, false) == Dates.Date(2015, 07, 27)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 2, false) == Dates.Date(2015, 07, 20)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 3, false) == Dates.Date(2015, 07, 13)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 4, false) == Dates.Date(2015, 07, 06)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 5, false) == Dates.Date(2015, 06, 29)

@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 1, true) == Dates.Date(2015, 07, 03)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 2, true) == Dates.Date(2015, 07, 10)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 3, true) == Dates.Date(2015, 07, 17)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 4, true) == Dates.Date(2015, 07, 24)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 5, true) == Dates.Date(2015, 07, 31)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 6, true) == Dates.Date(2015, 08, 07)

@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 1, false) == Dates.Date(2015, 07, 31)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 2, false) == Dates.Date(2015, 07, 24)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 3, false) == Dates.Date(2015, 07, 17)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 4, false) == Dates.Date(2015, 07, 10)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 5, false) == Dates.Date(2015, 07, 03)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 6, false) == Dates.Date(2015, 06, 26)

@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 1, true) == Dates.Date(2015, 07, 01)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 2, true) == Dates.Date(2015, 07, 08)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 3, true) == Dates.Date(2015, 07, 15)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 4, true) == Dates.Date(2015, 07, 22)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 5, true) == Dates.Date(2015, 07, 29)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 6, true) == Dates.Date(2015, 08, 05)

@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 1, false) == Dates.Date(2015, 07, 29)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 2, false) == Dates.Date(2015, 07, 22)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 3, false) == Dates.Date(2015, 07, 15)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 4, false) == Dates.Date(2015, 07, 08)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 5, false) == Dates.Date(2015, 07, 01)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 6, false) == Dates.Date(2015, 06, 24)

@test_throws AssertionError BusinessDays.findweekday(Dates.Wednesday, 2015, 07, -1, false)
@test_throws AssertionError BusinessDays.findweekday(Dates.Wednesday, 2015, 07, -1, true)

# Test this only on 64bit or higher systems
len = typemax(UInt32) + 1
if len > typemax(UInt32)
    d0 = Dates.Date(1950, 2, 1)
    d1 = d0 + Dates.Day(len)
    @test_throws AssertionError BusinessDays._create_bdays_cache_arrays(bhc, d0, d1)
end

# Create HolidayCalendar instances
hc_australiaasx = BusinessDays.AustraliaASX()
hc_australiaact = BusinessDays.Australia(:ACT)
hc_australiansw = BusinessDays.Australia(:NSW)
hc_australiant  = BusinessDays.Australia(:NT)
hc_australiaqld = BusinessDays.Australia(:QLD)
hc_australiasa  = BusinessDays.Australia(:SA)
hc_australiatas = BusinessDays.Australia(:TAS)
hc_australiawa  = BusinessDays.Australia(:WA)
hc_australiavic = BusinessDays.Australia(:VIC)
hc_brazil       = BusinessDays.Brazil()
hc_brazil_bmf   = BusinessDays.BrazilBMF()
hc_usa          = BusinessDays.USSettlement()
hc_uk           = BusinessDays.UKSettlement()
hc_usnyse       = BusinessDays.USNYSE()
hc_usgovbond    = BusinessDays.USGovernmentBond()
hc_canadatsx    = BusinessDays.CanadaTSX()
hc_canada       = BusinessDays.CanadaSettlement()

####################
# Calendar Tests
####################

usecache = false
include("calendar_tests.jl")

usecache = true
include("calendar_tests.jl")

####################
# Performance Tests
####################
include("perftests.jl")

BusinessDays.cleancache(hc_brazil)
BusinessDays.cleancache()

struct TestCalendar <: HolidayCalendar end
cc = TestCalendar()
@test_throws ErrorException isholiday(cc, Dates.Date(2015,1,1))
BusinessDays.isholiday(::TestCalendar, dt::Dates.Date) = dt == Dates.Date(2015,8,27)
@test isholiday(cc, Dates.Date(2015,8,26)) == false
@test isholiday(cc, Dates.Date(2015,8,27)) == true
BusinessDays.initcache(cc)
@test isholiday(cc, Dates.Date(2015,8,26)) == false
@test isholiday(cc, Dates.Date(2015,8,27)) == true
BusinessDays.cleancache(cc)
isholiday(:Brazil, Dates.Date(2016,2,1))
isholiday(:TestCalendar, Dates.Date(2016,2,1))
isholiday("TestCalendar", Dates.Date(2016,2,1))
isbday(:Brazil, Dates.Date(2016,2,1))
isbday(:TestCalendar, Dates.Date(2016,2,1))
isbday("TestCalendar", Dates.Date(2016,2,1))

sym_vec = [:Brazil, :UKSettlement]
BusinessDays.initcache(sym_vec)
BusinessDays.cleancache(sym_vec)

str_vec = ["Brazil", "UKSettlement", "Canada", "UnitedStates", "USNYSE", "USGovernmentBond", "CanadaTSX", "WeekendsOnly"]
BusinessDays.initcache(str_vec)
BusinessDays.cleancache(str_vec)

BusinessDays.initcache("UKSettlement", Dates.Date(2000,1,1), Dates.Date(2000,5,2))
BusinessDays.initcache("UKSettlement", Dates.Date(2000,1,1), Dates.Date(2000,5,2)) # repeating initcache should work
BusinessDays.initcache("UKSettlement", Dates.Date(2000,1,1), Dates.Date(2000,1,1)) # single date cache should work
BusinessDays.cleancache("UKSettlement")

# equality and hash function for Australia
@test BusinessDays.Australia(:ACT) == BusinessDays.Australia(:ACT)
@test BusinessDays.Australia(:ACT) != BusinessDays.Australia(:NSW)
@test hash(BusinessDays.Australia(:ACT)) == hash(BusinessDays.Australia(:ACT))
@test hash(BusinessDays.Australia(:ACT)) != hash(BusinessDays.Australia(:NSW))

BusinessDays.initcache(BusinessDays.Australia(:ACT))
@test haskey(BusinessDays.CACHE_DICT, BusinessDays.Australia(:ACT))
@test !haskey(BusinessDays.CACHE_DICT, BusinessDays.Australia(:NSW))
BusinessDays.cleancache(BusinessDays.Australia(:ACT))

include("customcalendar-example.jl")

#########################
# GenericHolidayCalendar
#########################

d0 = Dates.Date(1980,1,1)
d1 = Dates.Date(2050,1,1)
gen_brazil = GenericHolidayCalendar(listholidays(hc_brazil, d0, d1))
@test isbday(gen_brazil, Dates.Date(2014, 12, 31)) == true # wednesday
@test isbday(gen_brazil, Dates.Date(2015, 01, 01)) == false # new year
@test isbday(gen_brazil, Dates.Date(2015, 01, 02)) == true # friday

@test advancebdays(gen_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5]) == advancebdays(hc_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5])
@test advancebdays(gen_brazil, Dates.Date(2015,9,1), 0:5) == advancebdays(hc_brazil, Dates.Date(2015,9,1), 0:5)
@test listholidays(gen_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30)) == listholidays(hc_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30))

@test tobday(gen_brazil, Dates.Date(2013, 02, 08)) == tobday(hc_brazil, Dates.Date(2013, 02, 08))
@test tobday(gen_brazil, Dates.Date(2013, 02, 09)) == tobday(hc_brazil, Dates.Date(2013, 02, 09))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(1900,2,1), Dates.Date(2000,2,1))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(2000,2,1), Dates.Date(2100,2,1))

d0_test = Dates.Date(1980,1,2)
d1_test = Dates.Date(2049,12,28)
@test bdays(hc_brazil, d0_test, d1_test) == bdays(gen_brazil, d0_test, d1_test)
@test bdayscount(hc_brazil, d0_test, d1_test) == bdayscount(gen_brazil, d0_test, d1_test)
println("a million with GenericHolidayCalendar...")
@time for i in 1:1000000 bdays(gen_brazil, d0_test, d1_test) end
BusinessDays.cleancache(gen_brazil)

# Start all over, but without cache
gen_brazil = GenericHolidayCalendar(listholidays(hc_brazil, d0, d1), d0, d1, false)

@test isbday(gen_brazil, Dates.Date(2014, 12, 31)) == true # wednesday
@test isbday(gen_brazil, Dates.Date(2015, 01, 01)) == false # new year
@test isbday(gen_brazil, Dates.Date(2015, 01, 02)) == true # friday

@test advancebdays(gen_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5]) == advancebdays(hc_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5])
@test advancebdays(gen_brazil, Dates.Date(2015,9,1), 0:5) == advancebdays(hc_brazil, Dates.Date(2015,9,1), 0:5)
@test listholidays(gen_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30)) == listholidays(hc_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30))

@test tobday(gen_brazil, Dates.Date(2013, 02, 08)) == tobday(hc_brazil, Dates.Date(2013, 02, 08))
@test tobday(gen_brazil, Dates.Date(2013, 02, 09)) == tobday(hc_brazil, Dates.Date(2013, 02, 09))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(1900,2,1), Dates.Date(2000,2,1))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(2000,2,1), Dates.Date(2100,2,1))

@test BusinessDays.needs_cache_update(gen_brazil, d0, d1) == false
BusinessDays.initcache(gen_brazil)
@test BusinessDays.needs_cache_update(gen_brazil, d0, d1) == true
@test isbday(gen_brazil, Dates.Date(2014, 12, 31)) == true # wednesday
@test isbday(gen_brazil, Dates.Date(2015, 01, 01)) == false # new year
@test isbday(gen_brazil, Dates.Date(2015, 01, 02)) == true # friday

# does nothing, because cache is already there
BusinessDays.initcache(gen_brazil)

# A GenericHolidayCalendar is defined by its set of holidays, dtmin, dtmax
dtmin = Dates.Date(2018,1,15)
dtmax = Dates.Date(2018,1,19)
gen_1 = GenericHolidayCalendar(Set([Dates.Date(2018,1,16), Dates.Date(2018,1,18)]), dtmin, dtmax)
gen_2 = GenericHolidayCalendar(Set([Dates.Date(2018,1,16), Dates.Date(2018,1,18)]), dtmin, dtmax)
@test gen_1 == gen_2

# On a set, they are the same element
cal_set = Set([gen_1, gen_2])
@test length(cal_set) == 1

# On a dict, they represent the same key
cal_dict = Dict(gen_1 => "hey")
@test cal_dict[gen_2] == "hey"

# Broadcast
@test BusinessDays.tobday.(BusinessDays.WeekendsOnly(), [Dates.Date(2019, 8, 30), Dates.Date(2019, 8, 31), Dates.Date(2019, 9, 1), Dates.Date(2019, 9, 2)]) == [Dates.Date(2019, 8, 30), Dates.Date(2019, 9, 2), Dates.Date(2019, 9, 2), Dates.Date(2019, 9, 2)]


text:
using Documenter, POMDPs

makedocs(
    modules = [POMDPs],
    format = :html,
    sitename = "POMDPs.jl"
)

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


text:
p = pwd()
cd(Pkg.dir("EiSCor/deps/eiscor-master"))
run(`make`)
is_apple()   ? run(`cp libeiscor.dylib.0.1.0 ../libeiscor.dylib`) : nothing
is_linux()   ? run(`cp libeiscor.so.0.1.0 ../libeiscor.so`) : nothing
is_windows() ? run(`cp libeiscor.dll.0.1.0 ../libeiscor.dll`) : nothing
cd(p)


text:
# Some additional checks compared to
# test_lord_cressie.jl

using Empirikos
using Test
using Hypatia
using Random

Random.seed!(1)
Zs_standard = StandardNormalSample.(randn(1000) .* sqrt(2))
true_prior = Normal(0, 1)

g_prior = MixturePriorClass( Normal.(-2:0.1:2, 1))

gauss_loc_02 = InfinityNormDensityBand(; a_min = -3.0, a_max = 3.0, α = 0.2)

# deepcopying below to make sure RNG status remains the same.
fitted_gauss_loc_02 = fit(gauss_loc_02, Zs_standard)

floc_interval = FLocalizationInterval(;flocalization = gauss_loc_02,
                                       solver = Hypatia.Optimizer,
                                       convexclass = g_prior)

floc_interval_prefitted = FLocalizationInterval(;flocalization = fitted_gauss_loc_02,
                                       solver = Hypatia.Optimizer,
                                       convexclass = g_prior)

target1 = PosteriorMean(StandardNormalSample(1.0))

@test target1(true_prior) == 0.5

# check if result is same with unfitted or prefitted flocalization

confint1 = confint(floc_interval, target1, Zs_standard)
confint1_prefitted = confint(floc_interval_prefitted, target1, Zs_standard)
confint1_prefitted2 = confint(floc_interval_prefitted, target1)

@test confint1 == confint1_prefitted
@test confint1 == confint1_prefitted2

# For KDE-F-Localization, check whether solutions with LinearVexity and ConvexVexity match up to numerical error.

confint1_prefitted_convexvexity = confint(floc_interval_prefitted, target1,
                                        Zs_standard, Empirikos.ConvexVexity())

@test confint1.α == confint1_prefitted_convexvexity.α
@test confint1.lower ≈ confint1_prefitted_convexvexity.lower rtol=0.005
@test confint1.upper ≈ confint1_prefitted_convexvexity.upper rtol=0.005

# For multiple targets, check if broadcasting works correctly.

target2 = PosteriorMean(StandardNormalSample(2.0))
@test target2(true_prior) == 1.0

confint2 = confint(floc_interval, target2, Zs_standard)

both_targets = [target1; target2]

both_cis = confint.(floc_interval, both_targets, Zs_standard)
@test both_cis[1] == confint1
@test both_cis[2] == confint2

both_cis_prefitted = confint.(floc_interval_prefitted, both_targets)
@test both_cis_prefitted == both_cis

both_cis_level = confint.(floc_interval, both_targets, Zs_standard; level=0.8)
@test both_cis == both_cis_level

both_cis_prefitted_level = confint.(floc_interval_prefitted, both_targets; level=0.8)
@test both_cis_prefitted_level == both_cis

# Check if level is double-checked when using `confint` command.
@test_throws ErrorException confint(floc_interval, target1; level = 0.95)
@test_throws ErrorException confint.(floc_interval, both_targets, Zs_standard; level=0.95)
@test_throws ErrorException confint.(floc_interval_prefitted, both_targets; level=0.1)

# When target is marginal density, check whether it works for linear functionals.

density_target = MarginalDensity(StandardNormalSample(0.0))
estimated_density_at_1 = fitted_gauss_loc_02.interp_kde.itp(0.0)
band_width = fitted_gauss_loc_02.C∞
true_density_at_1 = density_target(true_prior)

infinity_band_ub = estimated_density_at_1 + band_width
infinity_band_lb = estimated_density_at_1 - band_width

density_ci = confint(floc_interval_prefitted, density_target)

@test infinity_band_ub >= density_ci.upper
@test infinity_band_lb <= density_ci.lower

# now use a more-fine grained prior class so that we can get exact equality
floc_interval_discrete = FLocalizationInterval(;flocalization = fitted_gauss_loc_02,
                                       solver = Hypatia.Optimizer,
                                       convexclass = DiscretePriorClass(-3:0.01:3))

density_ci_discrete = confint(floc_interval_discrete, density_target)

@test infinity_band_ub ≈ density_ci_discrete.upper rtol=0.001
@test infinity_band_lb ≈ density_ci_discrete.lower rtol=0.001

# Intervals are nested for larger level.

gauss_loc_01 = InfinityNormDensityBand(; a_min = -3.0, a_max = 3.0, α = 0.1)

floc_interval_01 = FLocalizationInterval(;flocalization = gauss_loc_01,
                                       solver = Hypatia.Optimizer,
                                       convexclass = g_prior)
confint1_01 = confint(floc_interval_01, target1, Zs_standard)

@test confint1_01.α == 0.1
@test confint1_01.lower < confint1.lower
@test confint1_01.upper > confint1.upper


text:
@testset "BoundingBox" begin
    @test BoundingBox(1.0,2,3,4,5,6) === BoundingBox(1.0,2.0,3.0,4.0,5.0,6.0)

    @test eltype(BoundingBox(1,2,3,4,5,6)) === SVector{3, Int}

    # making bounding boxes by encapsulating geometry
    @test boundingbox(SVector(1,2,3)) === BoundingBox(1,2,3,1,2,3)
    @test boundingbox(SVector(1,2,3), SVector(4,5,6)) === BoundingBox(1,2,3,4,5,6)
    @test boundingbox(SVector(4,5,6), SVector(1,2,3)) === BoundingBox(1,2,3,4,5,6)
    @test boundingbox(BoundingBox(0,0,0,1,1,1), BoundingBox(1,1,1,2,2,2)) === BoundingBox(0,0,0,2,2,2)
    @test boundingbox(BoundingBox(0,0,0,1,1,1), BoundingBox(1.0,1.0,1.0,2.0,2.0,2.0)) === BoundingBox(0.0,0.0,0.0,2.0,2.0,2.0)
    @test boundingbox(BoundingBox(0,0,0,1,1,1), SVector(2,2,2)) === BoundingBox(0,0,0,2,2,2)
    @test boundingbox([SVector(1,2,0), SVector(0,1,2), SVector(2,0,1)]) === BoundingBox(0,0,0,2,2,2)
    @test boundingbox([BoundingBox(0,0,0,1,1,1), BoundingBox(1.0,1.0,1.0,2.0,2.0,2.0)]) === BoundingBox(0.0,0.0,0.0,2.0,2.0,2.0)

    # 27 octants
    @test SVector(1,2,3) ∈ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,2,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,2,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(1,-1,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,-1,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,-1,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(1,5,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,5,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,5,3) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))

    @test SVector(1,2,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,2,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,2,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(1,-1,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,-1,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,-1,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(1,5,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,5,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,5,-1) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))

    @test SVector(1,2,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,2,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,2,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(1,-1,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,-1,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,-1,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(1,5,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(-1,5,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))
    @test SVector(5,5,5) ∉ boundingbox(SVector(0,0,0), SVector(4,4,4))

    # non-empty intersection
    @test intersects(BoundingBox(0,0,0,1,1,1), BoundingBox(-1,-1,-1,2,2,2))
    @test intersects(BoundingBox(0,0,0,1,1,1), BoundingBox(0.25,0.25,0.25,0.75,0.75,0.75))
    @test intersects(BoundingBox(0,0,0,1,1,1), BoundingBox(0.5,0.5,0.5,1.5,1.5,1.5))
    @test intersects(BoundingBox(0,0,0,1,1,1), BoundingBox(-0.5,-0.5,-0.5,0.5,0.5,0.5))

    # isempty
    @test !isempty(BoundingBox(0,0,0, 1, 1, 1))
    @test  isempty(BoundingBox(0,0,0,-1, 1, 1))
    @test  isempty(BoundingBox(0,0,0, 1,-1, 1))
    @test  isempty(BoundingBox(0,0,0, 1, 1,-1))
    @test  isempty(BoundingBox(0,0,0,-1,-1, 1))
    @test  isempty(BoundingBox(0,0,0,-1, 1,-1))
    @test  isempty(BoundingBox(0,0,0, 1,-1,-1))
    @test  isempty(BoundingBox(0,0,0,-1,-1,-1))

    # pad
    @test pad(BoundingBox(1,2,3,4,5,6), 1) === BoundingBox(0,1,2,5,6,7)

    # wireframe
    @test wireframe(BoundingBox(0,0,0,1,1,1)) == [Line(SVector(0,0,0), SVector(1,0,0)),
                                                  Line(SVector(0,1,0), SVector(1,1,0)),
                                                  Line(SVector(0,0,1), SVector(1,0,1)),
                                                  Line(SVector(0,1,1), SVector(1,1,1)),
                                                  Line(SVector(0,0,0), SVector(0,1,0)),
                                                  Line(SVector(1,0,0), SVector(1,1,0)),
                                                  Line(SVector(0,0,1), SVector(0,1,1)),
                                                  Line(SVector(1,0,1), SVector(1,1,1)),
                                                  Line(SVector(0,0,0), SVector(0,0,1)),
                                                  Line(SVector(1,0,0), SVector(1,0,1)),
                                                  Line(SVector(0,1,0), SVector(0,1,1)),
                                                  Line(SVector(1,1,0), SVector(1,1,1))]
end