text:
import Base: getindex, zero, show, -, +, ==, !=, *, /, view

struct WannierFunction{N, T<:AbstractFloat} <: AbstractArray{SVector{N, Complex{T}}, 3}
	points::Array{Point{3, T}, 3}
	values::Array{SVector{N, Complex{T}}, 3}
end

function WannierFunction(filename::AbstractString, points::Array{Point3{T}, 3}=read_points_from_xsf(filename)) where {T <: AbstractFloat}
	re = read_values_from_xsf(T, filename)
	values = [SVector(complex(a)) for a in re]
	return normalize(WannierFunction(points, values))
end

function WannierFunction(filename_re::String, filename_im::String, points::Array{Point3{T}, 3} = read_points_from_xsf(filename_re)) where {T <: AbstractFloat}
	re, im = read_values_from_xsf.(T, (filename_re, filename_im))
	values = [SVector(Complex(a, b)) for (a, b) in zip(re, im)]
	return normalize(WannierFunction(points, values))
end

function WannierFunction(filename_up_re::String, filename_up_im::String, filename_down_re::String, filename_down_im::String, points::Array{Point3{T}, 3} = read_points_from_xsf(filename_up_re)) where {T <: AbstractFloat}
	up_re, up_im, down_re, down_im =
		read_values_from_xsf.(T, (filename_up_re, filename_up_im, filename_down_re, filename_down_im))
	values = [SVector(Complex(a, b), Complex(c, d)) for (a, b, c, d) in zip(up_re, up_im, down_re, down_im)]
	return normalize(WannierFunction(points, values))
end

WannierFunction(point_func::Function,  points::Array) =
    normalize(WannierFunction(points, point_func.(points)))

for f in (:size, :getindex, :setindex!)
	@eval @inline @propagate_inbounds Base.$f(x::WannierFunction, i...) =
		Base.$f(x.values, i...)
end

for f in (:length, :stride, :ndims, :axes, :strides)
	@eval @inline Base.$f(w::WannierFunction) = Base.$f(w.values)
end

Base.similar(x::WannierFunction,::Type{S}) where S =
  WannierFunction(x.points, similar(x.values, S))

Base.unsafe_convert(T::Type{<:Ptr}, x::WannierFunction) =
	unsafe_convert(T, x.values)

Base.Broadcast.broadcastable(w::WannierFunction) =
	w.values

#### LinearAlgebra overloads
function LinearAlgebra.adjoint(w::WannierFunction)
	out = WannierFunction(w.points, similar(w.values))
	adjoint!(out, w)
end

LinearAlgebra.adjoint!(w1::WannierFunction, w2::WannierFunction) =
	w1 .= adjoint.(w2)

function LinearAlgebra.dot(w1::WannierFunction{T}, w2::WannierFunction{T}) where {T}
    s = zero(T)
    for (v1, v2) in zip(w1.values, w2.values)
        s += v1' * v2
    end
    return real(s)
end

function LinearAlgebra.dot(v::Vector, wfs::Vector{<:WannierFunction})
    res = WannierFunction(wfs[1].points, zeros(eltype(wfs[1].values), size(wfs[1].values)))
    for ic in 1:length(v)
        res .+= v[ic] .* wfs[ic]
    end
    return res
end
LinearAlgebra.dot(wfs::Vector{<:WannierFunction}, v::Vector) =
    dot(v, wfs)

LinearAlgebra.norm(wfc::WannierFunction) =
	dot(wfc, wfc)

LinearAlgebra.normalize!(wfc::WannierFunction) =
	wfc ./= sqrt(norm(wfc))

same_grid(w1::WannierFunction, w2::WannierFunction) =
	w1.points === w2.points

function wan_op(op::Function, w1::W, w2::W) where {W <: WannierFunction}
	@assert same_grid(w1, w2) "Wannier functions are not defined on the same grid"
	return WannierFunction(w1.points, op(w1.values, w2.values))
end

+(w1::WannierFunction, w2::WannierFunction) = wan_op(+, w1, w2)
*(w1::WannierFunction, w2::WannierFunction) = wan_op(*, w1, w2)
-(w1::WannierFunction, w2::WannierFunction) = wan_op(-, w1, w2)

*(w1::WannierFunction, n::Number) = WannierFunction(w1.points, w1.values .* n)
*(n::Number, w1::WannierFunction) = WannierFunction(w1.points, n .* w1.values)
/(w1::WannierFunction, n::Number) = WannierFunction(w1.points, n ./ w1.values)
/(n::Number, w1::WannierFunction) = WannierFunction(w1.points, w1.values ./ n)

LinearAlgebra.dot(w1::WannierFunction, n::Number) = w1 * n
LinearAlgebra.dot(n::Number, w1::WannierFunction) = n * w1

function generate_wannierfunctions(k_filenames::Vector{String}, chk_info, wannier_plot_supercell::NTuple{3,Int}, wan_plot_list=1:chk_info.n_wann)
    num_kpts = length(chk_info.kpoints)
    U = permutedims(chk_info.U_matrix, (2, 1, 3))
    U_opt = permutedims(chk_info.U_matrix_opt,(2,1,3))
    tu = read_unk(k_filenames[1])
    nrx, nry, nrz = size(tu,1), size(tu,2), size(tu,3)
    supercell_xrange = -div(wannier_plot_supercell[1],2)*nrx : div(wannier_plot_supercell[1] + 1, 2)*nrx - 1
    supercell_yrange = -div(wannier_plot_supercell[2],2)*nry : div(wannier_plot_supercell[2] + 1, 2)*nry - 1
    supercell_zrange = -div(wannier_plot_supercell[3],2)*nrz : div(wannier_plot_supercell[3] + 1, 2)*nrz - 1

    nx, ny, nz = length.((supercell_xrange, supercell_yrange, supercell_zrange))

    nwfun = length(wan_plot_list)
    wfuncs_all = zeros(eltype(tu), nwfun, (wannier_plot_supercell .* (nrx, nry, nrz))...,size(tu, 5))
    n_wann = chk_info.n_wann
    r_wan      = zeros(eltype(tu), chk_info.n_wann, nrx, nry, nrz)

    p = Progress(length(chk_info.kpoints))
    @inbounds for ik = 1:length(chk_info.kpoints)
        k = chk_info.kpoints[ik]
        unk_all = read_unk(k_filenames[ik])
        for is = 1:size(tu, 5)
            u = U[wan_plot_list, :, ik]
            u_opt = U_opt[:, :, ik]
            inc_ids = findall(!iszero, chk_info.lwindow[:, ik])
            fill!(r_wan, 0.0)
            for ib in 1:chk_info.ndimwin[ik]
                iib = inc_ids[ib]
                Threads.@threads for nz in 1:nrz
                    for ny in 1:nry
                        for nx in 1:nrx
                            @simd for iw in 1:n_wann
                                r_wan[iw, nx, ny, nz] += u_opt[iw, ib] * unk_all[nx,ny,nz,iib,is]
                            end
                        end
                    end
                end
            end
            Threads.@threads for iisz in 1:nz
                isz = supercell_zrange[iisz]
                iz = mod1(isz, nrz)
                for iisy in 1:ny
                    isy = supercell_yrange[iisy]
                    iy = mod1(isy, nry)
                    for iisx in 1:nx
                        isx = supercell_xrange[iisx]
                        ix = mod1(isx, nrx)
                        scalfac = exp(2im*π*dot(k, Vec3((isx-1)/nrx, (isy-1)/nry, (isz-1)/nrz)))
                        for ib in 1:n_wann
                            rt = r_wan[ib, ix, iy, iz] * scalfac
                            for iw in 1:nwfun
                                wfuncs_all[iw, iisx, iisy, iisz, is] += u[iw, ib] * rt
                            end
                        end
                    end
                end
            end
        end
        next!(p)
    end
    wfuncs_all ./= num_kpts
    if size(tu, 5) == 1
        Threads.@threads for iw = 1:size(wfuncs_all, 1)
            tmaxx = 0.0
            cmod  = 1.0+0.0im
            for iisx in 1:nx
                for iisy in 1:ny
                    for iisz in 1:nz
                        w = wfuncs_all[iw, iisx, iisy, iisz, 1]
                        t = abs2(w)
                        if t > tmaxx
                            tmaxx = t
                            cmod = w
                        end
                    end
                end
            end
            cmod /= abs(cmod)
            @views wfuncs_all[:, :, :, iw, 1] ./= cmod
        end
    end
    str_cell = ustrip.(chk_info.cell)
    points = [str_cell * Point3((x-1)/nrx, (y-1)/nry, (z-1)/nrz) for x in supercell_xrange, y in supercell_yrange, z in supercell_zrange]
    if size(tu,5) == 1
        wfuncs_out = Vector{WannierFunction{1, eltype(wfuncs_all).parameters[1]}}(undef, size(wfuncs_all, 1))
        Threads.@threads for i=1:size(wfuncs_all, 1)
            wfuncs_out[i] = WannierFunction{1, eltype(wfuncs_all).parameters[1]}(points, map(x -> SVector(x), view(wfuncs_all,i, :, :, :, 1)))
        end
        return normalize!.(wfuncs_out)
    else
        wfuncs_out = Vector{WannierFunction{2, eltype(wfuncs_all).parameters[1]}}(undef, size(wfuncs_all, 1))
        Threads.@threads for i=1:size(wfuncs_all, 1)
            wfuncs_out[i] = WannierFunction{2, eltype(wfuncs_all).parameters[1]}(points, map(x -> SVector(x), zip(view(wfuncs_all, i, :, :, :, 1), view(wfuncs_all, i, :, :, :, 2))))
        end

        return wfuncs_out
    end
end

function generate_wannierfunctions(job::Job, supercell::NTuple{3,Int}, args...)
    if ismagnetic(job.structure) && Structures.iscolin(job.structure) && !any(Calculations.issoccalc, job.calculations)
        wfuncs = Vector{WannierFunction}[]
        for (is, s) in enumerate(("up", "down"))
            wan_calc  = getfirst(x -> eltype(x)==Wannier90&& x[:spin] == s, job.calculations)
            chk_info  = read_chk(joinpath(job, "$(wan_calc.name).chk"))
            unk_files = filter(x->occursin(".$is", x), searchdir(job, "UNK"))
            push!(wfuncs, generate_wannierfunctions(unk_files, chk_info, supercell, args...))
        end
        return (up=wfuncs[1], down=wfuncs[2])
    else
        wan_calc  = getfirst(x -> eltype(x)==Wannier90, job.calculations)
        chk_info  = read_chk(joinpath(job, "$(wan_calc.name).chk"))
        unk_files = searchdir(job, "UNK")
        return generate_wannierfunctions(unk_files, chk_info, supercell, args...)
    end
end
function bloch_sum(wfunc, kpoint)
    cell_boundaries = div.(size(wfunc.points), 3)
    bloch = WannierFunction(wfunc.points, copy(wfunc.values))
    dims = size(wfunc.values)
    for i1 in -1:1:1, i2 in -1:1:1, i3 in -1:1:1
        R_cryst = Vec3(i1, i2, i3)
        o1, o2, o3 = cell_boundaries .* R_cryst
        phase = ℯ^(2im*π*(R_cryst ⋅ kpoint))
        if i1 + i2 + i3 == 0
            continue
        end
        for j3 in 1:dims[3]
            oid3 = j3 - o3
            if !(0 < oid3 <= dims[3])
                continue
            end
            for j2 in 1:dims[2]
                oid2 = j2 - o2
                if !(0 < oid2 <= dims[2])
                    continue
                end
                for j1 in 1:dims[1]
                    oid1 = j1 - o1
                    if !(0 < oid1 <= dims[1])
                        continue
                    end
                    bloch.values[j1, j2, j3] += phase * wfunc.values[oid1, oid2, oid3]
                end
            end
        end
    end
    return bloch
end


text:
using StaticArrays
using Test

using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.Checkpoint
using ClimateMachine.ConfigTypes
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using ClimateMachine.Grids
using ClimateMachine.ODESolvers

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

Base.@kwdef struct AcousticWaveSetup{FT}
    domain_height::FT = 10e3
    T_ref::FT = 300
    α::FT = 3
    γ::FT = 100
    nv::Int = 1
end

function (setup::AcousticWaveSetup)(problem, bl, state, aux, localgeo, t)
    # callable to set initial conditions
    FT = eltype(state)

    λ = longitude(bl, aux)
    φ = latitude(bl, aux)
    z = altitude(bl, aux)

    β = min(FT(1), setup.α * acos(cos(φ) * cos(λ)))
    f = (1 + cos(FT(π) * β)) / 2
    g = sin(setup.nv * FT(π) * z / setup.domain_height)
    Δp = setup.γ * f * g
    p = aux.ref_state.p + Δp

    ts = PhaseDry_pT(bl.param_set, p, setup.T_ref)
    q_pt = PhasePartition(ts)
    e_pot = gravitational_potential(bl.orientation, aux)
    e_int = internal_energy(ts)

    state.ρ = air_density(ts)
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.ρe = state.ρ * (e_int + e_pot)
    return nothing
end

function main()
    FT = Float64

    # DG polynomial orders
    N = (4, 4)

    # Domain resolution
    nelem_horz = 4
    nelem_vert = 6
    resolution = (nelem_horz, nelem_vert)

    t0 = FT(0)
    timeend = FT(3600)
    # Timestep size (s)
    dt = FT(1800)

    setup = AcousticWaveSetup{FT}()
    T_profile = IsothermalProfile(param_set, setup.T_ref)
    orientation = SphericalOrientation()
    ref_state = HydrostaticState(T_profile)
    turbulence = ConstantDynamicViscosity(FT(0))
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        param_set;
        init_state_prognostic = setup,
        orientation = orientation,
        ref_state = ref_state,
        turbulence = turbulence,
        moisture = DryModel(),
        source = (Gravity(),),
    )

    ode_solver = ClimateMachine.MultirateSolverType(
        splitting_type = ClimateMachine.HEVISplitting(),
        fast_model = AtmosAcousticGravityLinearModel,
        implicit_solver_adjustable = true,
        slow_method = LSRK54CarpenterKennedy,
        fast_method = ARK2ImplicitExplicitMidpoint,
        timestep_ratio = 300,
    )
    driver_config = ClimateMachine.AtmosGCMConfiguration(
        "GCM Driver test",
        N,
        resolution,
        setup.domain_height,
        param_set,
        setup;
        solver_type = ode_solver,
        model = model,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )

    cb_test = 0
    result = ClimateMachine.invoke!(
        solver_config;
        user_info_callback = () -> cb_test += 1,
    )
    @test cb_test > 0
end

main()

ext: jl

max_forks_repo_name: jm-c/CLIMA


text:
function FromVk(def::StructDefinition{false})
    p = Dict(
        :category => :function,
        :name => :from_vk,
        :args => [:(T::Type{$(name(def))}), :(x::$(def.spec.name))],
        :short => true,
    )
    p[:body] = :(T($(filter(!isnothing, from_vk_call.(def.spec.members))...)))
    FromVk(def, p)
end


text:
"""
path related functions.
"""
module PATH
using Libdl
using Pkg
using ..Comonicon

function project(m::Module, xs...)
    path = pathof(m)
    path === nothing && return dirname(Pkg.project().path)
    return joinpath(dirname(dirname(path)), xs...)
end

project(xs...) = project(Comonicon, xs...)

sysimg() = "libcomonicon.$(Libdl.dlext)"
sysimg(name) = "lib$name.$(Libdl.dlext)"

"""
    default_exename()

Default Julia executable name: `joinpath(Sys.BINDIR, Base.julia_exename())`
"""
default_exename() = joinpath(Sys.BINDIR, Base.julia_exename())

"""
    default_julia_bin()

Return the default path to `.julia/bin`.
"""
default_julia_bin() = joinpath(first(DEPOT_PATH), "bin")

"""
    default_julia_fpath()

Return the default path to `.julia/completions`
"""
default_julia_fpath() = joinpath(first(DEPOT_PATH), "completions")

"""
    default_name(x)

Return the lowercase of `nameof(x)` in `String`.
"""
default_name(x) = lowercase(string(nameof(x)))

"""
    default_name(x)

Return the lowercase of a given package name. It will
ignore the suffix if it ends with ".jl".
"""
function default_name(x::String)
    if endswith(x, ".jl")
        name = x[1:end-3]
    else
        name = x
    end
    return lowercase(name)
end

end


text:
module PostProcess

using DataFrames
using JLD2
using PlotlyJS
using HypothesisTests
using Clustering

function ma(X, k)
    pad = zeros(eltype(X), k ÷ 2)
    return [pad; [mean(X[i:(i + k)]) for i in 1:(length(X) - k)]; pad]
end

function bname(fpath)
    return split(fpath, r"[\\/]")[end]
end

"""
    metadata(file::AbstractString)

Get metadata from file
"""
function metadata(file::AbstractString)
    return jldopen(f -> f["metadata"], file, "r")
end

function merge_result(file::AbstractString)
    io = jldopen(file, "r")
    df = DataFrame(; step=0:2880)
    foreach(keys(io)) do k
        if all(isdigit.(collect(k)))
            seed = parse(Int, k)
            df_seed = io[k]
            for name in names(df_seed)
                if name === "step"
                    continue
                end
                df[!, "$(name)_$(string(seed; pad=4))"] = df_seed[!, name]
            end
        end
    end
    columns = ["step"; sort(setdiff(names(df), ["step"]))]
    close(io)
    return select!(df, columns)
end

function PlotlyJS.plot(filepath::AbstractString, attr)
    f = jldopen(filepath, "r")
    traces = map(1:1000) do seed
        key = string(seed)
        df = f[key]
        scatter(; x=df.step, y=df[:, attr], name="Seed $key")
    end
    filebasename = splitpath(filepath)[end]
    title = replace(filebasename, "_" => "\n")
    layout = Layout(; title="$title", xaxis_title="Step", yaxis_title="$attr")
    close(f)
    return plot(traces, layout)
end

function plot_bph(filepath)
    return plot(filepath, :count_bph)
end
function plot_rice(filepath)
    return plot(filepath, :food)
end

function test_rice(filepath, p0)
    f = jldopen(filepath, "r")
    passed = map(1:1000) do seed
        key = string(seed)
        df = f[key]
        df.food[end] < df.food[begin] ÷ 2 - 100
    end
    close(f)
    return t = BinomialTest(count(passed), length(passed), p0)
end

function is_flower_effective(t; alpha=0.05)
    # H0: flower is not effective -> p = p0
    # Ha: flower is effective -> p < p0 -> tail = left
    # good -> reject H0 -> alpha > pvalue
    #a, b = confint(t; tail=:right, level=1-alpha)
    return alpha ≥ pvalue(t; tail=:left)
end

function test_rice(files::AbstractVector, p0; alpha=0.05)
    metadata = map(files) do file
        jldopen(f -> f["metadata"], file)
    end
    df = DataFrame(metadata)
    df.test = map(files) do file
        test_rice(file, p0)
    end
    df = transform(
        df,
        :envmap => ByRow(bname) => :envmap,
        :test => ByRow(t -> t.x) => :k,
        :test => ByRow(t -> t.n) => :n,
        :test => ByRow(t -> is_flower_effective(t; alpha=alpha)) => :accept,
    )
    return df
end

function batch_test_rice(dir::AbstractString, p0; alpha, clean=false)
    files = joinpath.(dir, readdir(dir))
    files = filter(files) do f
        isfile(f) && endswith(f, ".jld2")
    end
    df = test_rice(files, p0; alpha=alpha)
    if clean
        select(df, Not("test"))
    else
        df
    end
end

export JLD2, PlotlyJS, ProgressMeter, HypothesisTest, plot, savefig, test_rice

# end module
end
	smooth = 4
	# a2 = Axis(f[1,2])
	n = length(X)
	lines!(a1, X)

	lines!(a1, ma(X, smooth))
	# lines!(a1, fill(mean(df.count_bph), n))
	# lines!(a1, fill(quantile(df.count_bph, 0.25), n))
	# lines!(a1, fill(quantile(df.count_bph, 0.75), n))
	lines!(a1, fill(thres * std(X) + mean(X), n))
	peaks = peak_population(X, threshold=thres, smooth=smooth)
	scatter!(a1, peaks, X[peaks])
	lines!(a1, peaks, X[peaks])
	f
end

# ╔═╡ f8950add-031c-4af6-af56-be09d60f3f8f
let f = Figure()
	ax = Axis(f[1, 1])
	jldopen(file) do f
		map(1:1000) do sd
			df = f["$sd"]
			lines!(ax, df.step, df.count_bph)
		end
	end
	f
end

# ╔═╡ 6e98acef-274a-4e7d-8d69-96dbc116c333
let f = Figure()
	ax = Axis(f[1, 1])
	jldopen(file) do f
		map(1:1000) do sd
			df = f["$sd"]
			lines!(ax, df.step/24, df.count_bph)
		end
	end
	f
end

# ╔═╡ 9ea40c9e-d563-454e-987a-a17574620bfa
import PlotlyJS

# ╔═╡ 9c81ba94-5993-47ec-a57c-11ab6bb1c89c
dinhray = let f = jldopen(file)
	traces = map(1:1000) do s
		df = f["$s"]
		PlotlyJS.scatter(; x = df.step/24, y = df.count_bph, legend=false, name="Seed $s")
	end
	close(f)
	p = PlotlyJS.plot(traces, PlotlyJS.Layout(showlegend=false))
	PlotlyJS.savefig(p, "/tmp/dinh-ray.png")
end

# ╔═╡ b042ad4f-48e2-480b-93ba-ae6b4c94b6d2
LocalResource(dinhray)

# ╔═╡ 2820e1eb-80bb-4a1e-9d80-2bdd00f8bea7
DownloadButton(read(dinhray), splitpath(dinhray)[end])

# ╔═╡ de405835-0d61-4b68-bf94-ab7d796a69c6
md"# Batch peak data"

# ╔═╡ 68e35a9d-ab6e-420f-9c96-866b7fab75a8
function file2peaks(file)
	peaks = jldopen(file) do f
		map(1:1000) do seed
			df = f["$seed"]
			peak_population(df.count_bph)
		end
	end
	peaks = let lens = length.(peaks)
		s = std(lens)
		m = mean(lens)
		a, b = m - 3s, m + 3s
		filter(x -> a ≤ length(x) ≤ b, peaks)
	end
	npeaks = maximum(length.(peaks))
	peaks = reduce(vcat, peaks)

	cl = kmeans(transpose(peaks), npeaks)
	peaks = map(1:npeaks) do k
		peaks[k .== cl.assignments]
	end
	mpeaks = mean.(peaks)
	perm = sortperm(mpeaks)
	peaks[perm]
end

# ╔═╡ 0a70ede0-63b8-44cf-953e-cf081d7e5759
file2peaks(file)

# ╔═╡ 3005e698-d920-4767-bbcb-420c2739f591
md"# Peak across data files"

# ╔═╡ eb459b23-2f8b-4d9a-a349-5ddee45cc9bf
function parsefile(path)
	meta = jldopen(f -> f["metadata"], path)
	name = last(split(meta.envmap, r"[\\/]"))
	(meta..., envmap=name, path=path)
end

# ╔═╡ d8b94bd2-e0e7-4da3-9582-726f28269232
peakdf = let df = DataFrame(parsefile.(files))
	df
	df = transform(df, :path => ByRow(file2peaks) => :peaks)
end

# ╔═╡ 6cd527cd-0eb7-43c2-98a2-d2164e7d8187
function splitpeaks(df)
	function getstat(stat, x, n)
		stat(get(x, n, [missing]))
	end
	npeaks = maximum(length.(df.peaks))
	transformation = map(1:npeaks) do n
		[:peaks => ByRow(x -> getstat(mean, x, n)) => Symbol("peak_$n"),
		:peaks => ByRow(x -> getstat(std, x, n)) => Symbol("std_peak_$n")]
	end
	df = transform(df, transformation...)
	df = select(df, Not([:peaks, :path]))
end

# ╔═╡ 3e1f5482-76e3-4ce8-a211-d50c482e0394
statpeak = splitpeaks(peakdf)

# ╔═╡ 72b5b845-c7e2-4028-844b-deea0aec22ca
cleanpeak = let to_day = x -> @. round(x / 24, digits=1)
	cols = names(statpeak)
	cols = filter(x -> occursin("peak", x), names(statpeak))
	tfm = map(cols) do col
		col => to_day => col
	end
	transform(statpeak, tfm...)
end

# ╔═╡ 2b188539-5a31-4b89-84e9-382b95218a80
cleandf2 = let df = cleanpeak
	df = @where(df, (:pr_killed .== 0.15) .| (:pr_killed  .== 0))
	df = select(df,  [[3, 1, 4, 2]; 5:size(df, 2)])
	df = sort(df, [:pr_killed, :envmap, :init_position, :nb_bph_init], rev=true)
	df = transform(df, :envmap => ByRow(x -> match(r"(\dx\d)", x)[1]) => :envmap)
	cols = names(df)
	cols = map(enumerate(cols)) do (i, name)
		if i == 2
			raw"$L$"
		elseif i == 1
			raw"$\Pr_{\text{kill}}$"
		elseif i == 3
			raw"$P_0$"
		elseif i == 4
			raw"$N$"
		else
			name = replace(raw"$" * name * raw"$", "std" => raw"\sigma")
			name = replace(name, r"sigma_peak_(\d+)" => s"sigma(t_\1)")
			name = replace(name, raw"$sigma" => raw"$\sigma")
			name = replace(name, r"peak_" => raw"t_")
		end
	end
	df = rename(df, cols)
	n = size(df, 2)
	df = select(df, Not([n, n-1]))
	df
end

# ╔═╡ 9dad2144-45c7-46f2-959f-c42093aeb63f
md"# Latex"

# ╔═╡ 12c6567b-4029-463e-8b80-a8c0cb0a64d5
function formaltable(df)
	s = latextabular(df; latex=false)
	lines = split(s, "\n")
	tbegin = replace(lines[1], "" => "")
	# header
	thead = strip(replace(lines[2], "\\\\" => ""))
	# thead = join("\\textbf{" .* split(thead, r" *& *") .* "}", " & ")
	tbody = join(lines[3:end-2], "\n\\midrule\n")
	s = """
	$(tbegin)
	\\toprule
	$(thead) \\\\
	\\midrule
	\\midrule
	$(tbody)
	\\bottomrule
	\\end{tabular}
	"""
	s = replace(s, r"([^a-zA-Z_.])([0-9]+\.?[0-9]*)([^a-zA-Z_.])" => s"\1$\2$\3")
	# s = replace(s, r"([0-9]+\.[0-9]+)" => s"$\1$")
	s = replace(s, "1x2" => raw"\OneByTwo")
	s = replace(s, "3x3" => raw"\ThreeByThree")
	s = replace(s, "false" => raw"\Reject")
	s = replace(s, "true" => raw"\Accept")
	s = replace(s, "corner" => raw"\Corner")
	s = replace(s, "border" => raw"\Border")
	# s = replace(s, r"_" => raw"\_")
	s = replace(s, "#" => raw"\#")
	s = replace(s, "missing" => "--")
	clipboard(s)
	return s
end

# ╔═╡ 5cf687c3-7d2f-4751-af5a-9c38ebb2d4f1
s = formaltable(cleandf2)

# ╔═╡ b619d9cb-64d2-4e88-9421-7a202d71f691
with_terminal(() -> println(s))

# ╔═╡ 6847a684-f74c-49df-baf3-c71ba4e102f4
write("population-peak.tex", s)

# ╔═╡ Cell order:
# ╠═ed7b7fa2-ca12-11eb-3142-830f24a474d2
# ╟─20e4a147-4f6f-4c21-97f4-a62cc6d24745
# ╠═9023aa1e-509a-4ca9-b35b-2d2875dda584
# ╠═087175fd-bc09-4623-91a5-b779c0fc62c8
# ╠═28d52568-c5de-4d0c-b2fc-a1e3a3690042
# ╠═22f38eaa-8e10-4dcd-a469-635cbb390480
# ╠═a3ab9801-22f8-4ee1-8be6-c33acf9007bb
# ╠═86e45cf0-0290-4350-bfdb-056073630021
# ╠═e18bc7ed-15c3-4e1d-a90b-75fc2367ae1b
# ╠═fbfd993c-3897-4c5c-b504-d8ace62762de
# ╟─c9e69c07-2b3b-40c5-9aaa-88407040e882
# ╠═d9bcb11b-61c0-48d0-a90a-72e0e766c22e
# ╠═5639f9fe-660c-430f-a7d2-3d1aad4d6152
# ╠═62586022-2d8a-43ef-98f0-a039d446cd25
# ╠═6aa63ead-af1e-4e1d-b3e7-60d3a0f700d4
# ╠═f8950add-031c-4af6-af56-be09d60f3f8f
# ╠═6e98acef-274a-4e7d-8d69-96dbc116c333
# ╠═9ea40c9e-d563-454e-987a-a17574620bfa
# ╠═9c81ba94-5993-47ec-a57c-11ab6bb1c89c
# ╠═b042ad4f-48e2-480b-93ba-ae6b4c94b6d2
# ╠═2820e1eb-80bb-4a1e-9d80-2bdd00f8bea7
# ╟─de405835-0d61-4b68-bf94-ab7d796a69c6
# ╠═b3ed17c4-213a-4547-b4be-2fd43612b563
# ╠═68e35a9d-ab6e-420f-9c96-866b7fab75a8
# ╠═0a70ede0-63b8-44cf-953e-cf081d7e5759
# ╟─3005e698-d920-4767-bbcb-420c2739f591
# ╠═eb459b23-2f8b-4d9a-a349-5ddee45cc9bf
# ╠═d8b94bd2-e0e7-4da3-9582-726f28269232
# ╠═6cd527cd-0eb7-43c2-98a2-d2164e7d8187
# ╠═3e1f5482-76e3-4ce8-a211-d50c482e0394
# ╠═7280ad9b-5a95-44c0-9a48-2b90195c176f
# ╠═72b5b845-c7e2-4028-844b-deea0aec22ca
# ╠═ee13068b-f94a-4c4d-9749-faa27e4c1d4e
# ╠═2b188539-5a31-4b89-84e9-382b95218a80
# ╠═9dad2144-45c7-46f2-959f-c42093aeb63f
# ╠═5b8f13a4-fbc8-470e-99d6-33cbbfdb9030
# ╠═12c6567b-4029-463e-8b80-a8c0cb0a64d5
# ╠═5cf687c3-7d2f-4751-af5a-9c38ebb2d4f1
# ╠═b619d9cb-64d2-4e88-9421-7a202d71f691
# ╠═6847a684-f74c-49df-baf3-c71ba4e102f4

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
using Test
using MiscFunctions

@testset "Tests" begin

    x = [-8π, -2π, 2π, 0, 4π, 6π]
    m = cmean(x)
    sd = cstd(x)
    @test abs(m) < 1e-15
    @test sd ≈ 0

    array = [0.1, -8.6, 2.2, 1.97, 0.15, -1.93, -0.1]
    value = 1.94
    index = find_closest(value, array)
    @test index == 4

    value = -2.2
    index = find_closest(value, array)
    @test index == 6

    # Test binary boundary search
    X_target = -13
    fun = x -> x > X_target
    err = 0.1
    x_range = [-17, 84]
    X = binary_boundary_search(fun, x_range, err)
    @test abs(X_target - X) < err

    # Target is outside of the range to the left
    fun = x -> x > -100
    X = binary_boundary_search(fun, x_range, err)
    @test X == x_range[1]

    # Target is outside of the range to the right
    fun = x -> x > 100
    X = binary_boundary_search(fun, x_range, err)
    @test X == x_range[2]

    # r2_score
    x = [1, 2, 3]
    y = [3, 5, 1]
    R2 = r2_score(x, y)
    @test R2 ≈ -7.5

    x = [0, 4, 5]
    y = [1, 2, 3]
    R2 = r2_score(x, y)
    @test R2 ≈ 1-9/14

    x = [0, 1, -8]
    y = [0, 1, -8]
    R2 = r2_score(x, y)
    @test R2 ≈ 1

    # ptp
    x = [1 5 -1 7; 4 0 -5 0]
    @test ptp(x) ≈ 12
    @test ptp(x; dims=2) ≈ [8; 9]

    # smooth
    y = [-2, -5, 1, -1, 3, 0, 6]
    yy = smooth(y)
    @test all(yy .≈ [-2, -2, -0.8, -0.4, 1.8, 3, 6])

    yy = smooth(y, 1)
    @test all(yy .≈ [-2, -5, 1, -1, 3, 0, 6])

    yy = smooth(y, 100)
    @test all(yy .≈ [-2, -2, -0.8, 0.2857142857142857, 1.8, 3, 6])

    # find_all_combinations
    vectors = [[1, 0, 2]]
    combinations = find_all_combinations(vectors)
    @test combinations == hcat([1, 0, 2])

    vectors = [[1, 0, 2], [7, 6]]
    combinations = find_all_combinations(vectors)
    @test combinations == [1 7; 1 6; 0 7; 0 6; 2 7; 2 6]

    vectors = [[1, 0, 2], ['a'], [7, 6]]
    combinations = find_all_combinations(vectors)
    @test combinations == [1 'a' 7; 1 'a' 6; 0 'a' 7; 0 'a' 6; 2 'a' 7; 2 'a' 6]

end

max_stars_repo_name: vkumpost/MiscFunctions

max_issues_repo_name: vkumpost/MiscFunctions

max_forks_repo_name: vkumpost/MiscFunctions


text:
using Documenter
using Mercury

makedocs(
    sitename = "Mercury.jl",
    format = Documenter.HTML(prettyurls = false),
    pages = ["Introduction" => "index.md", "API" => "api.md"],
)

deploydocs(repo = "github.com/RoboticExplorationLab/Mercury.jl.git", devbranch = "main")


text:
"""
    Coulomb Potential
"""
mutable struct Coulomb <: Potential
    q_1::Float64
    q_2::Float64
    ε0::Float64
end

"""
    Coulomb(params::Dict)

Creates a Coulomb potential.
"""
function Coulomb(params::Dict)
    # Creates the Coulomb model
    q_1 = params["q_1"]
    q_2 = params["q_2"]
    ε0 = params["ε0"]
    return Coulomb(q_1, q_2, ε0)
end

"""
    potential_energy(p::Coulomb, r::Position, args...)

Calculates Coulomb potential energy.
"""
function potential_energy(p::Coulomb, r::Position, args...)
    return p.q_1 * p.q_2 / (4.0 * π * p.ε0 * norm(r))
end


text:
"""
    function negloglik(ts::TS; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand())::Real

Compute the negative log-likelihood of ts with parameters in kwargs.
"""
function negloglik(ts::TS; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand())::Real
    tst = [μ, ϕ, γ] .< 0
    any(tst) && (return Inf)

    times = first(ts.times) isa TimeType ? Dates.value.(DateTime.(ts.times)) ./ (1000*3600*24) : ts.times
    starttime = first(times)
    endtime = last(times)
    T = endtime - starttime

    vol = volfunc(times, ts, γ)
    term1 = sum(log.(μ .+ ϕ .* vol))
    term2 = μ * T + ϕ/γ * sum(1 .- exp.(-γ .* (endtime .- times)))

    return term2 - term1
end

"""
    negloglik(sepp::SEPPExpKern)::Real

Compute the negative log-likelihood of the time series in data with parmaters of the sepp model.
"""
function negloglik(sepp::SEPPExpKern)::Real
    ts = sepp.data
    isnothing(ts) && error("No data in model, can't compute log-likelihood")

    θ = params(sepp)
    return negloglik(ts ; θ...)
end

#=
function nlL1(mts::MarkedTimeSeries; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand(), δ::Real = 0)
    tst = [μ, ϕ, γ, δ] .< 0
    w = [:μ, :ϕ, :γ, :δ][tst]
    any(tst) && ((μ, ϕ, γ, δ) = abs.((μ, ϕ, γ, δ)) ; @warn string(string(["$symb " for symb in w]...), "must be positive or zero, taking absolute value"))

    times = first(mts.times) isa TimeType ? Dates.value.(Date.(mts.times)) : mts.times
    starttime = first(times)
    endtime = last(times)
    T = endtime - starttime

    marks = mts.marks

    vol = volfunc(times, mts, γ, δ)

    term1 = sum(log.(μ .+ ϕ .* vol))
    term2 = μ * T + ϕ/γ * sum((1 .+ δ .* marks) .* (1 .- exp.(-γ .* (endtime .- times))))

    return term2 - term1
end
=#

"""
    negloglik(mts::MarkedTimeSeries, markdens ; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand(), δ::Real = 0, ξ::Real = rand(), α::Real = rand(), β::Real = rand(), κ::Real = 1)::Real

Compute the negative log-likelihood of the marked time series with parmaters in kwargs.
"""
function negloglik(mts::MarkedTimeSeries, markdens::SupportedMarksDistributions, impact_func::Function = (d -> 1 + tanh(d)) ; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand(), δ::Real = 0, ξ::Real = rand(), α::Real = rand(), β::Real = rand(), κ::Real = 1)::Real
    tst = [μ, ϕ, γ, δ, β, α, κ] .< 0
    any(tst) && (return Inf)

    times = first(mts.times) isa TimeType ? Dates.value.(DateTime.(mts.times)) ./ (1000*3600*24) : mts.times
    starttime = first(times)
    endtime = last(times)
    T = endtime - starttime

    marks = mts.marks

    vol = volfunc(times, mts, γ, δ, impact_func)

    term1 = sum(log.(μ .+ ϕ .* vol))
    term2 = μ * T + ϕ/γ * sum((impact_func(δ .* marks)) .* (1 .- exp.(-γ .* (endtime .- times))))

    σ = β .+ α .* vol
    sig_marks = hcat(σ, marks)

    if markdens == Distributions.GeneralizedPareto
        mark_contrib = (sig_mark -> logpdf(markdens(0, sig_mark[1], ξ), sig_mark[2])).(eachrow(sig_marks))
    else        # EGPD case
        mark_contrib = (sig_mark -> logpdf(markdens(sig_mark[1], ξ, κ), sig_mark[2])).(eachrow(sig_marks))
    end

    term3 = sum(mark_contrib)

    return term2 - term1 - term3

end

"""
    negloglik(sempp::SEMPPExpKern)::Real

 Compute the negative log-likelihood of the marked time series in data with parmaters of the sempp model.
"""
function negloglik(sempp::SEMPPExpKern)::Real
    mts = sempp.data
    isnothing(mts) && error("No data in model, can't compute log-likelihood")

    θ = params(sempp)
    markdens = θ[:markdens]
    impact_func = θ[:impact_function]
    delete!(θ, :markdens)
    delete!(θ, :impact_function)
    return negloglik(mts, markdens, impact_func ; θ...)
end

"""
    fit!(sepp::SEPPExpKern)::Real

Fit a self exciting point process model (whithout marks) to the time series in data, based on MLE.

Note that if the time series has marks, this method ignores them, that is modelling the ground process as independent of the marks.
Returns the minimal Log-likelihood found.
"""
function fit! end

include("macro_def_fits.jl")

@def_fit_sepp [:s1]
@def_fit_sepp [:s1, :s2]
@def_fit_sepp [:s1, :s2, :s3]

function fit!(sepp::SEPPExpKern)
    return fit!(sepp, :μ, :ϕ, :γ)
end

@def_fit_sempp [:s1]
@def_fit_sempp [:s1, :s2]
@def_fit_sempp [:s1, :s2, :s3]
@def_fit_sempp [:s1, :s2, :s3, :s4]
@def_fit_sempp [:s1, :s2, :s3, :s4, :s5]
@def_fit_sempp [:s1, :s2, :s3, :s4, :s5, :s6]
@def_fit_sempp [:s1, :s2, :s3, :s4, :s5, :s6, :s7]
@def_fit_sempp [:s1, :s2, :s3, :s4, :s5, :s6, :s7, :s8]

function fit!(sempp::SEMPPExpKern)
    return fit!(sempp, :μ, :ϕ, :γ, :δ, :ξ, :α, :β)
end

#=
function fit!(sepp::SEPPExpKern)::Real
    ts = sepp.data
    isnothing(ts) && error("No data in model, can't fit")

    model = Model(Ipopt.Optimizer)
    θ = params(sepp)

    function to_min(μ, ϕ, γ)
        return negloglik(ts, μ = μ, ϕ = ϕ, γ = γ)
    end

    JuMP.register(model, :to_min, 3, to_min, autodiff=true)

    @variables(model, begin
        mu >= 0, (start = θ[:μ])
        phi >= 0, (start = θ[:ϕ])
        gamma >= 0, (start = θ[:γ])
    end)
    @NLobjective(model, Min, to_min(mu, phi, gamma))

    optimize!(model)

    sepp.μ = value(mu)
    sepp.ϕ = value(phi)
    sepp.γ = value(gamma)

    x = [sepp.μ, sepp.ϕ, sepp.γ]

    to_hess(θ) = to_min(θ...)

    hess = ForwardDiff.hessian(to_hess, x)

    sepp.cov_mat = inv(hess)

    return objective_value(model)
end

"""
    fit!(sempp::SEMPPExpKern, bounds::Union{Vector{<:Real}, Nothing} = nothing)::Real

Fit a self exciting marked point process model to the marked time series in data, based on MLE.

Returns the minimal Log-likelihood found.
"""
function fit!(sempp::SEMPPExpKern, bounds::Union{Vector{<:Real}, Nothing} = nothing)::Real # default xi >= 0

    mts = sempp.data
    isnothing(mts) && error("No data in model, can't fit")

    markdens = sempp.markdens

    if markdens == Distributions.GeneralizedPareto
        df = DataFrame(Marks = mts.marks)
        gml = gpfit(df, :Marks)
        sempp.β = exp(gml.θ̂[1])
        sempp.ξ = gml.θ̂[2]
    else
        egppower = EGPD.EGPpowerfit(mts.marks)
        sempp.ξ = egppower.ξ
        sempp.β = egppower.σ
        sempp.κ = egppower.κ
    end

    model = Model(Ipopt.Optimizer)
    θ = params(sempp)

    function to_min_GPD(μ, ϕ, γ, δ, ξ, α, β)
        return negloglik(mts, Distributions.GeneralizedPareto, μ = μ, ϕ = ϕ, γ = γ, δ = δ, ξ = ξ, α = α, β = β)
    end

    function to_min_EGPD(μ, ϕ, γ, δ, ξ, α, β, κ)
        return negloglik(mts, markdens, μ = μ, ϕ = ϕ, γ = γ, δ = δ, ξ = ξ, α = α, β = β, κ = κ)
    end

    if markdens == Distributions.GeneralizedPareto
        JuMP.register(model, :to_min_GPD, 7, to_min_GPD, autodiff=true)
    else
        JuMP.register(model, :to_min_EGPD, 8, to_min_EGPD, autodiff=true)
    end

    @variables(model, begin
        mu >= 0, (start = θ[:μ])
        phi >= 0, (start = θ[:ϕ])
        gamma >= 0, (start = θ[:γ])
        delta >= 0, (start = θ[:δ])
        alpha >= 0, (start = θ[:α])
        beta >= 0, (start = θ[:β])
    end)

    if markdens == EGPD.EGPpower
        @variable(model, kappa >= 0, start = θ[:κ])
    end

    if isnothing(bounds)
        @variable(model, xi >= 0, start = θ[:ξ])
    else
        @variable(model, bounds[1] <= xi <= bounds[2], start = θ[:ξ])
    end

    if markdens == Distributions.GeneralizedPareto
        @NLobjective(model, Min, to_min_GPD(mu, phi, gamma, delta, xi, alpha, beta))
    else
        @NLobjective(model, Min, to_min_EGPD(mu, phi, gamma, delta, xi, alpha, beta, kappa))
    end

    optimize!(model)

    sempp.μ = value(mu)
    sempp.ϕ = value(phi)
    sempp.γ = value(gamma)
    sempp.δ = value(delta)
    sempp.ξ = value(xi)
    sempp.α = value(alpha)
    sempp.β = value(beta)

    x = [sempp.μ, sempp.ϕ, sempp.γ, sempp.δ, sempp.ξ, sempp.α, sempp.β]

    if markdens == Distributions.GeneralizedPareto
        to_hess1(θ) = to_min_GPD(θ...)
        hess = ForwardDiff.hessian(to_hess1, x)
    else
        sempp.κ = value(kappa)
        push!(x, sempp.κ)
        to_hess2(θ) = to_min_EGPD(θ...)
        hess = ForwardDiff.hessian(to_hess2, x)
    end

    sepp.cov_mat = inv(hess)

    return objective_value(model)
end
=#


text:
module ReactiveMPMessageTest

using Test
using ReactiveMP
using Distributions
using Random

import InteractiveUtils: methodswith
import Base: methods
import Base.Iterators: repeated, product
import ReactiveMP: materialize!
import ReactiveMP: mirrorlog, xtlog
import SpecialFunctions: loggamma

@testset "Message" begin

    @testset "Default methods" begin
        data    = PointMass(1)

        for clamped in [ true, false ], initial in [ true, false ]
            msg = Message(data, clamped, initial)
            @test getdata(msg)    === data
            @test is_clamped(msg) === clamped
            @test is_initial(msg) === initial
            @test materialize!(msg) === msg
            @test occursin("Message", repr(msg))
        end

    end

    @testset "multiply_messages" begin
        dist1 = NormalMeanVariance(randn(), rand())
        dist2 = NormalMeanVariance(randn(), rand())

        @test getdata(Message(dist1, false, false) * Message(dist2, false, false))    == prod(ProdAnalytical(), dist1, dist2)
        @test getdata(Message(dist2, false, false) * Message(dist1, false, false))    == prod(ProdAnalytical(), dist2, dist1)

        for (left_is_initial, right_is_initial) in product(repeated([ true, false ], 2)...)
            @test is_clamped(Message(dist1, true, left_is_initial) * Message(dist2, false, right_is_initial)) == false
            @test is_clamped(Message(dist1, false, left_is_initial) * Message(dist2, true, right_is_initial)) == false
            @test is_clamped(Message(dist1, true, left_is_initial) * Message(dist2, true, right_is_initial))  == true
            @test is_clamped(Message(dist2, true, left_is_initial) * Message(dist1, false, right_is_initial)) == false
            @test is_clamped(Message(dist2, false, left_is_initial) * Message(dist1, true, right_is_initial)) == false
            @test is_clamped(Message(dist2, true, left_is_initial) * Message(dist1, true, right_is_initial))  == true
        end

        for (left_is_clamped, right_is_clamped) in product(repeated([ true, false ], 2)...)
            @test is_initial(Message(dist1, left_is_clamped, true) * Message(dist2, right_is_clamped, true)) == !(left_is_clamped && right_is_clamped)
            @test is_initial(Message(dist2, left_is_clamped, true) * Message(dist1, right_is_clamped, true)) == !(left_is_clamped && right_is_clamped)
            @test is_initial(Message(dist1, left_is_clamped, false) * Message(dist2, right_is_clamped, false)) == false
            @test is_initial(Message(dist2, left_is_clamped, false) * Message(dist1, right_is_clamped, false)) == false
        end

        @test is_initial(Message(dist1, true, true) * Message(dist2, true, true)) == false
        @test is_initial(Message(dist1, true, true) * Message(dist2, true, false)) == false
        @test is_initial(Message(dist1, true, false) * Message(dist2, true, true)) == false
        @test is_initial(Message(dist1, false, true) * Message(dist2, true, false)) == true
        @test is_initial(Message(dist1, true, false) * Message(dist2, false, true)) == true
        @test is_initial(Message(dist2, true, true) * Message(dist1, true, true)) == false
        @test is_initial(Message(dist2, true, true) * Message(dist1, true, false)) == false
        @test is_initial(Message(dist2, true, false) * Message(dist1, true, true)) == false
        @test is_initial(Message(dist2, false, true) * Message(dist1, true, false)) == true
        @test is_initial(Message(dist2, true, false) * Message(dist1, false, true)) == true
    end

    @testset "Statistics" begin

        distributions = [
            PointMass(0.5),
            Gamma(10.0, 2.0),
            NormalMeanVariance(-10.0, 10.0),
            Wishart(4.0, [ 2.0 -0.5; -0.5 1.0 ]),
            MvNormalMeanPrecision([ 2.0, -1.0 ], [ 7.0 -1.0; -1.0 3.0 ]),
            Bernoulli(0.5),
            Categorical([ 0.8, 0.2 ])
        ]

        # Here we get all methods defined for a particular type of a distribution
        dists_methods = map(d -> methodswith(eval(nameof(typeof(d)))), distributions)

        methods_to_test = [
            Distributions.mean,
            Distributions.median,
            Distributions.mode,
            Distributions.shape,
            Distributions.scale,
            Distributions.rate,
            Distributions.var,
            Distributions.std,
            Distributions.cov,
            Distributions.invcov,
            Distributions.logdetcov,
            Distributions.entropy,
            Distributions.params,
            Base.precision,
            Base.length,
            Base.ndims,
            Base.size,
            mean_cov,
            mean_invcov,
            mean_precision,
            weightedmean_cov,
            weightedmean_invcov,
            weightedmean_precision,
            probvec,
            weightedmean
        ]

        for (distribution, distribution_methods) in zip(distributions, dists_methods), method in methods_to_test
            T       = typeof(distribution)
            message = Message(distribution, false, false)
            # Here we check that a specialised method for a particular type T exist
            ms = methods(method, (T, ))
            if !isempty(ms) && all(m -> m ∈ distribution_methods, ms)
                @test method(message) == method(distribution)
            end
        end

        fn_mean_functions = (inv, log, xtlog, mirrorlog, loggamma)

        for distribution in distributions, fn_mean in fn_mean_functions
            F       = typeof(fn_mean)
            T       = typeof(distribution)
            message = Message(distribution, false, false)
            # Here we check that a specialised method for a particular type T exist
            ms = methods(mean, (F, T, ), ReactiveMP)
            if !isempty(ms)
                @test mean(fn_mean, message) == mean(fn_mean, distribution)
            end
        end

        _getpoint(rng, distritubution) = _getpoint(rng, variate_form(distritubution), distritubution)
        _getpoint(rng, ::Type{ <: Univariate }, distribution) = 10rand(rng)
        _getpoint(rng, ::Type{ <: Multivariate }, distribution) = 10 .* rand(rng, 2)

        distributions2   = [
            Gamma(10.0, 2.0),
            NormalMeanVariance(-10.0, 1.0),
            MvNormalMeanPrecision([ 2.0, -1.0 ], [ 7.0 -1.0; -1.0 3.0 ]),
            Bernoulli(0.5),
            Categorical([ 0.8, 0.2 ])
        ]

        methods_to_test2 = [
            Distributions.pdf,
            Distributions.logpdf,
        ]

        rng = MersenneTwister(1234)

        for distribution in distributions2, method in methods_to_test2
            message = Message(distribution, false, false)

            for _ in 1:3
                point = _getpoint(rng, distribution)
                @test method(message, point) === method(distribution, point)
            end
        end

    end

end

end


text:
using BigArrays.GSDicts
using Test

# test storage utility first
#include(joinpath(dirname(@__FILE__), "google_cloud/storage_util.jl"))

@testset "\ntest gsdict io" begin
    kv  = GSDict("gs://seunglab/jpwu/test/image")

    println("test fakekey error capture")
    try
        kv["fakekey"]
    catch err
        if isa(err, KeyError)
            println("get normal KeyError: $(err)")
        else
            println("get an error: ", err)
            println("error type: ", typeof(err))
            rethrow()
        end
    end

    a = rand(UInt8, 50)

    kv["test"] = a
    b = kv["test"]
    b = reinterpret(UInt8, b)

    println("make sure that the value saved in the cloud is the same with local")
    @test all(a .== b)
    @test true == haskey(kv, "test")
    @info("delete the file in google cloud storage")
    delete!(kv, "test")

    println("test no such key error...")
    try
        kv["test"]
    catch err
        @test isa(err, KeyError)
    end

end


text: ≈(x, y; tol=0.01) = isapprox(x, y, atol=tol)


text:
function kernel(i)
    getfield((1,2,3), i, false)
end

code_llvm(stdout, kernel, Tuple{Int})

using TinyCassette
struct Ctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{Ctx,
                typeof(kernel),
                Int})


text:
# mps.jl
"""
    abstract type AbstractMPS{S<:EuclideanSpace}

The abstract type of all 1D tensor networks, includes 1D MPS and 1D MPO.
"""
abstract type AbstractMPS{S<:EuclideanSpace} end

"""
    field(::Type{<:AbstractMPS{S}}) where {S<:EuclideanSpace}
    field(mps::AbstractMPS)

Return the field type over which an MPS instance or type is defined.
"""
TX.field(::Type{<:AbstractMPS{S}}) where {S<:EuclideanSpace} = field(S)
TX.field(mps::AbstractMPS) = field(typeof(mps))

"""
    spacetype(::Type{<:AbstractMPS{S}}) where {S<:EuclideanSpace}
    spacetype(mps::AbstractMPS)

Return the ElementarySpace type associated to an MPS instance or type.
"""
TX.spacetype(::Type{<:AbstractMPS{S}}) where {S<:EuclideanSpace} = S
TX.spacetype(mps::AbstractMPS) = spacetype(typeof(mps))

"""
    sectortype(::Type{<:AbstractMPS{S}}) where {S<:EuclideanSpace}
    sectortype(mps::AbstractMPS) -> Type{<:Sector}

Return the type of sector over which the MPS instance or type is defined.
"""
TX.sectortype(::Type{<:AbstractMPS{S}}) where {S<:EuclideanSpace} = sectortype(S)
TX.sectortype(mps::AbstractMPS) = sectortype(typeof(mps))

"""
    get_chi(mps::AbstractMPS)

Return the bond dimensions of the MPS on bonds 1:N+1.
"""
function get_chi(mps::AbstractMPS)
    chi_list = append!([dim(codomain(mps.Ts[i],1)) for i in 1:length(mps.Ts)],
                        [dim(domain(mps.Ts[end]))])
    return chi_list
end

"""
    entanglement_entropy(mps::AbstractMPS)
"""
function entanglement_entropy(mps::AbstractMPS)
    S_list = [ ]
    for i in 1:length(mps.Ss)
        entropy = 0.0
        for c in blocksectors(mps.Ss[i])
            t = block(mps.Ss[i],c)
            entropy += -sum([(t[j,j])^2*log((t[j,j])^2) for j in 1:size(t,1)])
        end
        push!(S_list,entropy)
    end
    return S_list
end

"""
    expectation_value_single_site(mps::AbstractMPS)
"""
function expectation_value_single_site(mps::AbstractMPS)

end

# Base.getindex(mps::AbstractMPS, i::Integer) = mps.sitetensors[i]
#
# @inline function Base.setindex!(mps::AbstractMPS, v, indices::Vararg{Int})
#     data = t[]
#     @boundscheck checkbounds(data, indices...)
#     @inbounds data[indices...] = v
#     return v
#
# @inline function Base.iterate(psi::AbstractMPS, ::Val{i} = Val(1)) where {i}
#     if i > length(psi)
#         return nothing
#     else
#         return psi.sitetensors[i], Val(i+1)
#     end
# end
#
# Base.eltype(psi::AbstractMPS) = eltype(psi.sitetensors[1])
#
# Base.IteratorEltype(::Type{<:AbstractMPS}) = Base.HasEltype()
# Base.IteratorSize(::Type{<:AbstractMPS}) = Base.HasLength()
#
# function Base.hash(psi::AbstractMPS, h::UInt)
#     h = hash(psi.canonical_form,h)
#     h = hash(psi.schmidtvalues, h)
#     h = hash(psi.sitetensors, h)
#     return h
# end
#
# function Base.:(==)(psi1::AbstractMPS, psi2::AbstractMPS)
#     psi1.canonical_form == psi.canonical_form || return false
#     psi1.schmidtvalues == psi.schmidtvalues || return false
#     psi1.sitetensors == psi2.sitetensors || return false
#     return true
# end

# function Base.isapprox(psi1::AbstractMPS, psi2::AbstractMPS;
#                 atol::Real=0, rtol::Real=Base.rtoldefault(eltype(psi1), eltype(psi2), atol))
#     d = mapreduce(norm,+,ntuple(x -> psi1[x]-psi2[x],length(psi1)))
#     if isfinite(d)
#         return d <= max(atol, rtol*max(norm(t1), norm(t2)))
#     else
#         return false
#     end
# end


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

import Base: size, length

function NURBS(i, p, u, t)
    p == 0 && return t[i] <= u <= t[i+1] ? 1.0 : 0.0
    anom = u-t[i]
    adenom = t[i+p]-t[i]
    a = isapprox(adenom, 0.0) ? 0.0 : anom/adenom
    bnom = t[i+p+1]-u
    bdenom = t[i+p+1]-t[i+1]
    b = isapprox(bdenom, 0.0) ? 0.0 : bnom/bdenom
    result = a*NURBS(i,p-1,u,t) + b*NURBS(i+1,p-1,u,t)
    return result
end


text:
# # quadratic.jl
# # Source: PolynomialRoots.jl: solve_quadratic_eq

function solve_quadratic_roots(polycoeff::Vector{T}) where {T<:Real}

    Rts = Vector{Complex{T}}()

    if length(polycoeff) < 3
        @error "There are too few coefficients given."
        return Rts
    end

    (c, b, a) = polycoeff

    if isapprox(a,zero(T))
        @warn "The leading cubic coefficient is approximately zero. Returning empty array."
        return Rts
    end

    # This is copied from PolynomialRoots.jl:
    Δ = sqrt(Complex(b*b - 4*a*c))
    if real(conj(b)*Δ) >= 0
        x0 = -(b + Δ)/2
    else
        x0 = -(b - Δ)/2
    end
    if x0 == 0
        x1 = x0
    else
        x1 = c / x0 # Viete's formula
        x0 = x0 / a
    end
    # end copied section.
    push!(Rts,x0,x1)
    return Rts
end

max_stars_repo_name: jd-foster/HyperCubicRoots

max_issues_repo_name: jd-foster/HyperCubicRoots

max_forks_repo_name: jd-foster/HyperCubicRoots


text:
/**
 * This class models Crypto operations.
 */
class CryptoLibrary {
    //static int cl1TESTOUTPUT = 7;

    // All the interesting constraints are coded up in CEDataFlow
    public static int encrypt(int publicKey, int plaintext) {
 	// here is where the magic encrypting happens.
	// the result has no relation to the plaintext.
	return  ({L})plaintext+1;
    }

    public static int decrypt(int privateKey, int ciphertext) {
	return ciphertext-1;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
# https://github.com/docker-library/julia/pull/6
download("https://google.com")

# https://github.com/docker-library/julia/pull/9
if VERSION.major > 0 || (VERSION.major == 0 && VERSION.minor >= 7)
	# https://github.com/docker-library/julia/pull/21
	# https://github.com/JuliaLang/julia/tree/v0.7.0-beta2/stdlib/Pkg
	using Pkg
end
Pkg.add("JSON")

max_stars_repo_name: sortie/official-images

max_issues_repo_name: sortie/official-images

max_forks_repo_name: sortie/official-images


text:
using BenchmarkTools
using Random
using Test

@testset "algorithmic order" begin
    ##
    ## The goal of this problem is to get you thinking about algorithmic order.
    ##
    # There's no src/ code for this problem, write everything here.
    # Note that `@elapsed` is like time and is simple in that it just runs your operation once.
    # You may need to execute this file more than once in case compilation time complicates the interpretation.

    str = randstring(10^6)
    d = Dict{Char,Int}()
    t = @elapsed for c in ['a':'z'; 'A':'Z'; '0':'9']
        d[c] = count(==(c), str)
    end

    d2 = Dict{Char,Int}()
    t2 = @elapsed for c in str
        d2[c] = get(d2, c, 0) + 1
    end
    @test t2 < t

    # Since it's a discrete list of values, can you do even better with an array?
    # Hint: try `Int('z')` for ideas.
    counter = zeros(Int, 26*2+10)

    t3 = @elapsed for c in str
        if Int(c) > 96
            idx = (c-'a')+1
        elseif Int(c) > 64
            idx = (c-'A') + 27
        else
            idx = (c-'0') + 53
        end
        counter[idx]+=1
    end
    @test t3 < t2
end


text:
using MegTools
using AxisKeys
using Test

# BESA epoch data
test_path = joinpath(@__DIR__, "testData", "test_data.mat")
test_data = load_cont_epochs(test_path)

# BESA averaged data
BESA_av = joinpath(@__DIR__, "testData", "besa_av_singlecondition.mat")
BESA_av_data = load_besa_av(BESA_av)
# Brainstorm epoch data_path
BStest_path = joinpath(@__DIR__, "testData", "sub_001")
BStest_data = load_BSepochs(BStest_path)

@testset "IO_Meg.jl" begin
    # Testing BESA output
    @test test_data isa Dict
    @test length(collect(keys(test_data))) == 4
    for condition in collect(keys(test_data))
        @test size(test_data[condition])==(1501, 319, 100)  # Time, Channels, Trials
    end
    # Loading average
    @test BESA_av_data isa KeyedArray

    # Testing Brainstorm outputs
    @test BStest_data isa Dict
    @test length(collect(keys(BStest_data))) == 2 # There are only 2 conditions
    for condition in collect(keys(BStest_data))
        # In condition "7" one trial is rejected (005- as labelled in brainstormstudy.mat)
        @test size(BStest_data[condition])==(1501, 351, 2)  # Time, Channels, Trials
    end

end

@testset "analysis_recipes.jl" begin
    # Testing single conditon, average_across_trials
    @test ndims(test_data["1"]) == 3
    averaged_trials = average_across_trials(test_data["1"])
    @test ndims(averaged_trials) == 2
    @test size(test_data["1"])[1:2] == size(averaged_trials)
    # Subject (Dict) test
    subject_av = average_across_trials(test_data)
    @test  subject_av isa Dict
    @test ndims(subject_av["1"]) == 2
    @test size(test_data["1"])[1:2] == size(subject_av["1"])

    # Testing single condition, select_channels
    auditory_channels,left_labels,right_labels = select_channels(
        test_data["1"],
        paradigm="auditoryN1m",
    )
    @test size(auditory_channels)[2] ==  (length(left_labels) + length(right_labels))
    @test auditory_channels.channels == vcat(left_labels, right_labels)
    # Selecting custom channels
    custom_channels, custom_left, custom_right = select_channels(
        test_data["1"],
        paradigm="custom_channels",
        left_channels=["MEG0241", "MEG0231"],
        right_channels="MEG2641"
    )
    @test length(custom_channels.channels) == 3
    @test length([custom_left;custom_right]) == 3
        # Selecting custom channels
        custom_channels, custom_left, custom_right = select_channels(
            test_data["1"],
            paradigm="custom_channels",
            left_channels=["MEG0241", "MEG0231"],
            right_channels="MEG2641"
        )
        @test length(custom_channels.channels) == 3
        @test length([custom_left;custom_right]) == 3
    # Select from subject (Dict)
    auditory_channels,left_labels,right_labels = select_channels(
        test_data,
        paradigm="auditoryN1m",
    )
     # Select INCORRECT paradigm from subject (Dict)
    @test_throws MethodError auditory_channels,left_labels,right_labels = select_channels(
        test_data,
        paradigm="HamSandwich",
    )
    @test size(auditory_channels["1"])[2] ==  (length(left_labels) + length(right_labels))
    @test auditory_channels["1"].channels == vcat(left_labels, right_labels)
    # Same with averaged trial, which is of a single condition
    auditory_channels,left_labels,right_labels = select_channels(
        averaged_trials,
        paradigm="auditoryN1m",
    )
    @test size(auditory_channels)[2] ==  (length(left_labels) + length(right_labels))
    @test auditory_channels.channels == vcat(left_labels, right_labels)
    # Testing baseline_correction
    filtered = highlow_butterworth_filter(auditory_channels, 1000)
    baseline_corrected = baseline_correction(filtered)
    @test filtered ≠ baseline_corrected
    @test sum(baseline_corrected) < sum(filtered)
    # Single subject (Dict)
    baseline_corrected = baseline_correction(test_data)
    @test test_data["1"] ≠ baseline_corrected["1"]
    @test sum(baseline_corrected["1"]) < sum(test_data["1"])
    # Getting the baseline value (using averaged data (all trials))
    baseline_corrected, indv_baseline = baseline_correction(averaged_trials, output_baseline = true)
    @test length(indv_baseline) == length(averaged_trials.channels)
    # Now checking a dict of subjects
    baseline_corrected, indv_baseline = baseline_correction(test_data, output_baseline = true)
    @test size(indv_baseline["1"]) == (1,319,100) # (time, channels, trials)
    # Determine if robust with different number of total trials
    not_filtered =  test_data["1"]
    baseline_corrected = baseline_correction(not_filtered)
    not_filtered_1 =  test_data["1"][:,:,1]
    baseline_corrected_1 = baseline_correction(not_filtered_1)
    @test baseline_corrected[:,:,1] == baseline_corrected_1
    # Testing with input baseline_condiitons (averaging approach as in Zacharias 2011)
    baseline_candidates= ["4","1"]
    averaged_baseline = get_averaged_baseline(test_data, baseline_candidates)
    @test size(averaged_baseline)[1] == 1
    @test size(averaged_baseline)[3] == 1
    @test size(averaged_baseline)[2] == length(test_data["1"].channels)
    # Carrying out the baseline correction using the averaging method
    # Checking that the baselines are different
    av_baseline_corrected = baseline_correction(test_data, averaged_baseline)
    per_erf_baseline = baseline_correction(test_data)
    @test av_baseline_corrected ≠ per_erf_baseline
    # However, they are should be the same shape
    @test length(av_baseline_corrected) == length(per_erf_baseline)
    @test size(av_baseline_corrected["1"]) == size(per_erf_baseline["1"])
    # Testing find find_peaks
    a,_,_,b,_,_,c,d = find_peaks(averaged_trials, left_labels, right_labels)
    @test ndims(a) == 1
    @test length(a) == length(averaged_trials[:,1])
    @test ndims(b) == 1
    @test length(b) == length(averaged_trials[:,1])
    # Testing custom channel with find peaks

    custom_channels2, custom_left2, custom_right2 = select_channels(
        averaged_trials,
        paradigm="custom_channels",
        left_channels=["MEG0241", ],
        right_channels=["MEG2641", ]
    )
    a,_,_,b,_,_,c,d = find_peaks(custom_channels2, custom_left2, custom_right2)
    # Making sure we only have a single channel in the output
    @test size(a)[2] == 1
    @test length(a) == length(averaged_trials[:,1])
    @test size(b)[2] == 1
    @test length(b) == length(averaged_trials[:,1])
 # Testing with singular non-array input
    custom_channels3, custom_left3, custom_right3 = select_channels(
        averaged_trials,
        paradigm="custom_channels",
        left_channels="MEG0241",
        right_channels="MEG2641"
    )
    a,_,_,b,_,_,c,d = find_peaks(custom_channels3, custom_left3, custom_right3)
    @test ndims(a) == 1
    @test length(a) == length(averaged_trials[:,1])
    @test ndims(b) == 1
    @test length(b) == length(averaged_trials[:,1])

    # Testing find_peaks with incorrect dimentions
    @test_throws MethodError a,_,_,b,_,_,c,d = find_peaks(test_data, left_labels, right_labels)
    test_dict = Dict("fake_array"=>rand(10,10,1,10))
    @test_throws MethodError a,_,_,b,_,_,c,d = find_peaks(test_dict, left_labels, right_labels)
    # Testing mean amplitude of the m100 peaks
    _,l,_,r,_,_ = find_mean_amplitude(averaged_trials, left_labels, right_labels)
    @test length(l) == 1
    @test length(r) == 1
    #@test_throws MethodError find_mean_amplitude(test_data["1"], left_labels, right_labels)
    @test_logs (:error, "Dimentions more than 2, average this data first.")  find_mean_amplitude(test_data["1"], left_labels, right_labels)
    @test_throws MethodError find_mean_amplitude(test_dict, left_labels, right_labels)
    # Testing subject data will all conditions
    averaged_trials_all = average_across_trials(test_data)
    peak_means = find_mean_amplitude(averaged_trials_all, left_labels, right_labels)
    @test peak_means["1"]["left_mean_amplitude"] == l
    @test peak_means["1"]["right_mean_amplitude"] == r
    # Testing collection of mean amplitudes
    soi, left, right = collect_mean_amps(peak_means)
    @test length(soi) == length(left) == length(right)
    @test soi == sort(soi)
    # Single subject (Dict)
    subpeaks = find_peaks(subject_av, left_labels, right_labels)
    @test ndims(subpeaks["1"]["left_peak_erf"]) == 1
    @test length(subpeaks["1"]["left_peak_erf"]) == length(subject_av["1"][:,1])
    @test ndims(subpeaks["1"]["right_peak_erf"]) == 1
    @test length(subpeaks["1"]["right_peak_erf"]) == length(subject_av["1"][:,1])
    sois, left, right = collect_peaks(subpeaks)
    @test length(sois) ==  length(left) ==  length(right)
    @test typeof(sois[1]) == Float64

    # Testing will all averaged_trials
    #=TODO? a,b,c,d = find_peaks(
        highlow_butterworth_filter(test_data["1"], 1000),
        left_labels,
        right_labels
    )
    =#

end

@testset "filters.jl" begin
    # Testing single condition, highlow_butterworth_filter
    averaged_trials = average_across_trials(test_data["1"])
    auditory_channels,left_labels,right_labels = select_channels(
        averaged_trials,
        paradigm="auditoryN1m",
    )
    filtered = highlow_butterworth_filter(auditory_channels, 1000)
    filtered_nooffset = highlow_butterworth_filter(auditory_channels, 1000, offset=false)
    @test filtered ≠ auditory_channels
    @test filtered ≠ filtered_nooffset
    @test size(filtered) == size(auditory_channels)
    @test filtered.channels == vcat(left_labels, right_labels)
    # Subject (Dict)
    subject_av = average_across_trials(test_data)
    auditory_channels,left_labels,right_labels = select_channels(
        subject_av,
        paradigm="auditoryN1m",
    )
    filtered = highlow_butterworth_filter(auditory_channels, 1000)
    for (cond, data) in filtered
        @test filtered[cond]  ≠ auditory_channels[cond] && size(filtered[cond]) == size(auditory_channels[cond])
        @test filtered[cond].channels == vcat(left_labels, right_labels)
    end
    filtered_nooffset = highlow_butterworth_filter(auditory_channels, 1000, offset=false)
    for (cond, data) in filtered
        @test filtered_nooffset[cond]  ≠ filtered[cond] && size(filtered_nooffset[cond]) == size(filtered[cond])
        @test filtered_nooffset[cond].channels == vcat(left_labels, right_labels)
    end

end

@testset "trigger_keys.jl" begin
    triggers = load_trigger_values("regsoi")
    @test triggers isa Dict
    @test length(triggers) == 20
    try
        load_trigger_values("MyUnknownProject")
    catch
        @test 1==1 #statating that an error is occuring
    end
    triggers = load_trigger_values("soistream")
    @test triggers isa Dict
    @test length(triggers) == 10
end

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text:
# addprocs(3)
using Test
# using Compat
# using IncrementalInference

@testset "out of module evalPotential..." begin
    include("TestModuleFunctions.jl")
end

include("testStateMachine.jl")

include("testCompareVariablesFactors.jl")

@testset "Ensure memory return is working properly..." begin
    include("typeReturnMemRef.jl")
end

include("basicGraphsOperations.jl")

include("testPartialFactors.jl")

@testset "basic Bayes tree construction" begin
    include("testBayesTreeiSAM2Example.jl")
end
@testset "Ensure converter types can be run from extending namespaces..." begin
    include("saveconvertertypes.jl")
end
@testset "packing converters work..." begin
    include("testgraphpackingconverters.jl")
end
include("testNLsolve.jl")

@testset "generic root finding by numeric solve of residual functions..." begin
    include("testNumericRootGenericRandomized.jl")
end
@testset "GenericWrapParam functors..." begin
    include("testCommonConvWrapper.jl")
end
@testset "with simple local constraint examples Odo, Obsv2..." begin
    include("testlocalconstraintexamples.jl")
end

include("testFactorMetadata.jl")

include("testExplicitMultihypo.jl")

include("TestCSMMultihypo.jl")

include("testMultiHypo2Door.jl")

include("testMultithreaded.jl")

@testset "partial constraints..." begin
    include("testpartialconstraint.jl")
end
@testset "null hypothesis..." begin
    include("testnullhypothesis.jl")
end
@testset "standardized multihypothesis..." begin
    include("testmultihypothesisapi.jl")
end
@testset "with local Graphs.jl dictionary and arrays only (multicore)..." begin
    include("fourdoortest.jl")
end

@testset "saving to and loading from FileDFG" begin
    saveFolder = "/tmp/dfg_test"
    saveDFG(fg, saveFolder)
    retDFG = GraphsDFG{SolverParams}(params=SolverParams())
    retDFG = loadDFG(saveFolder, IncrementalInference, retDFG)
    @test symdiff(ls(fg), ls(retDFG)) == []
    @test symdiff(lsf(fg), lsf(retDFG)) == []
end

@warn "must return testExpandedJLD.jl to testing -- currently skipped since jld2 files cannot be loaded."
# include("testExpandedJLD.jl")

include("testTexTreeIllustration.jl")

#


text:
module NonuniformResampling1D

include("WindowFunction.jl")
include("window_functions.jl")
include("range_utilities.jl")

export nuresample

"""
    nuresample(xin, yin, xout, smoothing_function; kwargs...)

Take data `yin` sampled at uniformly spaced locations `xin` and resample at
nonuniform locations `xout`.

# Positional Arguments

- `xin` is an `AbstractRange` of input locations
- `yin` are the input values
- `xout` are the increasing output locations
- `smoothing_function` is a [`WindowFunction`](@ref). The width of the window
  function is scaled proportional to the distance to the nearest neighbor on the
  left and right side of the output point. At the start and end of `xout` a
  symmetrical window is used. The default is a rectangular window without
  overlap.

# Keyword Arguments

- `upsampling_function`: a [`WindowFunction`](@ref) that is used to perform
  upsampling when the density of input points is not high enough compared to the
  density of output points. The default is Lanczos3.
- `required_points_per_slice`: the number of input points that must be in a
  slice (left or right half of a window) before smoothing. If the number of
  input points is lower than this value, upsampling is performed to create
  `required_points_per_slice` many points before applying the smoothing
  function.
  The default is `round(4 * smoothing_function.width)`, which means 4 points are
  required per unit width to the nearest neighbor.
  The minimum value is 1, in which case upsampling is only performed if no input
  point falls in a slice.

# Examples

```jldoctest
julia> nuresample(1:9, 1:9, [4.2, 6.2])
2-element Vector{Float64}:
 4.201106036510037
 6.201106036510037
```

```jldoctest
julia> nuresample(1:9, 1:9, [4.2, 6.2],
        rect_window(.5), # moving average without overlap
        required_points_per_slice = 1, # no upsampling
        upsampling_function = lanczos_window(2)) # Lanczos2 interpolation if a slice were empty
2-element Vector{Float64}:
 4.5
 6.5
```

# Output Type

Returns an `Array{Float64, 1}`. Other output types are currently unsupported.
"""
function nuresample(xin::AbstractRange, yin, xout, smoothing_function=rect_window(); kwargs...)
    nuresample(StepRangeLen(xin), yin, xout, smoothing_function; kwargs...)
end

function nuresample(xin::StepRangeLen, yin, xout,
    smoothing_function::WindowFunction = rect_window();
    required_points_per_slice::Integer=Int(round(4 * smoothing_function.width)),
    upsampling_function::WindowFunction=lanczos_window()
    )
    # validate inputs
    @assert required_points_per_slice >= 1 "required_points_per_slice must at least be 1"
    @assert smoothing_function.width > 0 "smoothing_function must have a width bigger than 0. "*
        "Maybe use upsampling_function as smoothing_function."
    @assert upsampling_function.width >= 1 "upsampling_function width must at least be 1 to ensure "*
        "there's at least one input point to each side of the upsampled point"
    # allocate
    yout = Array{Float64, 1}(undef, length(xout))

    # first unit
    out_ind = 1
    right_unit_width = calculate_right_unit_width(xout, out_ind)
    left_unit_width = right_unit_width

    while true
        yout[out_ind] = interpolate_point(xin, yin, xout[out_ind], left_unit_width, right_unit_width,
        smoothing_function,
        required_points_per_slice=required_points_per_slice,
        upsampling_function=upsampling_function)

        # prepare next point
        left_unit_width = right_unit_width # next left unit width is the last right unit width
        out_ind += 1
        out_ind > length(xout) && break

        if out_ind < length(xout) # keep previous slice width on last element
            right_unit_width = calculate_right_unit_width(xout, out_ind)
        end
    end

    yout
end

function calculate_right_unit_width(xout, out_ind)
    right_unit_width = xout[out_ind + 1] - xout[out_ind]
    right_unit_width > 0 || throw(ArgumentError(
        "xout must be increasing everywhere. Violated between index $out_ind "*
        "and $(out_ind + 1): $(xout[out_ind]) >= $(xout[out_ind+1])"
    ))
    right_unit_width
end

struct Slice
        width
        start
        stop
        start_ind
        stop_ind
        upsampling_required
end

function Slice(unit_width, window::WindowFunction, xpoint, xin, left::Bool, required_points_per_slice)
    # the unit width is the width to the neighboring point
    # now calculate slice width, which is the width in which points will be
    # considered according to the window funciton width
    width = unit_width*window.width
    if left
        start = xpoint - width
        stop = xpoint
    else
        start = xpoint
        stop = xpoint + width
    end
    @assert start > xin[1] - Float64(xin.step) "not enough points at the beginning of the input for slice from $start to $stop"
    @assert stop < xin[end] + Float64(xin.step) "not enough points at the end of the input for slice from $start to $stop"

    #find relevant input indices
    start_ind = find_first_above_or_equal(start, xin)
    stop_ind = if left find_last_below(stop, xin) else find_last_below_or_equal(stop, xin) end

    input_points = stop_ind - start_ind + 1
    upsampling_required = input_points < required_points_per_slice

    Slice(
        width,
        start,
        stop,
        start_ind,
        stop_ind,
        upsampling_required,
    )
end

function interpolate_point(xin, yin, xpoint, left_unit_width, right_unit_width, window::WindowFunction;
    required_points_per_slice=1, upsampling_function=missing)

    left = Slice(left_unit_width, window, xpoint, xin, true, required_points_per_slice)
    right = Slice(right_unit_width, window, xpoint, xin, false, required_points_per_slice)

    if left.upsampling_required || right.upsampling_required
        # upsample both slices, even if only required for one of them, to keep balance between left and right
        # while still maintaining continuity when an output point is moved slightly over an input point
        @assert !ismissing(upsampling_function) "interpolation required but no upsampling_function given"
        upsample_step = min(left.width, right.width)/required_points_per_slice
        left_x, left_y = upsample_prepare_input(left, xin, yin, upsample_step, upsampling_function)
        right_x, right_y = upsample_prepare_input(right, xin, yin, upsample_step, upsampling_function)
    else
        # neither of the slices needs upsampling
        left_x, left_y = prepare_input(left, xin, yin)
        right_x, right_y = prepare_input(right, xin, yin)
    end

    left_val_acc, left_win_acc = slice_weighted_mean(xpoint, left_unit_width, left_x, left_y, window)
    right_val_acc, right_win_acc = slice_weighted_mean(xpoint, right_unit_width, right_x, right_y, window)

    return (left_val_acc + right_val_acc)/(left_win_acc + right_win_acc)
end

function prepare_input(sc::Slice, xin, yin)
    input_x = xin[sc.start_ind:sc.stop_ind]
    input_y = @view yin[sc.start_ind:sc.stop_ind]
    return (input_x, input_y)
end

function upsample_prepare_input(slice::Slice, xin, yin, upsample_step, upsampling_function)
    input_x = range(slice.start + upsample_step/2, step=upsample_step, stop=slice.stop)
    in_step = Float64(xin.step)
    input_y = [interpolate_point(xin, yin, up_x, in_step, in_step, upsampling_function) for up_x in input_x]
    return (input_x, input_y)
end

function slice_weighted_mean(xpoint, unit_width, xin_points, yin_points, window)
    val_acc::Float64 = 0.
    win_acc::Float64 = 0.
    for i in eachindex(xin_points)
        x = xin_points[i]
        y = yin_points[i]
        rel_pos = abs(x - xpoint)/unit_width
        win_value::Float64 = window(rel_pos)
        win_acc += win_value
        val_acc += y*win_value
    end
    @assert win_acc > 0 "Slice weight should be bigger zero"
    return (val_acc, win_acc)
end

end # module


text:
# The Hilbert space spanned by the spherical harmonics obey the
# angular momentum algebra, and contractions can be expressed as
# matrix elements of spherical tensors, which can be calculated using
# Wigner 3-j symbols.

using WIGXJPF
using UnicodeFun

import Base: transpose, |, identity, float, ==, show

abstract AngularSpace <: HilbertSpace

type AngBra <: Bra
    j::Rational
    m::Rational
    AngBra(j,m) = abs(m) > j ? error("Invalid bra, |$(m)| > $(j)") : new(j, m)
end
bra(::Type{AngularSpace}) = AngBra
==(b1::Bra, b2::Bra) = (b1.j == b2.j) && (b1.m == b2.m)

type AngKet <: Ket
    j::Rational
    m::Rational
    AngKet(j,m) = abs(m) > j ? error("Invalid ket, |$(m)| > $(j)") : new(j, m)
end
ket(::Type{AngularSpace}) = AngKet
==(k1::Ket, k2::Ket) = (k1.j == k2.j) && (k1.m == k2.m)

transpose(b::AngBra) = AngKet(b.j, b.m)
transpose(k::AngKet) = AngBra(k.j, k.m)

type SphericalTensor <: Operator
    k::Integer # Rank
    q::Integer
    SphericalTensor(k, q) = abs(q) > k ? error("Invalid tensor, |$(q)| > $(k) ") : new(k, q)
end
==(C1::SphericalTensor, C2::SphericalTensor) =
    (C1.k == C2.k) && (C1.q == C2.q)

op(::Union{Type{AngularSpace},AngBra,AngKet}) = SphericalTensor
identity(::Union{Type{AngularSpace},SphericalTensor,Type{SphericalTensor}}) =
    SphericalTensor(0, 0)
space(::Union{AngBra,AngKet,SphericalTensor}) = AngularSpace

type AngScalar <: Scalar
    b::AngBra
    C::SphericalTensor
    k::AngKet
end

scalar(::Union{Type{AngularSpace},AngBra,AngKet}) = AngScalar

function float(s::AngScalar)
    full = wig3j(s.b.j, s.C.k, s.k.j,
                 -s.b.m, s.C.q, s.k.m)
    red = wig3j(s.b.j, s.C.k, s.k.j,
                0, 0, 0)
    (-1)^(2s.b.j-s.b.m)*√((2s.b.j+1)*(2s.k.j+1))*full*red
end

const ells = "spdfghiklmnopqrtuvwxyz"

function partial_wave_str(ell, m)
    rat(v) = den(v) == 1 || num(v) == 0 ?
        string(num(v)) : to_fraction(num(v), den(v))

    if ell == 0
        "s"
    elseif ell == 1
        if m == 0
            "p"
        elseif m == 1
            "p⁺"
        else
            "p⁻"
        end
    elseif den(ell) == 1 && den(m) == 1
        @sprintf("%s;%i", ells[num(ell)+1], m)
    else
        @sprintf("%s;%s", rat(ell), rat(m))
    end
end

show(stream::IO, b::AngBra) =
    write(stream::IO, @sprintf("〈%s|", partial_wave_str(b.j, b.m)))
show(stream::IO, k::AngKet) =
    write(stream::IO, @sprintf("|%s〉", partial_wave_str(k.j, k.m)))
show(stream::IO, C::SphericalTensor) =
    write(stream::IO, @sprintf("C%s%s",
                               to_subscript(C.q),
                               to_superscript("($(C.k))")))

function show(stream::IO, Ck::Tuple{SphericalTensor, AngKet})
    show(stream, Ck[1])
    show(stream, Ck[2])
end

function show(stream::IO, bC::Tuple{AngBra, SphericalTensor})
    show(stream, bC[1])
    show(stream, bC[2])
end

function show(stream::IO, s::AngScalar)
    show(stream, s.b)
    if s.C.k == 0 && s.C.q == 0
        write(stream, partial_wave_str(s.k.j, s.k.m))
        write(stream, "〉")
    else
        show(stream, s.C)
        show(stream, s.k)
    end
    write(stream::IO, @sprintf(" ≈ %0.5g", float(s)))
end

export AngularSpace


text:
## Model Generators for 2D Gaussian PSFs

using SpecialFunctions

"""
    PSF_airy2D <: PSF

Contains the parameter used to calculate an Airy Pattern PSF

The Airy PSF is

I(r)=ν²/(4π)(2*J₁(ν*r)/(ν*r))²

where

ν=πD/(λf)=2*π*nₐ/λ

!!! note
    The Gaussian approximation is σ = 0.42*π/ν
"""
struct PSF_airy2D <: PSF
    ν::Float32
end

# midpoint rule approximation to besselj1 integral
function mybesselj1(x::Float32, m::Int32=Int32(10))
    j1 = 0f0;
    for k = 1:(m - 1)
        j1 += sin(x * sin(pi / (2.0f0 * m) * (k + 0.5f0))) * sin(pi / (2.0f0 * m) * (k + .5f0))
    end
    return j1 / m
end

function airy_amplitude(r::Float32, ν::Float32)
    w = r * ν
    w = max(w, 1f-5)
    return (4π)^-(1 / 2) * ν * (2 * besselj1(w) / w)
end

function genmodel!(s::StateFlatBg, sz::Int32, psf::PSF_airy2D, model::Array{Float32,2})

    for jj = 1:sz
        for ii = 1:sz
            model[ii,jj] = s.bg
        end
    end

    for nn = 1:s.n
        y=s.y[nn]
        x=s.x[nn]
        photons=s.photons[nn];
         for jj = 1:sz
            for ii = 1:sz
                r = sqrt((ii - y)^2 + (jj - x)^2)
                model[ii,jj] += photons * airy_amplitude(r, psf.ν)^2
            end
        end
    end

end

function max2int(psf::PSF_airy2D)
    return 4π / psf.ν^2
end


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(The BTG Program)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#################################################################################
#copy this box into a file to use the BTG program
using LinearAlgebra
using Test
using Distributions
using DataFrames
using CSV
using Polynomials
using PyPlot
include("covariatefun.jl")
include("validation/det_loocv.jl")
include("validation/loocv.jl")
include("quadrature/quadrature.jl")
include("transforms/transforms.jl")
include("priors/priors.jl")
include("bayesopt/incremental.jl")
include("kernels/kernel.jl")
include("datastructs.jl")
include("options.jl")
include("computation/buffers0.jl") #datastruct, kernel, incremental, quadrature
include("model0.jl") #buffers, datastructs, several auxiliary
include("computation/finitedifference.jl")
include("computation/tdist.jl") #model0 and buffer0
include("statistics.jl") #posterior statistics computation
include("prediction.jl")
#################################################################################

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg