text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Xorg_xkbcomp"
version = v"1.4.2"

# Collection of sources required to build xkbcomp
sources = [
    ArchiveSource("https://www.x.org/archive/individual/app/xkbcomp-$(version).tar.bz2",
                  "6dd8bcb9be7e85bd7294abe261b8c7b0539d2fc93e41b80fb8bd013767ce8424"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/xkbcomp-*
./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p) || Sys.isfreebsd(p)]

products = [
    ExecutableProduct("xkbcomp", :xkbcomp)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    BuildDependency("Xorg_util_macros_jll"),
    Dependency("Xorg_libxkbfile_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: devmotion/Yggdrasil

max_issues_repo_name: devmotion/Yggdrasil

max_forks_repo_name: devmotion/Yggdrasil


text:
using JuliaRunClient
ctx = Context()
nb = self()

initParallel()

NWRKRS = 2
println("scale up to $NWRKRS")

@result setJobScale(ctx, nb, NWRKRS)
waitForWorkers(NWRKRS)

using StatsBase

@everywhere using StatsBase

@everywhere function montecarlo(howmany, data_generator, bins)
    h  = Histogram(bins)
    for i=1:howmany
        push!(h, data_generator() )
    end
    return h.weights
end

w = @parallel (+) for i=1:nworkers()
  montecarlo(100000, randn, -3:.1:3)
end;

using Plots

@everywhere  function tracywidom_sample(β=2,n=200)
  h=n^(-1/3)
  x=[0:h:10;]
  N=length(x)
  d=(-2/h^2 .-x) +  2/sqrt(h*β)*randn(N) # diagonal
  e=ones(N-1)/h^2                   # subdiagonal
  eigvals(SymTridiagonal(d,e))[N]
end

plot()
for β = [1,2,4]
 bins = -4:.05:0.95
 w=
  @parallel (+) for i=1:nworkers()
      montecarlo(10000,()->tracywidom_sample(β), -4:.05:1)
  end;
plot!(bins, w/sum(w)*bins.step.hi)
end

plot!()

# Scale down
@result setJobScale(ctx, self(), 0)

nworkers()

ext: jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Test, OrdinaryDiffEq, HDF5, ModelConstructors
include(joinpath(dirname(@__FILE__), "../../../src/includeall.jl"))

m = Li2020()
stategrid, _ = initialize_nojump!(m)
ode_f, ode_callback = eqcond_nojump(m)
tspan = (stategrid[:w][1], stategrid[:w][end])
θ = parameters_to_named_tuple(map(x -> m.parameters[m.keys[x]], get_setting(m, :nojump_parameters)))
prob = ODEProblem(ode_f, get_setting(m, :boundary_conditions)[:p][1], tspan, θ, tstops = stategrid[:w][2:end - 1])
sol = solve(prob, get_setting(m, :ode_integrator),
            reltol = get_setting(m, :ode_reltol),
            abstol = get_setting(m, :ode_abstol), callback = ode_callback)

true_soln = h5read(joinpath(dirname(@__FILE__), "../../reference/models/li2020/li2020_nojump_eqm.h5"), "p_DP5")
@testset "Baseline solution for no jump equilibrium" begin
    @test @test_matrix_approx_eq sol.u true_soln
end


text:
module Image

export Dimension, Pixel, PPM, saveimage, pixel, examplefilename
export dimension4K

using Dates
using CthulhuVision.Light

struct Dimension
    width::Int
    height::Int
end

dimension4K() = Dimension(3840, 2160)

struct Pixel
    x::Int
    y::Int
end

struct PPM
    dimension::Dimension
    pixels::Dict{Pixel, RGB}

    PPM(dimension::Dimension) = new(dimension, Dict{Pixel, RGB}())
end

pixel(image::PPM, px::Pixel, rgb::RGB) = image.pixels[px] = rgb

function saveimage(image::PPM, filename::String)
    open(filename, "w") do io
        write(io, "P3\n")
        write(io, "$(image.dimension.width) $(image.dimension.height)\n")
        write(io, "255\n")

        for y = image.dimension.height - 1:-1:0
            for x = 0:image.dimension.width - 1
                rgb = image.pixels[Pixel(x, y)]
                r = trunc(Int, 255.99f0*sqrt(rgb.r))
                g = trunc(Int, 255.99f0*sqrt(rgb.g))
                b = trunc(Int, 255.99f0*sqrt(rgb.b))
                write(io, "$r $g $b\n")
            end
        end
    end
end

function examplefilename(name::AbstractString) :: String
    date = Dates.format(now(), "yyyymmdd-HHMMSS")
    "$(name)-$(date).ppm"
end

end


text:
using KissCaches
using KissCaches: deletecached!
using Test

function self(args...; kw...)
    (args, kw)
end

function test_cached(cache, f, args...; kw...)
    @test !iscached(cache, f, args...; kw...)
    @cached cache f(args...; kw...)
    @test iscached(cache, f, args...; kw...)
    deletecached!(cache, f, args...; kw...)
    @test !iscached(cache, f, args...; kw...)
end

@testset "test" begin
    cache = Dict()
    @test !iscached(cache, sin, 1)
    @cached cache sin(1)
    @test iscached(cache, sin, 1)
    empty!(cache)
    @test !iscached(cache, sin, 1)
    @test !iscached(cache, cos, 1)

    test_cached(Dict(), sin, 1)
    test_cached(Dict(), self)
    test_cached(Dict(), self, 1)
    test_cached(Dict(), self, x=:x)
    test_cached(Dict(), self, 1, x=1)
    test_cached(Dict(), self, 1, x=1, y="hi", z=["a", Any])
    test_cached(Dict(), self, 1, x=1, y="hi", z=["a", Any])
end

@testset "GLOBAL_CACHE" begin
    @test !iscached(KissCaches.GLOBAL_CACHE, self)
    @cached self()
    @test iscached(KissCaches.GLOBAL_CACHE, self)
    for (args, kw) in [
            ((1,), NamedTuple()),
            ((), (a=1,)),
        ]
        @test !iscached(KissCaches.GLOBAL_CACHE, self, args...; kw...)
        @cached self(args...; kw...)
        @test iscached(KissCaches.GLOBAL_CACHE, self, args...; kw...)
    end
end

@testset "interesting behaviour" begin
    # This behaviour is not necesserily desired.
    # These tests are there so that we don't forget
    # about it

    cache = Dict()
    @cached cache self(1.0)
    @test iscached(cache, self, UInt(1))
    @test iscached(cache, self, Int(1))
    @test iscached(cache, self, Float64(1))
    @test iscached(cache, self, Bool(1))
    @test iscached(cache, self, Complex(1))
    @test !iscached(cache, self, fill(1, (1,1)))

    @cached cache self("hi")
    @test iscached(cache, self, "hi")
    @test !iscached(cache, self, :hi)

    @cached cache self(1:1)
    @test iscached(cache, self, 1:1)
    @test iscached(cache, self, [1])
    @test iscached(cache, self, [true])
    @test iscached(cache, self, [1.0])
end


text:
# ---
# title: 907. Sum of Subarray Minimums
# id: problem907
# author: Indigo
# date: 2021-07-12
# difficulty: Medium
# categories: Array, Stack
# link: <https://leetcode.com/problems/sum-of-subarray-minimums/description/>
# hidden: true
# ---
#
# Given an array of integers `A`, find the sum of `min(B)`, where `B` ranges
# over every (contiguous) subarray of `A`.
#
# Since the answer may be large, **return the answer modulo`10^9 + 7`.**
#
#
#
# **Example 1:**
#
#
#
#     Input: [3,1,2,4]
#     Output: 17
#     Explanation: Subarrays are [3], [1], [2], [4], [3,1], [1,2], [2,4], [3,1,2], [1,2,4], [3,1,2,4].
#     Minimums are 3, 1, 2, 4, 1, 1, 2, 1, 1, 1.  Sum is 17.
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 30000`
#   2. `1 <= A[i] <= 30000`
#
#
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,1,2,4]
#     Output: 17
#     Explanation:
#     Subarrays are [3], [1], [2], [4], [3,1], [1,2], [2,4], [3,1,2], [1,2,4], [3,1,2,4].
#     Minimums are 3, 1, 2, 4, 1, 1, 2, 1, 1, 1.
#     Sum is 17.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [11,81,94,43,3]
#     Output: 444
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 3 * 104`
#   * `1 <= arr[i] <= 3 * 104`
#
#
## @lc code=start
using LeetCode

function sum_subarray_mins(arr::Vector{Int})
    csum, res, md = 0, 0, Int(1e9) + 7
    stk = Int[]
    for i in eachindex(arr)
        while !isempty(stk) && arr[stk[end]] >= arr[i]
            tp = pop!(stk)
            ntop = isempty(stk) ? 0 : stk[end]
            csum += (arr[i] - arr[tp]) * (tp - ntop);
        end
        csum += arr[i]
        push!(stk, i)
        res += csum
    end
    return res % md
end
## @lc code=end


text:
function ODEState(arg0::jdouble, arg1::Vector{jdouble})
    return ODEState((jdouble, Vector{jdouble}), arg0, arg1)
end

function ODEState(arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return ODEState((jdouble, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function get_complete_state(obj::ODEState)
    return jcall(obj, "getCompleteState", Vector{jdouble}, ())
end

function get_complete_state_dimension(obj::ODEState)
    return jcall(obj, "getCompleteStateDimension", jint, ())
end

function get_number_of_secondary_states(obj::ODEState)
    return jcall(obj, "getNumberOfSecondaryStates", jint, ())
end

function get_primary_state(obj::ODEState)
    return jcall(obj, "getPrimaryState", Vector{jdouble}, ())
end

function get_primary_state_dimension(obj::ODEState)
    return jcall(obj, "getPrimaryStateDimension", jint, ())
end

function get_secondary_state(obj::ODEState, arg0::jint)
    return jcall(obj, "getSecondaryState", Vector{jdouble}, (jint,), arg0)
end

function get_secondary_state_dimension(obj::ODEState, arg0::jint)
    return jcall(obj, "getSecondaryStateDimension", jint, (jint,), arg0)
end

function get_time(obj::ODEState)
    return jcall(obj, "getTime", jdouble, ())
end


text:
#####
##### Timescale for diffusion across one cell
#####

using Oceananigans.Grids: topology, min_Δx, min_Δy, min_Δz

function min_Δxyz(grid, ::ThreeDimensionalFormulation)
    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    Δz = min_Δz(grid)
    return min(Δx, Δy, Δz)
end

function min_Δxyz(grid, ::HorizontalFormulation)
    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    return min(Δx, Δy)
end

min_Δxyz(grid, ::VerticalFormulation) = min_Δz(grid)

cell_diffusion_timescale(model) = cell_diffusion_timescale(model.closure, model.diffusivity_fields, model.grid)
cell_diffusion_timescale(::Nothing, diffusivities, grid) = Inf

maximum_numeric_diffusivity(κ::Number) = κ
maximum_numeric_diffusivity(κ::NamedTuple) = maximum(κ)
maximum_numeric_diffusivity(κ::NamedTuple{()}) = 0 # tracers=nothing means empty diffusivity tuples
maximum_numeric_diffusivity(::Nothing) = 0

# As the name suggests, we give up in the case of a function diffusivity
maximum_numeric_diffusivity(κ::Function) = 0

function cell_diffusion_timescale(closure::ScalarDiffusivity{TD, Dir}, diffusivities, grid) where {TD, Dir}
    Δ = min_Δxyz(grid, formulation(closure))
    max_κ = maximum_numeric_diffusivity(closure.κ)
    max_ν = maximum_numeric_diffusivity(closure.ν)
    return min(Δ^2 / max_ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::ScalarBiharmonicDiffusivity{Dir}, diffusivities, grid) where {Dir}
    Δ = min_Δxyz(grid, formulation(closure))
    max_κ = maximum_numeric_diffusivity(closure.κ)
    max_ν = maximum_numeric_diffusivity(closure.ν)
    return min(Δ^4/ max_ν, Δ^4 / max_κ)
end

function cell_diffusion_timescale(closure::SmagorinskyLilly, diffusivities, grid)
    Δ = min_Δxyz(grid, formulation(closure))
    min_Pr = closure.Pr isa NamedTuple{()} ? 1 : minimum(closure.Pr) # Innocuous value is there's no tracers
    max_νκ = maximum(diffusivities.νₑ.data.parent) * max(1, 1/min_Pr)
    return Δ^2 / max_νκ
end

function cell_diffusion_timescale(closure::AnisotropicMinimumDissipation, diffusivities, grid)
    Δ = min_Δxyz(grid, formulation(closure))
    max_ν = maximum(diffusivities.νₑ.data.parent)
    max_κ = diffusivities.κₑ isa NamedTuple{()} ? Inf : max(Tuple(maximum(κₑ.data.parent) for κₑ in diffusivities.κₑ)...)
    return min(Δ^2 / max_ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::TwoDimensionalLeith, diffusivities, grid)
    Δ = min_Δxyz(grid, ThreeDimensionalFormulation())
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

# Vertically-implicit treatment of vertical diffusivity has no time-step restriction
cell_diffusion_timescale(::ConvectiveAdjustmentVerticalDiffusivity{<:VerticallyImplicitTimeDiscretization},
                         diffusivities, grid) = Inf

cell_diffusion_timescale(closure::Tuple, diffusivities, grid) =
    min(Tuple(cell_diffusion_timescale(c, diffusivities, grid) for c in closure)...)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ aa854f02-46bf-4a2b-b1e9-5fb52b99925c
begin
    import Pkg
    Pkg.activate(Base.current_project())
    Pkg.instantiate()
end

# ╔═╡ 8b4426b8-588b-49ad-83c5-a79ed698d704
using PlutoUI

# ╔═╡ 3c07c33b-7f91-4e2d-a27a-09a8924f328c
using DataFrames

# ╔═╡ b59dd2c0-4c3e-4c3d-a4d4-f8da9dcbe286
using CSV

# ╔═╡ 76e2e659-c708-4a08-8263-def6c3dfa930
using Dates

# ╔═╡ 533d8e9e-40ff-4740-bbd9-8199cf0bcef6
using TemplateMatching

# ╔═╡ f07ab63f-dfb4-493c-a40c-6b3363a858c3
using Plots

# ╔═╡ 356cc4c2-e565-4847-9286-cfc1f835654d
using StatsBase

# ╔═╡ a7fcfdec-db16-4a40-b2f6-5c9c6f856ffc
using Printf

# ╔═╡ 33ff9417-b0b4-4f60-8c81-22c1181b37aa
using MultivariateStats

# ╔═╡ 3ec5251a-2e3e-4344-bb1f-fcd79daf9714
using LinearAlgebra

# ╔═╡ 9406c22e-fb52-461c-9b48-22bceb1b056c
md"## Reading templates catalogue"

# ╔═╡ 2fc4fdfe-81b2-45f1-9ed4-b8ad5bb86ee3
md"""Path of the CSV templates catalogue $(@bind templatespath TextField(default="../data/2021-01-12_20-25-30/cat.csv"))"""

# ╔═╡ 60ad8eb3-a313-4dce-82bd-a74a4d6842fd
begin
	import Base.+, Base.show

	+(x::Dates.UTInstant{Microsecond}, d::Microsecond) = Dates.UTInstant(x.periods + d)

	function show(io::IO, ::MIME"text/plain", x::Dates.UTInstant{Microsecond})
		date = DateTime(Dates.UTM(round(Int, 1e-3 * x.periods.value)))
		show(io, MIME"text/plain"(), date)
	end
end

# ╔═╡ 74063783-3454-4e4b-95d5-0427edbe9380
function DateTimeMicrosecond(y, m, d, h, mi, s, us)
	rata = us + 1_000_000 * (s + 60mi + 3600h + 86400 * Dates.totaldays(y, m, d))
	Dates.UTInstant{Microsecond}(Microsecond(rata))
end

# ╔═╡ 171de4ad-485a-4f92-818a-50cfbfa2330b
starttime = DateTimeMicrosecond(2021, 01, 12, 20, 25, 30, 0)

# ╔═╡ c76f6541-b0e2-4aba-9533-8ab6c23685bb
samplefreq = 2 # MHz

# ╔═╡ 56c57ae9-5254-4dca-a793-bfdeee9a0ba5
templates = let
	columns = [:Year, :Month, :Day, :Hour, :Minute, :Second, :North, :East, :Up]
	df = CSV.read(templatespath, DataFrame, select=columns)
	sec = floor.(Int, df.Second)
	usec = round.(Int, 1e6 * (df.Second - sec))
	templates = DataFrame()
	templates.datetime = DateTimeMicrosecond.(df.Year, df.Month, df.Day, df.Hour, df.Minute, sec, usec)
	templates.origin_sample = map(x -> round(Int, x.value * samplefreq), templates.datetime .- starttime)
	templates.origin_time = templates.origin_sample ./ samplefreq
	templates.north = df.North
	templates.east = df.East
	templates.up = df.Up
	templates
end

# ╔═╡ 57b35604-727b-4eac-bcbd-ea302e4c79a2
md"## Reading template matched catalogue"

# ╔═╡ b62d7318-3a6d-4442-b6cc-3166b87dff8a
md"""Path of the CSV template-match catalogue $(@bind cataloguepath TextField(default="../data/2021-01-12_20-25-30/augmented_catalogue.csv"))"""

# ╔═╡ 6c04792e-7438-4e36-a3cf-4d81da6a772c
catalogue = CSV.read(cataloguepath, DataFrame)

# ╔═╡ 102e4534-6be4-4436-be36-69726a393253
md" ## Statistics of the detections "

# ╔═╡ 2ce75a18-b5d6-4e6f-81ce-875527b14092
selection = let
	df = copy(catalogue)
	# Skip matches without enough channels
	dropmissing!(df)
	# Skip matches with low cc
	filter!(row -> row.crosscorrelation > 0.7 && row.nch > 6, df)
	# Select one template per match
	sort!(df, :origin_time)
	tokeep = TemplateMatching.selectbypeaksdistance(df.origin_time, df.crosscorrelation .* df.nch, 600)
	df = df[tokeep, :]
	# Drop funny relocalization residuals
	df[.!TemplateMatching.mad_test(df.multilateration_residual), :]
end

# ╔═╡ 08f3bb2c-c80e-43ff-814b-e17ae5a912dc
scatter(selection[:, :origin_time],
		templates[selection.template, :origin_time],
		zcolor=selection.crosscorrelation,
		c=:heat,
		markersize=4exp.(selection.relative_magnitude),
		ylabel="Template origin [μs]",
		xlabel="Detection origin [μs]",
		colorbar=:right,
		colorbar_title="\nCross-correlation",
   		right_margin = 3Plots.mm,
		legend=nothing,
		dpi=500)

# ╔═╡ 567cbd00-b375-417f-a7f5-8ac9a4bc4c7f
nontemplates = @. !(selection.crosscorrelation > 0.999 && selection.relative_magnitude ≈ 0.0)

# ╔═╡ feebcd5f-91c4-4a70-a064-e642e3172614
histogram(selection[nontemplates, :crosscorrelation], label=nothing, title="Cross-correlation")

# ╔═╡ 476101af-0108-4628-9ac1-f233a456a869
summarystats(selection[nontemplates, :crosscorrelation])

# ╔═╡ 3be237b0-9f19-4e9a-aa13-c6355771b8e8
histogram(selection[nontemplates, :relative_magnitude], label=nothing, title="Relative magnitude")

# ╔═╡ 7453767a-e8f5-4446-b198-fbe48c992aa5
summarystats(selection.relative_magnitude)

# ╔═╡ df8156eb-0061-4523-ba02-d7bb7eb4d9cb
histogram(selection[nontemplates, :nch], label=nothing, title="Num channels")

# ╔═╡ 97925230-7397-4cec-bce7-f15b3479b791
summarystats(selection[nontemplates, :nch])

# ╔═╡ 1f4c95b0-851b-4a91-bba6-3240cf312bf2
histogram(selection[nontemplates, :multilateration_residual], label=nothing, title="Multilateration residual")

# ╔═╡ 28efc2c9-ce88-4041-a84f-eda2845085e8
summarystats(selection[nontemplates, :multilateration_residual])

# ╔═╡ 2f98cf67-87ee-4f12-a507-b2259d780ffd
md"## Animations"

# ╔═╡ 6dddcebd-d730-470e-a831-154ee3102e0d
let nframes = 500, size_r = 0.8, alpha_r = 1e-6, data_len = Int(4e6)
	@gif for n = range(1, data_len, nframes)
		selection_upto = selection[selection.origin_time .<= n / samplefreq, :]
		if !isempty(selection_upto)
	    	scatter(selection_upto[!, :east],
					selection_upto[!, :north],
					selection_upto[!, :up],
					markersize=2.5,
					xlim=(0, 25),
					ylim=(0, 25),
					zlim=(0, 25),
					title=@sprintf("%.2f ms", 1e-3n / samplefreq),
					legend=nothing)
		end
	end
end

# ╔═╡ 1d576253-f719-4f21-8cc2-e794399a6ad5
let nframes = 500, size_r = 0.8, alpha_r = 1e-6, data_len = Int(4e6)
	@gif for n = range(1, data_len, nframes)
		selection_upto = selection[selection.peak_sample .<= n, :]
		if !isempty(selection_upto)
			count_upto = combine(groupby(selection_upto, :template), nrow => :count)
			lastsample_upto = combine(groupby(selection_upto, :template), :peak_sample => maximum)
			alphas = @. exp10(alpha_r * (lastsample_upto[:, :peak_sample_maximum] - n))
	    	scatter(templates[count_upto[!, :template], :east],
					templates[count_upto[!, :template], :north],
					templates[count_upto[!, :template], :up],
					markersize=size_r .* count_upto[:, :count],
					markeralpha=alphas,
					xlim=(0, 25),
					ylim=(0, 25),
					zlim=(0, 25),
					title=@sprintf("%.2f ms", 1e-3n / samplefreq),
					legend=nothing)
		end
	end
end

# ╔═╡ 0d237fe0-35bb-4624-a49a-6086c5d9ed09
md"## Principal Component Analysis"

# ╔═╡ 0131b25b-8b5c-4d1e-a520-830f5f20f345
X = transpose(Matrix(selection[!, [:north, :east, :up]]))

# ╔═╡ 09731759-c6c6-43ba-bc99-493327efa342
M = fit(PCA, X, maxoutdim=2)

# ╔═╡ 8c740959-ee9b-4139-a31b-fdd7ccc7f563
xlim, ylim = extrema(predict(M, X), dims=2)

# ╔═╡ 2837928e-692a-42ee-a064-07a3042d3d7b
let nframes = 500, data_len = Int(4e6)
	@gif for n = range(1, data_len, nframes)
		selection_upto = selection[selection.origin_time .<= n / samplefreq, :]
		if !isempty(selection)
			X = transpose(Matrix(selection_upto[!, [:north, :east, :up]]))
			Y = predict(M, X)
	    	scatter(Y[1, :],
					Y[2, :],
					xlim=xlim,
					ylim=ylim,
					title=@sprintf("%.2f ms", 1e-3n / samplefreq),
					legend=nothing)
		end
	end
end

# ╔═╡ 3212d37b-e498-4711-a067-14dca93a91a8
let
	Y = predict(M, X)
	plot(fit(Histogram, (Y[1, :], Y[2, :]), nbins=64))
end

# ╔═╡ 04d6e989-3fd0-4efd-a570-0ae34164cf4f
let
	Y = predict(M, X)
	Z = norm.(eachcol(X .- reconstruct(M, Y)))
	histogram(Z)
end

# ╔═╡ 432b6d5c-edef-402d-b207-a27b2dba18d1
md"## Candidate Repeaters"

# ╔═╡ c893eaba-5393-47ec-aff3-57352e515de4
repeaters = sort(filter(row -> row.crosscorrelation > 0.95, selection[nontemplates, :]), :template)

# ╔═╡ 9817148b-401b-45a8-ad0b-657502d6cb2f
scatter(repeaters.north, repeaters.east, repeaters.up)

# ╔═╡ 947f37b9-f8f9-482b-aca6-b28c90e1f5da
let
	templates = .!nontemplates
	df = vcat(repeaters, selection[templates, :])
	clusters = groupby(df, :template)
	errors = combine(clusters, ([:north, :east, :up] .=> std))
	map(col -> mean(filter(!isnan, errors[!, col])), [:north_std, :east_std, :up_std])
end

# ╔═╡ Cell order:
# ╠═aa854f02-46bf-4a2b-b1e9-5fb52b99925c
# ╠═8b4426b8-588b-49ad-83c5-a79ed698d704
# ╟─9406c22e-fb52-461c-9b48-22bceb1b056c
# ╠═3c07c33b-7f91-4e2d-a27a-09a8924f328c
# ╠═b59dd2c0-4c3e-4c3d-a4d4-f8da9dcbe286
# ╠═76e2e659-c708-4a08-8263-def6c3dfa930
# ╠═2fc4fdfe-81b2-45f1-9ed4-b8ad5bb86ee3
# ╠═74063783-3454-4e4b-95d5-0427edbe9380
# ╠═171de4ad-485a-4f92-818a-50cfbfa2330b
# ╠═60ad8eb3-a313-4dce-82bd-a74a4d6842fd
# ╠═c76f6541-b0e2-4aba-9533-8ab6c23685bb
# ╠═56c57ae9-5254-4dca-a793-bfdeee9a0ba5
# ╟─57b35604-727b-4eac-bcbd-ea302e4c79a2
# ╠═b62d7318-3a6d-4442-b6cc-3166b87dff8a
# ╠═533d8e9e-40ff-4740-bbd9-8199cf0bcef6
# ╠═f07ab63f-dfb4-493c-a40c-6b3363a858c3
# ╠═6c04792e-7438-4e36-a3cf-4d81da6a772c
# ╟─102e4534-6be4-4436-be36-69726a393253
# ╠═356cc4c2-e565-4847-9286-cfc1f835654d
# ╠═2ce75a18-b5d6-4e6f-81ce-875527b14092
# ╠═08f3bb2c-c80e-43ff-814b-e17ae5a912dc
# ╠═567cbd00-b375-417f-a7f5-8ac9a4bc4c7f
# ╟─feebcd5f-91c4-4a70-a064-e642e3172614
# ╟─476101af-0108-4628-9ac1-f233a456a869
# ╟─3be237b0-9f19-4e9a-aa13-c6355771b8e8
# ╟─7453767a-e8f5-4446-b198-fbe48c992aa5
# ╟─df8156eb-0061-4523-ba02-d7bb7eb4d9cb
# ╟─97925230-7397-4cec-bce7-f15b3479b791
# ╟─1f4c95b0-851b-4a91-bba6-3240cf312bf2
# ╟─28efc2c9-ce88-4041-a84f-eda2845085e8
# ╟─2f98cf67-87ee-4f12-a507-b2259d780ffd
# ╠═a7fcfdec-db16-4a40-b2f6-5c9c6f856ffc
# ╠═6dddcebd-d730-470e-a831-154ee3102e0d
# ╠═1d576253-f719-4f21-8cc2-e794399a6ad5
# ╟─0d237fe0-35bb-4624-a49a-6086c5d9ed09
# ╠═33ff9417-b0b4-4f60-8c81-22c1181b37aa
# ╠═0131b25b-8b5c-4d1e-a520-830f5f20f345
# ╠═09731759-c6c6-43ba-bc99-493327efa342
# ╠═8c740959-ee9b-4139-a31b-fdd7ccc7f563
# ╠═2837928e-692a-42ee-a064-07a3042d3d7b
# ╠═3ec5251a-2e3e-4344-bb1f-fcd79daf9714
# ╠═3212d37b-e498-4711-a067-14dca93a91a8
# ╠═04d6e989-3fd0-4efd-a570-0ae34164cf4f
# ╠═432b6d5c-edef-402d-b207-a27b2dba18d1
# ╠═c893eaba-5393-47ec-aff3-57352e515de4
# ╠═9817148b-401b-45a8-ad0b-657502d6cb2f
# ╠═947f37b9-f8f9-482b-aca6-b28c90e1f5da


text:
module Elements

export Element, Singleton, Group, kind, @kind, circulation, flux, streamfunction,
       complexpotential, conftransform, inverse_conftransform, jacobian, image,
       blobradius

using ..Properties

import ..Utils: ComplexDual, Dual, ComplexGradientConfig, seed!, seed,
                vector_mode_gradient, vector_mode_jacobian, checktag

abstract type Element end

#== Trait definitions ==#

abstract type Singleton end
abstract type Group end

function kind end

macro kind(element, k)
    esc(quote
        Elements.kind(::$element) = $k
        Elements.kind(::Type{$element}) = $k
        end)
end

@kind ComplexF64 Singleton

kind(::ComplexDual) = Singleton
kind(::Dual) = Singleton
kind(::AbstractArray{T}) where {T <: Union{Element, ComplexF64, Dual, ComplexDual}} = Group
kind(::Tuple) = Group

# Convenience functions to define wrapper types
# e.g. vortex sheets as a wrapper around vortex blobs
unwrap_src(e) = unwrap(e)
unwrap_targ(e) = unwrap(e)
unwrap(e) = e

## Actual Definitions of Properties

"""
    Elements.position(src::Element)

Returns the complex position of a potential flow element.
This is a required method for all `Element` types.

# Example

```jldoctest
julia> point = Vortex.Point(1.0 + 0.0im, 1.0);

julia> Elements.position(point)
1.0 + 0.0im

julia> points = Vortex.Point.([1.0im, 2.0im], 1.0);

julia> Elements.position.(points)
2-element Array{Complex{Float64},1}:
 0.0 + 1.0im
 0.0 + 2.0im
```
"""
@property begin
    signature = position(src::Source)
    stype = ComplexF64
end

"""
    Elements.blobradius(src)

Returns the blob radius of a regularized potential flow element.
If the element is singular, returns zero.
"""
@property begin
    signature = blobradius(src::Source)
    stype = Float64
end

"""
    Elements.circulation(src)

Returns the total circulation contained in `src`.

# Example

```jldoctest
julia> points = Vortex.Point.([1.0im, 2.0im], [1.0, 2.0]);

julia> blobs = Vortex.Blob.([1.0im, 2.0im], [1.0, 2.0], 0.1);

julia> Elements.circulation(points[1])
1.0

julia> Elements.circulation(points)
3.0

julia> Elements.circulation((points, blobs))
6.0

julia> Elements.circulation.(points)
2-element Array{Float64,1}:
 1.0
 2.0

julia> Elements.circulation.((points, blobs))
(3.0, 3.0)

julia> Elements.circulation(Source.Point(rand(), rand()))
0.0

julia> Elements.circulation(Source.Blob(rand(), rand(), rand()))
0.0
```
"""
@property begin
    signature = circulation(src::Source)
    reduce = (+)
    stype = Float64
end

"""
    Elements.flux(src)

Returns the flux through a unit circle induced by `src`.

# Example

```jldoctest
julia> points = Source.Point.([1.0im, 2.0im], [1.0, 2.0]);

julia> blobs = Source.Blob.([1.0im, 2.0im], [1.0, 2.0], 0.1);

julia> Elements.flux(points[1])
1.0

julia> Elements.flux((points, blobs))
6.0

julia> Elements.flux.(points)
2-element Array{Float64,1}:
 1.0
 2.0

julia> Elements.flux.((points, blobs))
(3.0, 3.0)

julia> Elements.flux(Vortex.Point(rand(), rand()))
0.0

julia> Elements.flux(Vortex.Blob(rand(), rand(), rand()))
0.0
```
"""
@property begin
    signature = flux(src::Source)
    reduce = (+)
    stype = Float64
end

@doc raw"""
    Elements.impulse(src)

Return the aerodynamic impulse of `src` about (0,0):
```math
P := \int \boldsymbol{x} \times \boldsymbol{\omega}\,\mathrm{d}A.
```
This is a required method for all vortex types.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(1.0im, -π, 0.1));

julia> Elements.impulse(sys[1])
3.141592653589793 + 0.0im

julia> Elements.impulse(sys)
0.0 + 0.0im
```
"""
@property begin
    signature = impulse(src::Source)
    reduce = (+)
    stype = ComplexF64
end

@doc raw"""
    Elements.angularimpulse(src)

Return the aerodynamic angular impulse of `src` about (0,0):
```math
\Pi := \frac{1}{2}\int \boldsymbol{x} \times (\boldsymbol{x}\times\boldsymbol{\omega})\,\mathrm{d}A.
```
This is a required method for all vortex types.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(1.0im, -π, 0.1));

julia> Elements.angularimpulse(sys[1])
3.141592653589793 + 0.0im

julia> Elements.angularimpulse(sys)
0.0 + 0.0im
```
"""
@property begin
    signature = angularimpulse(src::Source)
    reduce = (+)
    stype = ComplexF64
end

@property begin
    signature = streamfunction(targ::Target, src::Source)
    preallocator = allocate_streamfunction
    stype = Float64
end

@property begin
    signature = complexpotential(targ::Target, src::Source)
    preallocator = allocate_complexpotential
    stype = ComplexF64
end

@doc raw"""
    Elements.conftransform(src,body)

Return the conformally transformed position of `src` via the transform defined
by `body`. In this context, the position(s) in `src` is interpreted as
coordinates in the circle plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = ConformalBody([1/4,0,1/4]);

julia> Elements.conftransform(sys,b)
(0.0 + 0.0im, 0.0 + 0.375im)
```
"""
@property begin
    signature = conftransform(src::Source,b)
    preallocator = allocate_conftransform
    stype = ComplexF64
end

@doc raw"""
    Elements.inverse_conftransform(src,body)

Return the inverse conformally transformed position of `src` via the transform defined
by `body`. In this context, the position(s) in `src` is interpreted as
coordinates in the physical plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = ConformalBody([1/4,0,1/4]);

julia> Elements.inverse_conftransform(sys,b)
(0.0 + 0.0im, 0.0 + 0.375im)
```
"""
@property begin
    signature = inverse_conftransform(src::Source,b)
    preallocator = allocate_inv_conftransform
    stype = ComplexF64
end

@doc raw"""
    Elements.jacobian(src,body)

Return the Jacobian of the conformal transform at the position of `src` via the transform defined
by `body`. In this context, the position(s) in `src` is interpreted as
coordinates in the circle plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = ConformalBody([1/4,0,1/4]);

julia> Elements.jacobian(sys,b)
(0.5 + 0.0im, 0.0 + 0.5im)

```
"""
@property begin
    signature = jacobian(src::Source,b)
    preallocator = allocate_jacobian
    stype = ComplexF64
end

@doc raw"""
    Elements.image(src,body)

Return the image position of `src` in the circle-plane representation of `body`.
In this context, the position(s) in `src` is interpreted as
coordinates in the circle plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = PowerBody([1/4,0,1/4],zero(ComplexF64),0.0)

julia> Elements.image(sys,b)
(0.0 + 1.0im, 0.0 + 0.5im)

```
"""
@property begin
    signature = image(src::Source,b)
    preallocator = allocate_image
    stype = ComplexF64
end

"""
    seed_strength(T,v::Vector{Element},i::Int)

Given a collection `v` of points or blobs, create a copy of the collection with
the `i`th element's strength replaced by a unit dual. The entire output
collection has positions of `Dual` type, but only the replaced element has unit
partials; the others have partials equal to zero. `T` is a `Tag`, and can be set
to `Nothing`.
"""
function seed_strength end

"""
    property_type(::Element)

Return the type of the properties of the element.
"""
function property_type end

ComplexGradientConfig(f::F,v::Vector{<:Element},w...) where {F} =
      ComplexGradientConfig(f,position(v),w...)

"""
    seed_zeros(v::Vector{Element},cfg::ComplexGradientConfig)

Given a collection `v` of points or blobs, create a copy of the collection with
the element's positions and strengths replaced by complex duals with zeros.
`cfg` specifies the type and size of duals. The use of complex duals (in spite of the real-valued
strength) ensures correct treatment for gradient calculation, since some
calculations will be complex-valued.
"""
function seed_zeros(strength::F,v::Vector{<:Element},cfg::ComplexGradientConfig) where F
  posduals = convert.(eltype(cfg.duals),position(v))
  strduals = convert.(eltype(cfg.duals),complex(strength.(v)))
  return posduals, real(strduals)
end

"""
    seed_position(v::Vector{Element},cfg::ComplexGradientConfig)

Given a collection `v` of points or blobs, create a copy of the collection with
the element's positions replaced by complex duals with unit seeds and with
the element's strengths replaced by complex duals with zeros. `cfg` specifies
the type and size of duals. The length of the partials in these duals is twice
the number of elements, to enable computation of gradients with respect to
all positions.
"""
function seed_position(strength::F,v::Vector{<:Element},cfg::ComplexGradientConfig) where F
  posduals = copy(cfg.duals)
  strduals = copy(cfg.duals)
  seed!(posduals,position(v),cfg.rseeds,cfg.iseeds)
  seed!(strduals,complex(strength.(v)))
  return posduals, real(strduals)
end

"""
    seed_strength(v::Vector{Element},cfg::ComplexGradientConfig)

Given a collection `v` of points or blobs, create a copy of the collection with
the element's strengths replaced by complex duals with unit seeds and with
the element's positions replaced by complex duals with zeros. `cfg` specifies
the type and size of duals. The use of complex duals (in spite of the real-valued
strength) ensures correct treatment for gradient calculation, since some
calculations will be complex-valued.
"""
function seed_strength(strength::F,v::Vector{<:Element},cfg::ComplexGradientConfig) where F
  posduals = copy(cfg.duals)
  strduals = copy(cfg.duals)
  seed!(posduals,position(v))
  seed!(strduals,complex(strength.(v)),cfg.rseeds,cfg.iseeds)
  return posduals, real(strduals)
end

function gradient_position(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    vector_mode_gradient(f,v,cfg,vector_mode_dual_eval_position)
end

function gradient_strength(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    dz, dzstar = vector_mode_gradient(f,v,cfg,vector_mode_dual_eval_strength)
    return dz .+ dzstar
end

@inline vector_mode_dual_eval_position(f::F, v::Vector{<:Element},
          cfg::ComplexGradientConfig) where {F} = f(seed_position(v,cfg))

@inline vector_mode_dual_eval_strength(f::F, v::Vector{<:Element},
          cfg::ComplexGradientConfig) where {F} = f(seed_strength(v,cfg))

@inline vector_mode_dual_eval_param(f::F, v::Tuple{Vector{<:Element},T},
          cfg::ComplexGradientConfig) where {F,T} = f(seed_zeros(v[1],cfg),seed(v[2],cfg))

function jacobian_position(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    vector_mode_jacobian(f,v,cfg,vector_mode_dual_eval_position)
end

function jacobian_strength(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    dz, dzstar = vector_mode_jacobian(f,v,cfg,vector_mode_dual_eval_strength)
    return dz .+ dzstar
end

function jacobian_param(f,v::Tuple{Vector{<:Element},S},
        cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, ComplexF64[v[2]])) where {S,T}
    checktag(T, f, v[2])
    dz, dzstar = vector_mode_jacobian(f,v,cfg,vector_mode_dual_eval_param)
    return real(dz .+ dzstar)
end

end


text:
# Alignment Anchor
# ================
#
# Sequence alignment anchor type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Alignment operation with anchoring positions.
"""
struct AlignmentAnchor
    seqpos::Int
    refpos::Int
    alnpos::Int
    op::Operation
end

seqpos(anc::AlignmentAnchor) = anc.seqpos
refpos(anc::AlignmentAnchor) = anc.refpos
alnpos(anc::AlignmentAnchor) = anc.alnpos

function Base.show(io::IO, anc::AlignmentAnchor)
    print(io, "AlignmentAnchor(", anc.seqpos, ", ", anc.refpos,
                               ", ", anc.alnpos, ", '", anc.op, "')")
end


text:
module ViewHelper

using Genie, Genie.Flash, Genie.Router

export output_flash

function output_flash(flashtype::String = "danger") :: String
  flash_has_message() ? """
$(flash())
""" : ""
end

end


text:
py"""
import numpy
import torch

def tensor2numpy(x):
    return x.detach().to("cpu").numpy()

def numpy2tensor(x):
    return torch.from_numpy(x)

def cal_param(model):
    return sum(p.numel() for p in model.parameters() if p.requires_grad)

def load_ckpt(ckpt_path):
    return torch.load(ckpt_path, map_location= 'cpu')

def test_model(model, img):
    with torch.no_grad():
        model.eval()
        return model.forward(img)
"""


text:
struct Experience{T<:Real, V<:AbstractArray{T}, A, F}
    s::V
    a::A
    r::F
    s′::V
    done::Bool
end

abstract type ReplayMemoryBuffer end

mutable struct PriorityReplayMemoryBuffer{T, S<:Sequence} <: ReplayMemoryBuffer
    capacity::Int
    experience::CircularBuffer{Experience}
    priorities::CircularBuffer{T}
    currentMax::T
    batch_size::Int
    α::Float32
    β::S
    ϵ::Float32
end

# Constructor for Empty Memory
function PriorityReplayMemoryBuffer(n::Int, batch_size::Int; β::D, α=5f-1, ϵ=1f-3) where {D<:Sequence}
    PriorityReplayMemoryBuffer(n, CircularBuffer{Experience}(n), CircularBuffer{Float32}(n), ϵ, batch_size, α, β, ϵ)
end

# Utility functions
Base.length(mem::PriorityReplayMemoryBuffer) = length(mem.experience)
Base.size(mem::PriorityReplayMemoryBuffer) = length(mem)

# Memory Control
function addexp!(mem::PriorityReplayMemoryBuffer, Exp::Experience)
    push!(mem.experience, Exp)
    push!(mem.priorities, mem.currentMax)
end

function addexp!(mem::PriorityReplayMemoryBuffer, s::AbstractArray{T}, a::A,
                 r::F, s′::AbstractArray{T}, d::Bool) where {T, A, F}
    addexp!(mem, Experience(s, a, convert(Float32, r), s′, d))
end

# Parameter Control
update_β!(mem::PriorityReplayMemoryBuffer) = step!(mem.β)

# Update priorities for selected indicies - updated while training
function update_priorities!(mem::PriorityReplayMemoryBuffer, ids::Vector{T}, td_errs::V ) where {T<:Int, V <: AbstractArray}
    mem.priorities[ids] = (abs.(td_errs) .+ mem.ϵ).^mem.α
    mem.currentMax = maximum(mem.priorities)
end

# Sample from the buffer
function StatsBase.sample(mem::PriorityReplayMemoryBuffer)
    ids = sample(1:length(mem), Weights(mem.priorities) , mem.batch_size, replace=false)
    s = hcat((mem.experience[i].s for i in ids)...)
    r = hcat((mem.experience[i].r for i in ids)...)
    s′ = hcat((mem.experience[i].s′ for i in ids)...)
    d = hcat((mem.experience[i].done for i in ids)...)
    weights = mem.priorities[ids] ./ sum(mem.priorities)
    weights = (length(mem)*weights).^(-value(mem.β))
    weights = weights ./ maximum(weights)

    # Actions need to be converted to Cartesian indices so that they address
    # into the correct place
    a = [CartesianIndex(0,0) for i in ids]
    for (i, idx) in enumerate(ids)
        @inbounds a[i] = CartesianIndex(mem.experience[idx].a, i)
    end
    return (s, a, r, s′, d, ids, weights)
end


text:
#= to check : Raviart-Thomas assembling =#

function getBaseRT0(E::Float64, xE::Array{Float64, 2})
    return [
        x -> (x - xE[:, 3])/(2*E),
        x -> (x - xE[:, 1])/(2*E),
        x -> (x - xE[:, 2])/(2*E)
    ]
end

function melemRT0(vtx::Array{Float64, 2})
    @assert size(vtx, 2) == 3
    # initialize
    E = polygonArea(vtx)
    # div part
    DIV = ones(Float64, 3, 3) / E
    # mass part
    RT0 = getBaseRT0(E, vtx)
    MASS = zeros(Float64, 3, 3)
    quad = quadGL2d(3, vtx)
    for j=1:3
        for i=1:3
            MASS[i,j] = integrate(quad, x -> dot(RT0[i](x),RT0[j](x)))
        end
    end
    return DIV, MASS
end

function buildRT0matrices(msh::mesh, conn::Connectivity)
    e2dof = conn.e2dof
    ndof  = conn.ndof
    # 1st: preallocate
    A = zeros(Float64, ndof, ndof)
    B = zeros(Float64, ndof, ndof)
    # 2nd: main loop
    for e=1:length(msh.elt)
        @inbounds verts  = msh.vtx[:, msh.elt[e]]
        AE, BE = melemRT0(verts)
        nloc = length(e2dof[e])
        for i = 1:nloc
            @inbounds absI = abs(e2dof[e][i])
            @inbounds sgnI = sign(e2dof[e][i])
            for j=1:nloc
                @inbounds absJ = abs(e2dof[e][j])
                @inbounds sgnJ = sign(e2dof[e][j])
                @inbounds A[absI, absJ] += sgnI*sgnJ*(AE[i,j] + BE[i,j])
                @inbounds B[absI, absJ] += sgnI*sgnJ*BE[i,j]
                #
            end
        end
    end
    return A, B
end


text:
@testset "FlowrateProblem" begin

    @testset "constant solution" begin
    prob = FlowrateProblem(DiffusionEquation{2}(identity),
                           i=0.1, Qb=0)

    θ = solve(prob)

    ϕ = range(1e-6, 20, length=100)

    @test all(@inferred θ.(ϕ) .≈ θ.i)
    @test all(@inferred θ.(ϕ) .≈ θ.b)
    @test all(@inferred d_dϕ.(θ, ϕ) .≈ 0)
    @test θ.iterations == 0
    end

    @testset "Qb" begin
    Qb = 2
    height = 10
    θi = 0.1

    prob = FlowrateProblem(DiffusionEquation{2}(identity),
                           i=θi, Qb=Qb, height=height)

    θ = solve(prob)

    t = [1e-6, 1, 1.5, 5, 7.314]

    @test all(@inferred flux.(θ, :b, t) ≈ Qb./(2π.*rb.(θ, t).*height))
    @test isapprox(θ.i, θi, atol=1e-3)
    @test θ.iterations > 0
    @test all(isnan.(θ.(0,t)))
    end

end


text:
function _write_input{K<:Field, N}(si, F::Array{Polynomial{K,N},1}; perm=1:N)
    println(si, N, " ", length(F))
    for f in F
        if length(f)==0 # print as if it is a constant term = 0
            println(si, 1)
            for k=1:N
                print(si, "0 ")
            end
            println(si, 0)
        else
            println(si, length(f))
            h = lcm([den(c) for (i,c) in f])
            for (i,c) in f
                for k=1:N
                    print(si, i[perm[k]], " ")
                end
                println(si, num(c*h))
            end
        end
    end
end

function _read_output(so, K::DataType; perm=nothing)
    N = parse(Int, readuntil(so,' '))
    if perm==nothing
        perm=1:N
    end
    ip = invperm(perm)
    n_poly = parse(Int, readuntil(so,'\n')[1:end-1])
    F = Polynomial{K,N}[]
    for p=1:n_poly
        n_terms = parse(Int, readuntil(so,'\n')[1:end-1])
        f = Polynomial{K,N}([((UInt16[parse(UInt16, readuntil(so,' '))
                                      for i=1:N][ip]...),
                    convert(K,parse(BigInt, readuntil(so,'\n')[1:end-1])))
                    for t=1:n_terms])
        push!(F,f)
    end
    F
end

function groebner_fgb{K<:Field, N}(F::Polynomial{K,N} ...)
    call_fgb = joinpath(dirname(@__FILE__),
               "..", "deps", "bin", "call_fgb")
    (so,si,pr) = readandwrite(`$call_fgb`)

    _write_input(si, [F...])

    G = _read_output(so, K)

    close(so)
    close(si)
    close(pr)

    G
end

function call_fgb{K<:Field, N}(F::Array{Polynomial{K,N},1};
              perm=collect(1:N),
              verbose::Integer=0,
              n_threads::Integer=1,
              k_elim::Integer=0,
              compute::Integer=1,
              nb::Integer=0,
              n_output_max::Integer=100000,
              force_elim::Bool=false,
              index::Integer=1000000)
    @assert n_threads>0
    @assert k_elim>=0 && k_elim <=N
    @assert compute >=1 && compute <=12
    @assert !force_elim || k_elim>0
    @assert n_output_max>0
    @assert index>0
    @assert length(perm)==N && sort(perm)==collect(1:N)

    opt = [string("-v", verbose),
           string("-t", n_threads),
           string("-k", k_elim),
           string("-o", n_output_max),
           string("-c", compute),
           string("-n", nb),
           string("-m", index),
           string("-e", force_elim?1:0)];
    call_fgb = joinpath(dirname(@__FILE__),
               "..", "deps", "bin", "call_fgb")
    cmd = `$call_fgb $opt`

    (so,si,pr) = readandwrite(cmd)

    _write_input(si, F; perm=perm)

    G = _read_output(so, K; perm=perm)

    close(so)
    close(si)
    close(pr)

    G
end

function _adapt_vars_args(vars1::Vector, vars2::Vector)
    k=length(vars2)
    vars = vcat(vars1, vars2)
    N=length(vars)
    for v in vars
        @assert sum(v.index)==1
    end
    perm = Int[findfirst(v.index,1) for v in vars]
    @assert length(perm)==N && sort(perm)==collect(1:N)
    (k, perm)
end

function fgb_gbasis{K<:Field, N}(F::Array{Polynomial{K,N},1},
                                 vars1::Vector,
                                 vars2::Vector;
                                 verbose::Integer=0,
                                 n_threads::Integer=1,
                                 n_output_max::Integer=100000,
                                 index::Integer=1000000)
    (k, perm) = _adapt_vars_args(vars1, vars2)
    call_fgb(F; perm=perm, k_elim=k, verbose=verbose,
             n_threads=n_threads, n_output_max=n_output_max, index=index)
end

function fgb_gbasis_elim{K<:Field, N}(F::Array{Polynomial{K,N},1},
                                 vars1::Vector,
                                 vars2::Vector;
                                 verbose::Integer=0,
                                 n_threads::Integer=1,
                                 n_output_max::Integer=100000,
                                 index::Integer=1000000)
    (k, perm) = _adapt_vars_args(vars1, vars2)
    call_fgb(F; perm=perm, k_elim=k, force_elim=true, verbose=verbose,
             n_threads=n_threads, n_output_max=n_output_max, index=index)
end

function call_giac{K<:Field, N}(F::Array{Polynomial{K,N},1};
              perm=collect(1:N), greduce::Bool=false)
    @assert length(perm)==N && sort(perm)==collect(1:N)

    call_giac = joinpath(dirname(@__FILE__),
               "..", "deps", "bin", "call_giac")
    opt = [ string("-r", greduce?1:0)]
    cmd = `$call_giac $opt`

    (so,si,pr) = readandwrite(cmd)

    _write_input(si, F; perm=perm)

    G = _read_output(so, K; perm=perm)

    close(so)
    close(si)
    close(pr)

    G
end

function giac_gbasis{K<:Field, N}(F::Array{Polynomial{K,N},1}, vars::Vector)
    (k, perm) = _adapt_vars_args(vars, [])
    call_giac(F; perm=perm)
end

function giac_greduce{K<:Field, N}(f::Polynomial{K,N}, G::Array{Polynomial{K,N},1}, vars::Vector)
    (k, perm) = _adapt_vars_args(vars, [])
    r = call_giac(vcat(f, G); perm=perm, greduce=true)
    r[1]/first(values(r[2]))
end


text:
using LearnConvection
using Plots

D=32
N=4

train = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2"]
test = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2"]

# problem  = Slack("TKE"; parameters=TKEMassFlux.TKEParameters())
problem  = Sequential("dT")

train = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2", "general_strat_12_profiles.jld2", "general_strat_16_profiles.jld2"]
test = ["general_strat_28_profiles.jld2", "general_strat_32_profiles.jld2"]

k = 3
logγ = -1.7
distance = euclidean_distance
kernel   = get_kernel(k, logγ, 0.0, distance)

𝒟_train  = LearnConvection.Data.data(train, problem; D=D, N=N);
𝒟_test   = LearnConvection.Data.data(test, problem; D=D, N=N);

all_errors = zeros(D,D)

for stencil_size=1:D
    𝒢 = LearnConvection.GaussianProcess.model(𝒟_train; kernel = kernel, stencil_size = stencil_size)
    all_errors[stencil_size,:] = disparity_vector(𝒢, 𝒟_test)
end

p = plot(all_errors', yscale=:log10, legend=false, xlabel="Stencil size", ylabel="Error", title="Slack(TKE): Error for each gridpoint")

ns = 1:4:32
p = plot(size=(900,600))
for n in ns
    label = 𝒟_train.zavg[n]
    plot!(all_errors[:,n], linewidth=3, yscale=:log10, c=cgrad(:acton)[n*8], xlabel="Stencil size", ylabel="Error", label=label, title="Slack(TKE): Error at different depths")
end
p
png(p,"error_diff_depths2.png")

ns = 1:4:32
p = plot(size=(900,600))
for n in ns
    label = 𝒟_train.zavg[n]
    plot!(collect(17:32), all_errors[17:32,n], linewidth=3, c=cgrad(:acton)[n*8], xlabel="Stencil size", ylabel="Error", label=label, title="Slack(TKE): Error at different depths")
end
p
png(p,"error_diff_depths_linear.png")

mean_errors = [sum(all_errors[stencil_size,:])/D for stencil_size in 1:D]
s = plot(mean_errors, yscale=:log10, legend=false, xlabel="Stencil size", ylabel="Error", title="Slack(TKE): Mean error")
png(s,"mean_error.png")
s2 = plot(collect(16:32), mean_errors[16:32], legend=false, xlabel="Stencil size", ylabel="Error", title="Slack(TKE): Mean error")
png(s2,"mean_error_linear.png")

# png(anim, "error_vs_stencil_size.gif")

# anim = animate_profile(𝒢, 𝒟_test)
anim = animate_profile(𝒢, 𝒟_train)

anim = animate_profile_and_model_output(𝒢, 𝒟_train)
gif(anim, "all.gif")

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
@testset "Differential Private Kalman Filtering Tests" begin
  @testset "Static Input Block Design" begin
    ϵ = log(3); δ = 0.05
    k_priv = gaussianMechConstant(ϵ,δ)

    nusers = 10; ρ = ones(1,nusers)
    Ls = ones(1,1,nusers); As = 0.9*ones(1,1,nusers); Cs = ones(1,1,nusers)
    Vs = 0.01*ones(1,1,nusers); Ws = 0.2*ones(1,1,nusers)
    Vinvs = zero(Vs); Winvs = zero(Ws)
    for i=1:nusers
    	Vinvs[:,:,i] = inv(Vs[:,:,i])
    	Winvs[:,:,i] = inv(Ws[:,:,i])
    end
    (D, P_val, X_val, Ω_val, M) =
      staticInputBlock_DPKF_ss(Ls, As, Cs, Winvs, Vs, Vinvs, ρ, k_priv)
    cost = evaluateKFperf(D, Ls, As, Cs, Ws, Vs, ρ, k_priv)[1]

    @test sum(D) ≈ nusers atol=1e-6 # D should be all ones
    @test cost ≈ 1.6244783332772452
    @test tr(X_val) ≈ cost atol=1e-6

    # LQG test - need to finish, writing evaluation function for LQG
    Bs = ones(1,1,nusers)
    Q = 1.0 * Matrix(I,nusers,nusers); R = 0.5
    (Dlqg, perflqg, Plqg, X_val_lqg, Mlqg) =
      staticInputBlock_DPLQG_ss(Q, R, As, Bs, Cs, Ws, Winvs, Vs, Vinvs, ρ, k_priv)
    costlqg = evaluateLQGperf(Dlqg, Q, R, As, Bs, Cs, Ws, Vs, ρ, k_priv)[1]
    @test perflqg ≈ costlqg atol=1e-6
  end
end


text:
using Profile

input = "318946572"
# input = "389125467" # test pattern
orig_cups = parse.(Int, collect(input))

function play(cups, n)
    cup_idx = 1
    min_cup, max_cup = extrema(cups)
    picked = Array{Int64, 1}(undef, 3)

    for i in 1:n
        if (i % 100_000 == 0)
            println(i)
        end
        # println("-- move $(i) --")
        # println("cups: $(join(cups, " "))")
        cup = cups[1]
        # println("current cup: $(cup)")
        pick_idxs = 2:4
        picked = cups[pick_idxs]
        # remove picks from cups
#        filter!(c -> c ∉ picked, cups)
        # println("pick up: $(join(picked, " ")) ($(join(pick_idxs, " ")))")
        destination = cup - 1
        if destination < min_cup
            destination = max_cup
        end
        while destination in picked
            destination -= 1
            if destination < min_cup
                destination = max_cup
            end
        end
        dst_idx = findfirst(x -> x == destination, cups)
        # println("destination cup: $(destination) ($(dst_idx))")
        # insert cups
        cups[2:dst_idx-3] = cups[5:dst_idx]
        # println(cups)
        cups[(dst_idx-2):(dst_idx)] = picked
        # println(cups)
        cups[1:end-1] = cups[2:end]
        cups[end] = cup
    end

    return first_idx = findfirst(x -> x == 1, cups)
end

function play2(cups, n)
    list = Array{Int64, 1}(undef, length(cups))
    for i in 1:(length(cups)-1)
        list[cups[i]] = cups[i+1]
    end
    list[cups[end]] = cups[1]
    # println(list)
    cup = cups[1]
    min_cup, max_cup = extrema(cups)
    for i in 1:n
        # println("-- move $(i) --")
        # println("cup: ", cup)
        picked = [list[cup], list[list[cup]], list[list[list[cup]]]]
        # println("pick up: ", join(picked, ", "))
        destination = cup - 1
        if destination < min_cup
            destination = max_cup
        end
        while destination in picked
            destination -= 1
            if destination < min_cup
                destination = max_cup
            end
        end
        # println("destination: ", destination)
        list[cup] = list[picked[end]]
        list[picked[end]] = list[destination]
        list[destination] = picked[1]
        cup = list[cup]
    end

    idx = 1
    for i in 1:(length(cups))
        cups[i] = list[idx]
        idx = list[idx]
    end

    return cups
end

cups = copy(orig_cups)
i = play(cups, 100)
println(join(vcat(cups[(i+1):end], cups[1:(i-1)])))

cups = vcat(copy(orig_cups), collect((maximum(orig_cups)+1):1_000_000))
play2(cups, 10_000_000)
println(cups[1]*cups[2])


text:
export bm_solver

"""
`history_check(hist::Dict, g::Int)` sees if the guess `g` is consistent
with the previous answers stored in `hist`.
"""
function history_check(history::Dict{Int,Tuple{Int,Int}}, g::Int)::Bool
    for d in keys(history)
        if history[d] != bm_count(g, d)
            return false
        end
    end
    return true
end

"""
`history_count(h)` returns the number of possible codes that
are consistent with the history dictionary `h`.
"""
function history_count(history::Dict{Int,Tuple{Int,Int}})::Int
    length([c for c in 0:9999 if history_check(history, c)])
end

"""
`bm_solver(c::Int, verbose::Bool=true, delay::Int=0)` playes Bullets and Maggots to try
to find code `c`. Returns the number of guesses. It doesn't cheat!

* `verbose` controls printing during execution.
* `delay` indicates how many steps we should wait before using the history; `0` means no delay
"""
function bm_solver(c::Int, verbose::Bool = true, first_guess = guess_any)::Int
    if !code_check(c)
        throw(bad_code_message(c))
        return 0
    end
    if verbose
        println("Trying to find secret code $(string4(c))")
    end

    history = Dict{Int,Tuple{Int,Int}}()   # place to hold past guesses/results
    result_server(g::Int) = bm_count(g, c)
    bm_solver_engine(result_server, verbose, history, first_guess)
end

function bm_solver_engine(
    result_server::Function,
    verbose::Bool,
    history::Dict{Int,Tuple{Int,Int}},
    first_guess::Function = guess_any,
)::Int
    all_codes = collect(0:9999)            # randomized list of possible answers
    shuffle!(all_codes)

    steps = 0

    g = first_guess()
    cnt = result_server(g)
    history[g] = cnt
    steps += 1
    if verbose
        println(
            "Guess $steps:\t",
            string4(g),
            "\t",
            cnt,
            "\tsearch space size is now ",
            history_count(history),
        )
    end

    if cnt == (4, 0)
        if verbose
            println("Solved! Code is $(string4(g))")
        end
        return steps
    end

    for g in all_codes  # current guess (g) code

        if !history_check(history, g)
            continue
        end

        cnt = result_server(g)

        steps += 1
        history[g] = cnt

        if verbose
            println(
                "Guess $steps:\t",
                string4(g),
                "\t",
                cnt,
                "\tsearch space size is now ",
                history_count(history),
            )
        end

        if cnt == (4, 0)
            if verbose
                println("Solved! Code is $(string4(g))")
            end
            return steps
        end
    end
    return -1  # signal failure
end


text:
module ShaderManager

using ..GraphicsManager

API = nothing

function init()
	global API = GRAPHICSDRIVER().ShaderManager
end

create(name::Symbol, path::String, listen::Union{Void,Function}=nothing) = API.create(name, path, listen)

setListenerOnShaderPropertyUpdate(name::Symbol, listen::Tuple) = API.setListenerOnShaderPropertyUpdate(name, listen)

getShaderProperties(program::AbstractGraphicsShaderProgram, buffer::AbstractGraphicsData) = API.getShaderProperties(program, buffer)
render(buffer::AbstractGraphicsData, data::AbstractGraphicsData) = API.render(buffer, data)

watchShaderReload(f::Function) = API.watchShaderReload(f)
linkTo(program) = API.linkTo(program)

getCurrent() = API.getCurrent()
reset() = API.reset()
unlink() = API.unlink()

end #ShaderManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
@testset "memory" begin

let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# CUDA.available_memory()
    #=@test b == =# CUDA.total_memory()
end

# dummy data
T = UInt32
N = 5
data = rand(T, N)
nb = sizeof(data)

# buffers are untyped, so we use a convenience function to get a typed pointer
# we prefer to return a device pointer (for managed buffers) to maximize CUDA coverage
typed_pointer(buf::Union{Mem.Device, Mem.Unified}, T) = convert(CuPtr{T}, buf)
typed_pointer(buf::Mem.Host, T)                       = convert(Ptr{T},   buf)

# allocations and copies
for srcTy in [Mem.Device, Mem.Host, Mem.Unified],
    dstTy in [Mem.Device, Mem.Host, Mem.Unified]

    dummy = Mem.alloc(srcTy, 0)
    Mem.free(dummy)

    src = Mem.alloc(srcTy, nb)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N)

    dst = Mem.alloc(dstTy, nb)
    unsafe_copyto!(typed_pointer(dst, T), typed_pointer(src, T), N)

    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)

    @test data == ref

    if isa(src, Mem.Device) || isa(src, Mem.Unified)
        Mem.set!(typed_pointer(src, T), zero(T), N)
    end

    # test the memory-type attribute
    if isa(src, Mem.Device)
        @test CUDA.memory_type(typed_pointer(src, T)) == CUDA.MEMORYTYPE_DEVICE
    elseif isa(src, Mem.Host)
        @test CUDA.memory_type(convert(Ptr{T}, src)) == CUDA.MEMORYTYPE_HOST
    elseif isa(src, Mem.Unified)
        # unified memory can reside in either place
        # FIXME: does this depend on the current migration, or on the configuration?
        @test CUDA.memory_type(convert(CuPtr{T}, src)) == CUDA.MEMORYTYPE_HOST ||
              CUDA.memory_type(convert(CuPtr{T}, src)) == CUDA.MEMORYTYPE_DEVICE ||
        @test CUDA.memory_type(convert(CuPtr{T}, src)) == CUDA.memory_type(convert(Ptr{T}, src))
    end

    # test the is-managed attribute
    if isa(src, Mem.Unified)
        @test CUDA.is_managed(convert(Ptr{T}, src))
        @test CUDA.is_managed(convert(CuPtr{T}, src))
    else
        @test !CUDA.is_managed(typed_pointer(src, T))
    end
    # Test conversion to Ptr throwing an error
    if isa(src, Mem.Device)
        @test_throws ArgumentError convert(Ptr, src)
    end

    Mem.free(src)
    Mem.free(dst)
end

# pointer attributes
let
    src = Mem.alloc(Mem.Device, nb)

    attribute!(typed_pointer(src, T), CUDA.POINTER_ATTRIBUTE_SYNC_MEMOPS, 0)

    Mem.free(src)
end

# asynchronous operations
let
    src = Mem.alloc(Mem.Device, nb)

    @test_throws ArgumentError unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true, stream=CuDefaultStream())

    Mem.set!(typed_pointer(src, T), zero(T), N; async=true, stream=CuDefaultStream())

    Mem.free(src)
end

# pinned memory
let
    # can only get GPU pointer if the pinned buffer is mapped
    src = Mem.alloc(Mem.Host, nb)
    @test_throws ArgumentError convert(CuPtr{T}, src)
    Mem.free(src)

    # create a pinned and mapped buffer
    src = Mem.alloc(Mem.Host, nb, Mem.HOSTALLOC_DEVICEMAP)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    @test CUDA.memory_type(cpu_ptr) == CUDA.MEMORYTYPE_HOST
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    @test CUDA.memory_type(gpu_ptr) == CUDA.MEMORYTYPE_HOST
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.free(src)
    # NOTE: don't free dst, it's just a mapped pointer
end

# pinned memory with existing memory
if attribute(device(), CUDA.DEVICE_ATTRIBUTE_HOST_REGISTER_SUPPORTED) != 0
    # can only get GPU pointer if the pinned buffer is mapped
    @test_throws ArgumentError Mem.register(Mem.Host, pointer(data), 0)
    src = Mem.register(Mem.Host, pointer(data), nb)
    @test_throws ArgumentError convert(CuPtr{T}, src)
    Mem.unregister(src)

    # register a pinned and mapped buffer
    src = Mem.register(Mem.Host, pointer(data), nb, Mem.HOSTREGISTER_DEVICEMAP)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.unregister(src)
    # NOTE: don't unregister dst, it's just a mapped pointer
end

# unified memory
let
    src = Mem.alloc(Mem.Unified, nb)

    @test_throws BoundsError Mem.prefetch(src, 2*nb; device=CUDA.DEVICE_CPU)
    # FIXME: prefetch doesn't work on some CI devices, unsure why.
    @test_skip Mem.prefetch(src, nb; device=CUDA.DEVICE_CPU)
    Mem.advise(src, Mem.ADVISE_SET_READ_MOSTLY)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.free(src)
end

# 3d memcpy
let
    # TODO: use cuMemAllocPitch (and put pitch in buffer?) to actually get benefit from this

    data = collect(reshape(1:27, 3, 3, 3))

    dst = Mem.alloc(Mem.Device, sizeof(data))
    Mem.unsafe_copy3d!(typed_pointer(dst, Int), Mem.Device, pointer(data), Mem.Host, length(data))

    check = zeros(Int, size(data))
    Mem.unsafe_copy3d!(pointer(check), Mem.Host, typed_pointer(dst, Int), Mem.Device, length(data))

    @test check == data

    Mem.free(dst)
end
let
    # copying an x-z plane of a 3-D array

    T = Int
    nx, ny, nz = 4, 4, 4
    data = collect(reshape(1:(nx*nz), nx, nz))
    dst = Mem.alloc(Mem.Device, nx*ny*nz*sizeof(data))

    # host to device
    Mem.unsafe_copy3d!(typed_pointer(dst, T), Mem.Device, pointer(data), Mem.Host,
                       nx, 1, nz;
                       dstPos=(1,2,1),
                       srcPitch=nx*sizeof(T), srcHeight=1,
                       dstPitch=nx*sizeof(T), dstHeight=ny)

    # copy back
    check = zeros(T, size(data))
    Mem.unsafe_copy3d!(pointer(check), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       nx, 1, nz;
                       srcPos=(1,2,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=1)

    @test all(check .== data)

    # copy back into a 3-D array
    check2 = zeros(T, nx, ny, nz)
    Mem.unsafe_copy3d!(pointer(check2), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       nx, 1, nz;
                       srcPos=(1,2,1),
                       dstPos=(1,2,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=ny)
    @test all(check2[:,2,:] .== data)
end

let
    # copying an y-z plane of a 3-D array

    T = Int
    nx, ny, nz = 4, 4, 4
    data = collect(reshape(1:(ny*nz), ny, nz))
    dst = Mem.alloc(Mem.Device, nx*ny*nz*sizeof(data))

    # host to device
    Mem.unsafe_copy3d!(typed_pointer(dst, T), Mem.Device, pointer(data), Mem.Host,
                       1, ny, nz;
                       dstPos=(2,1,1),
                       srcPitch=1*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=ny)

    # copy back
    check = zeros(T, size(data))
    Mem.unsafe_copy3d!(pointer(check), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       1, ny, nz;
                       srcPos=(2,1,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=1*sizeof(T), dstHeight=ny)

    @test all(check .== data)

    # copy back into a 3-D array
    check2 = zeros(T, nx, ny, nz)
    Mem.unsafe_copy3d!(pointer(check2), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       1, ny, nz;
                       srcPos=(2,1,1),
                       dstPos=(2,1,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=ny)
    @test all(check2[2,:,:] .== data)
end

@testset "pin" begin
    buf = Array{UInt8}(undef, 500)
    Mem.pin(view(buf, 100:200))
    Mem.pin(view(buf, 50:250))   # should subsume the previous one
    Mem.pin(view(buf, 225:230))  # shouldn't do anything
    Mem.pin(view(buf, 250:300))  # should extend
    Mem.pin(view(buf, 25:50))    # should extend
    Mem.pin(view(buf, 1:100))    # should partially extend
    Mem.pin(view(buf, 250:350))  # should partially extend
    Mem.pin(view(buf, 1:350))    # do nothing
    Mem.pin(buf)
end

end


text:
using OrdinaryDiffEq
using GLMakie
using ForwardDiff

#=

    2. Higher dimensions

    We consider the first order ODE
    in ℝⁿ × ℝᵐ where we have

    dS/dt = -S*I
    dI/dt = S*I - I
    dR/dt = I
=#

function lorenz!(du, u, p, t)
     du[1] = 10.0*(u[2]-u[1])
     du[2] = u[1]*(28.0-u[3]) - u[2]
     du[3] = u[1]*u[2] - (8/3)*u[3]
    nothing  # no need to return anything
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)

prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob, saveat = 0.01)

begin
    f = Figure()

    ax = Axis3(f[1,1], title = "Lorenz model")
    lines!(sol[1,:], sol[2,:], sol[3,:])

    ax.elevation = pi/8
    ax.azimuth = -1pi/8

    f
end

display(f)

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
@testset "Meshing" begin
    c = Mesh(translate([4,5,6])Cuboid([1,2,3]))
    @test_broken length(c.vertices) == 7638
    @test_broken length(c.faces) == 15272

    @testset "Radiused Shelled Box" begin
        c2 = Cuboid([4,4.0,4])
        cyl1 = rotate(pi/6, [1,0,0])translate([2,2,0])Cylinder(1.0,10.0)

        u = Shell(0.5)RadiusedCSGUnion(1,c2, cyl1)
        u2 = CSGDiff(u, Cuboid([2,2,2]))

        m = Mesh(u2)
        @test_broken length(m.vertices) == 86406
        @test_broken length(m.faces) == 172678
    end
end


text:
using Test
using Clp
using LinearAlgebra
using SteinDiscrepancy:
    SteinInverseMultiquadricKernel,
    SteinRectangularDomainMetaKernel,
    stein_discrepancy,
    ksd,
    gsd,
    wassersteindiscrete

# import data
include("data/test_data.jl")

function uniform_gradlogp(x::Array{Float64,1})
    zeros(size(x))
end

function gaussian_gradlogp(x::Array{Float64,1})
    -x
end
# define the covariance terms for the Riemannian-Langevin diffusion
function volatility_covariance(x::Array{Float64,1})
    p = length(x)
    sqrt(1 + sum(x.^2)) * I
end
function grad_volatility_covariance(x::Array{Float64,1})
    x / sqrt(1 + sum(x.^2))
end

imqkernel = SteinInverseMultiquadricKernel()
imqrectkernel = SteinRectangularDomainMetaKernel(imqkernel, [0., 0.], [1., 1.])
solver = Clp.ClpSolver(LogLevel=4)

# Graph Stein discrepancy bounded test
@testset "Univariate Graph discrepancy test" begin
    res = gsd(points=UNIFORM_TESTDATA[:,1],
              gradlogdensity=uniform_gradlogp,
              solver=solver,
              supportlowerbounds=[0.0],
              supportupperbounds=[1.0])

    @test res.objectivevalue[1] ≈ 0.15 atol=1.0e-5
end
@testset "Multivariate Graph discrepancy test" begin
    res = stein_discrepancy(points=UNIFORM_TESTDATA,
                            gradlogdensity=uniform_gradlogp,
                            solver=solver,
                            method="graph",
                            supportlowerbounds=zeros(size(UNIFORM_TESTDATA,2)),
                            supportupperbounds=ones(size(UNIFORM_TESTDATA,2)))

    @test res.objectivevalue[1] ≈ 0.308 atol=1.0e-5
    @test res.objectivevalue[2] ≈ 0.31 atol=1.0e-5
end
@testset "Riemannian Graph discrepancy test" begin
    res = stein_discrepancy(points=GAUSSIAN_TESTDATA,
                            gradlogdensity=gaussian_gradlogp,
                            operator="riemannian-langevin",
                            solver=solver,
                            method="graph",
                            volatility_covariance=volatility_covariance,
                            grad_volatility_covariance=grad_volatility_covariance)

    @test res.objectivevalue[1] ≈ 6.19106 atol=1.0e-5
    @test res.objectivevalue[2] ≈ 7.05982 atol=1.0e-5
end

# Kernel Stein discrepancy test
@testset "Univariate Kernel discrepancy test" begin
    res = stein_discrepancy(points=GAUSSIAN_TESTDATA[:,1],
                            gradlogdensity=gaussian_gradlogp,
                            method="kernel",
                            kernel=imqkernel)

    @test res.discrepancy2 ≈ 0.269298 atol=1.0e-5
end
@testset "Multivariate Kernel discrepancy test" begin
    res = ksd(points=GAUSSIAN_TESTDATA,
              gradlogdensity=gaussian_gradlogp,
              kernel=imqkernel)

    @test res.discrepancy2 ≈ 0.807566 atol=1.0e-5
end

# Wassersteindiscrete test
@testset "Wasserstein bivariate discrete test" begin
    (emd, numnodes, numedges, status) =
        wassersteindiscrete(xpoints=GAUSSIAN_TESTDATA,
                            ypoints=UNIFORM_TESTDATA,
                            solver=solver)

    @test emd ≈ 1.56 atol=1e-5
    @test status == :Optimal
end

# Stein Meta Kernel test
@testset "Rectangular domain meta kernel test" begin
    res = ksd(points=UNIFORM_TESTDATA,
              gradlogdensity=uniform_gradlogp,
              kernel=imqrectkernel)

    @test res.discrepancy2 ≈ 0.0064761153 atol=1.0e-8
end


text:
# Routines to compute surface velocities and their jumps

@inline function surface_velocity_jump!(Δus::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    @unpack surfacevel_cache, motions, bodies = sys
    @unpack Vb = surfacevel_cache
    assign_velocity!(Vb,bodies,motions,t)
    surface_velocity_jump!(Δus,Vb,sys)
    return Δus
  end

@inline surface_velocity_jump!(Δus::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,InternalFlow}) where {NX,NY,N,MT,FS} =
                               (Δus .= -us; Δus)

@inline surface_velocity_jump!(Δus::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,ExternalFlow}) where {NX,NY,N,MT,FS} =
                               (Δus .= us; Δus)

@inline surface_velocity_jump!(Δus::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,ExternalInternalFlow}) where {NX,NY,N,MT,FS} =
                               (Δus .= 0.0; Δus)

@inline function surface_velocity!(ūs::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    @unpack surfacevel_cache, motions, bodies = sys
    @unpack Vb = surfacevel_cache
    assign_velocity!(Vb,bodies,motions,t)
    surface_velocity!(ūs,Vb,sys)
    return ūs
  end

@inline surface_velocity!(ūs::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,ExternalInternalFlow}) where {NX,NY,N,MT,FS} =
                               (ūs .= us; ūs)

@inline surface_velocity!(ūs::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,SD}) where {NX,NY,N,MT,FS,SD} =
                               (ūs .= 0.5*us; ūs)

@inline function relative_surface_velocity!(ūsr::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
  @unpack surfacevel_cache, motions, bodies = sys
  @unpack Vb = surfacevel_cache
  assign_velocity!(Vb,bodies,motions,t)
  surface_velocity!(ūsr,Vb,sys)
  ūsr .-= Vb
end


text:
export
    Detection,
    check!, update!, set_device_point!, draw, highlight_color

using Printf

# Detection is for detecting if a device point is within
# the node's polygon.
mutable struct Detection
    local_point::Geometry.Point{Float64}
    device_point::Geometry.Point{Float64}
    inside_node::Bool
    inside_color::Rendering.Palette
    outside_color::Rendering.Palette

    function Detection(inside_color::Rendering.Palette, outside_color::Rendering.Palette)
        o = new()

        o.device_point = Geometry.Point{Float64}()
        o.local_point = Geometry.Point{Float64}()
        o.inside_node = false
        o.inside_color = inside_color
        o.outside_color = outside_color

        o
    end
end

# Typically called from a Nodes.draw
function check!(detection::Detection, node::Ranger.AbstractNode, world::Ranger.World)
    # Map device/mouse coords to local-space of node.
    Nodes.map_device_to_node!(world,
        Int32(detection.device_point.x), Int32(detection.device_point.y),
        node, detection.local_point)

    detection.inside_node
end

function highlight_color(detection::Detection, inside::Bool)
    if inside
        return detection.inside_color
    else
        return detection.outside_color
    end
end

# Typically called from a Nodes.update
function update!(detection::Detection, node::Ranger.AbstractNode)
    detection.inside_node = Geometry.is_point_inside(node.polygon, detection.local_point)
end

function update!(detection::Detection, polygon::Geometry.Polygon)
    detection.inside_node = Geometry.is_point_inside(polygon, detection.local_point)
end

# Typically called from a Nodes.io_event
function set_device_point!(detection::Detection, x::Float64, y::Float64)
    Geometry.set!(detection.device_point, x, y)
end

# Draw debug local-space coords
function draw(detection::Detection, context::Rendering.RenderContext)
    color = if detection.inside_node
        detection.inside_color
    else
        detection.outside_color
    end
    Rendering.set_draw_color(context, color)
    text = @sprintf("L: %2.4f, %2.4f", detection.local_point.x, detection.local_point.y)
    Rendering.draw_text(context, 10, 70, text, 2, 2, false)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module BeliefUpdaters

using RPOMDPs
import RPOMDPs: Updater, update, initialize_belief, pdf, mode, updater, iterator
import Base: rand, rand!, mean, ==

using RPOMDPToolbox.ordered_states
using StatsBase

export
    VoidUpdater
include("void.jl")

export
    DiscreteBelief,
    DiscreteUpdater,
    uniform_belief,
    product                     # Remove because deprecated
include("discrete.jl")

export
    PreviousObservationUpdater,
    FastPreviousObservationUpdater,
    PrimedPreviousObservationUpdater

include("previous_observation.jl")

end


text:
abstract type PipelineStep end

"""
    run(step::PipelineStep, context, data)

Applies the operation `step` to `data`
"""
function run end

"""
    run!(buf, step::PipelineStep, context, data)

Applies the operation `step` inplace to `buf`. `buf` is mutated.
"""
function run! end

"""
    invert(step::PipelineStep, context, data)

Applies the inverse of the operation `step` to `data`
"""
function invert end

"""
    invert!(buf, step::PipelineStep, context, data)

Applies the inverse of the operation `step` to `buf` inplace. `buf` is mutated,
"""
function invert! end


text:
function fixed_R_YSZ()
  return CubicSplineInterpolation(
        600 : 50 : 800
      , 1 ./reverse([
        0.044905404632136
        0.03005779931599
        0.01725995101764
        0.010115351717062
        0.005578168458813
      ])
    )
end

# temperature data are in degrees of Celsia
# resistanace data are in Ohm/cm
# ... a bit different interpretation of porosity
function TI_clank(label)
  if label == "R_LSM"

    return CubicSplineInterpolation(
        600 : 50 : 800
      ,
      # fitted to 3D with "clankova" porozita
      [0.004640194824206025, 0.0040162711411057365, 0.0033223061268364204, 0.0024328408017493925, 0.0017474026238125484]
      # fitted to 2D
      #[0.005764, 0.005022, 0.00413, 0.003021, 0.002185]
    )
  elseif label == "R_YSZ"
    # ok data from pure YSZ (porosity = 0.0)
    return fixed_R_YSZ()
 end
end

function TI(label)
  if label == "R_LSM"

    return CubicSplineInterpolation(
        600 : 50 : 800
      ,
      # fitted to 3D
      [0.006418221532022811, 0.005556938763855908, 0.004598968863545539, 0.00337112453237228, 0.0024249253299838836]
      # fitted to 2D
      #[0.005764, 0.005022, 0.00413, 0.003021, 0.002185]
    )
  elseif label == "R_YSZ"
    # ok data from pure YSZ (porosity = 0.0)
    return fixed_R_YSZ()
 end
end

function TI_2D(label)
  if label == "R_LSM"

    return CubicSplineInterpolation(
        600 : 50 : 800
      ,
      # fitted to 3D
      #[0.006418221532022811, 0.005556938763855908, 0.004598968863545539, 0.00337112453237228, 0.0024249253299838836]
      # fitted to 2D
      [0.005764, 0.005022, 0.00413, 0.003021, 0.002185]
    )
  elseif label == "R_YSZ"
    # ok data from pure YSZ (porosity = 0.0)
    return fixed_R_YSZ()
 end
end

function TI_por_LSM(T)
  #0.12 porosity LSM
  return CubicSplineInterpolation(
      600 : 50 : 800
    , 1 ./reverse([
    286.138977466499
    205.521051620173
    150.497871331355
    124.493586820521
    107.754096313306
    ])
  )(T)
end

TI(label, T) = TI(label)(T)
TI_2D(label, T) = TI_2D(label)(T)
TI_clank(label, T) = TI_clank(label)(T)

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


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

using Markdown
using InteractiveUtils

# ╔═╡ 34d2b0f6-ff65-11ea-05bc-c1a781f7bd7d
using Pkg, DrWatson

# ╔═╡ 34d2e742-ff65-11ea-35e7-6dba8cede1d5
begin
  using Distributions
  using StatsPlots
  using StatsBase
  using LaTeXStrings
  using CSV
  using DataFrames
  using LinearAlgebra
  using Random
	using StatisticalRethinking
end

# ╔═╡ a63d2808-ff64-11ea-037d-f58c33309a6d
md"## Clip-09-01s.jl"

# ╔═╡ 34d76984-ff65-11ea-36f9-c1fe6d3edf77
md"### snippet 9.1"

# ╔═╡ 34e18818-ff65-11ea-327d-37ab359d53e7
md"##### Metropolis algorithm."

# ╔═╡ 34e21e60-ff65-11ea-12d2-bd654991499b
function generate_walk(N::Int64)
  num_weeks = N
  positions = zeros(Int64, num_weeks);
  current = 10
  d = Uniform(0, 1)

  for i in 1:N
    positions[i] = current  # Record current position
    proposal = current + sample([-1, 1], 1)[1] # Generate proposal
    proposal = proposal < 1  ? 10 : proposal
    proposal = proposal > 10  ? 1 : proposal
    prob_move = proposal/current  # Move?
    current = rand(d, 1)[1] <  prob_move ? proposal : current
  end

  positions

end

# ╔═╡ 34ee57d4-ff65-11ea-10a4-e3ee18d27928
md"##### Generate the visits ."

# ╔═╡ 34eefa9a-ff65-11ea-269d-3b6d166ebf54
begin
	N = 100000
	walk = generate_walk(N)
end;

# ╔═╡ 34f96ab6-ff65-11ea-3f24-bfe0191b164d
# Plot the first 100 weeks and a histogram of weeks per iosland

begin
	fig = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
	fig[1] = plot(walk[1:100], leg=false, xlabel="Week",
		ylabel="Island", title="First 100 steps")
	fig[2] = histogram(walk, leg=false, xlabel="Island", ylabel="Number of weeks",
	  title="$N steps")
	plot(fig..., layout=(1, 2))
end

# ╔═╡ 34fa3306-ff65-11ea-395b-5b35170d0099
md"## End of clip-09-01s.jl"

# ╔═╡ Cell order:
# ╟─a63d2808-ff64-11ea-037d-f58c33309a6d
# ╠═34d2b0f6-ff65-11ea-05bc-c1a781f7bd7d
# ╠═34d2e742-ff65-11ea-35e7-6dba8cede1d5
# ╟─34d76984-ff65-11ea-36f9-c1fe6d3edf77
# ╟─34e18818-ff65-11ea-327d-37ab359d53e7
# ╠═34e21e60-ff65-11ea-12d2-bd654991499b
# ╟─34ee57d4-ff65-11ea-10a4-e3ee18d27928
# ╠═34eefa9a-ff65-11ea-269d-3b6d166ebf54
# ╠═34f96ab6-ff65-11ea-3f24-bfe0191b164d
# ╟─34fa3306-ff65-11ea-395b-5b35170d0099


text:
include("test_common.jl")
using Random
using Suppressor

@testset "executables" begin
    @test read(node`--version`, String) |> strip == "v14.17.3"
    @test read(npm`--version`, String) |> strip == "6.14.13"
    @test @capture_out(begin
        npx("--version")
    end)|> strip == "6.14.13"
    ps = [
        getproperty(NPM, Symbol(replace(command, '-'=>'_')))("--help"; wait=false)
        for command in NPM.COMMANDS
    ]
    wait.(ps)
    @suppress NPM.init("-y")
    @suppress NPM.install("canvas") do p
        wait(p)
    end
    @test NPM.is_installed("canvas")
    @test NPM.is_installed("canvas"; is_global=true) isa Bool
end


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rastername, rasterpath, rasterurl

@testset "EarthEnv HabitatHeterogeneity" begin
    using RasterDataSources: rastername, rasterurl, rasterpath

    @test rastername(EarthEnv{HabitatHeterogeneity}, :Dissimilarity; res="1km") == "Dissimilarity_1km.tif"
    hh_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "EarthEnv", "HabitatHeterogeneity")
    @test rasterpath(EarthEnv{HabitatHeterogeneity}) == hh_path
    @test rasterpath(EarthEnv{HabitatHeterogeneity}, :Dissimilarity; res="1km") == joinpath(hh_path, "1km", "Dissimilarity_1km.tif")

    @test rasterurl(EarthEnv{HabitatHeterogeneity}, :Dissimilarity; res="1km") ==
        URI(scheme="https", host="data.earthenv.org", path="/habitat_heterogeneity/1km/Dissimilarity_01_05_1km_uint32.tif")
    raster_path = joinpath(hh_path, "25km", "Dissimilarity_25km.tif")
    @test getraster(EarthEnv{HabitatHeterogeneity}, :dissimilarity; res="25km") == raster_path
    @test getraster(EarthEnv{HabitatHeterogeneity}, (:Dissimilarity,)) == (dissimilarity=raster_path,)
    @test getraster(EarthEnv{HabitatHeterogeneity}, [:dissimilarity]) == (dissimilarity=raster_path,)
    @test isfile(raster_path)
    files = getraster(EarthEnv{HabitatHeterogeneity})
    @test all(map(isfile, files))
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
import Distributions
import LinearAlgebra
import Statistics
import StatsBase

ncol(df::DataFrames.AbstractDataFrame) = size(df, 2)

ncol(m::AbstractMatrix) = size(m, 2)

function numeric_columns(
        dat,
        )
    p = ncol(dat)
    is_numeric = falses(p)
    for j = 1:p
        typ = eltype(dat[:, j])
        if typ <: Real && typ ≠ Bool
            is_numeric[j] = true
        end
    end
    res = findall(is_numeric)
    res
end

function fill_diagonal!(
        X,
        diag_elems,
        )
    p = size(X, 2)
    for i = 1:p
        X[i, i] = diag_elems[i]
    end
end

function rose_real(
        X,
        n,
        ids_class,
        ids_generation,
        h_mult = 1,
        )
    X_array = convert(Array, X)
    p = size(X_array, 2)
    n_new = length(ids_generation)
    cons_kernel = (4/((p+2) * n))^(1/(p+4))

    if p ≠ 1
        # sd_mat = eye(p)
        sd_mat = Matrix(1.0LinearAlgebra.I, p, p)
        sd_vect = Statistics.std(X_array[ids_class, :];dims = 2,)
        fill_diagonal!(sd_mat, sd_vect)
        H = h_mult * cons_kernel * sd_mat
    else
        H = h_mult * cons_kernel * Statistics.std(X_array[ids_class, :])
    end
    X_new_num = randn(n_new, p) * H
    X_new_num = X_new_num + X_array[ids_generation, :]
    return X_new_num
end

function rose_sampling(
        X,
        y,
        prop,
        indcs_maj,
        indcs_min,
        y_majority,
        y_minority,
        h_mult_maj,
        h_mult_min,
        )
    n = size(X, 1)
    n_minority = sum(rand(Distributions.Binomial(1, prop), n))
    n_majority = n - n_minority

    indcs_maj_new = StatsBase.sample(indcs_maj, n_majority, replace = true)
    indcs_min_new = StatsBase.sample(indcs_min, n_minority, replace = true)

    numeric_cols = numeric_columns(X)

    # Create  X
    indcs = vcat(indcs_maj_new, indcs_min_new)
    X_new = X[indcs, :]
    if length(numeric_cols) > 0
        majority_rose_real_results = rose_real(
            X[:, numeric_cols],
            length(indcs_maj),
            indcs_maj,
            indcs_maj_new,
            h_mult_maj,
            )
        for k = 1:length(numeric_cols)
            col = numeric_cols[k]
            X_new[1:n_majority, col] = majority_rose_real_results[:, k]
        end
        minority_rose_real_results = rose_real(
            X[:, numeric_cols],
            length(indcs_min),
            indcs_min,
            indcs_min_new,
            h_mult_min,
            )
        for k = 1:length(numeric_cols)
            col = numeric_cols[k]
            X_new[(n_majority + 1):n, col] = minority_rose_real_results[:, k]
        end
    end

    # Create y
    y_new = similar(y)
    y_new[1:n_majority] .= y_majority
    y_new[(n_majority + 1):n] .= y_minority

    result = (X_new, y_new,)
    return result
end

"""
    classlabel(y)
Given a column from a DataFrames.DataFrame, this function returns the majority/minority class label.
"""
function classlabel(
        y::Array{T, 1},
        labeltype = :minority,
        ) where T
    count_dict = StatsBase.countmap(y)
    labels = collect(keys(count_dict))
    counts = collect(values(count_dict))
    if length(labels) > 2
        error("There are more than two classes in the target variable.")
    elseif length(labels) < 2
        error("There is only one class in the target variable.")
    end
    func = (labeltype == :majority) ? :argmax : :argmin
    indx = eval(func)(counts)
    res = collect(labels)[indx]
    res
end

function rose(
        dat::DataFrames.DataFrame,
        y_column::Symbol,
        prop::Float64 = 0.5,
        h_mult_maj = 1,
        h_mult_min = 1,
        )
    majority_label = classlabel(dat[y_column], :majority)
    minority_label = classlabel(dat[y_column], :minority)

    indcs_maj = findall(dat[y_column] .== majority_label)
    indcs_min = findall(dat[y_column] .== minority_label)
    p = size(dat, 2)
    y_indx = findfirst(names(dat) .== y_column)
    X_indcs = setdiff(1:p, y_indx)
    X = dat[:, X_indcs]
    y = dat[y_column]
    y_majority = majority_label
    y_minority = minority_label
    X_new, y_new = rose_sampling(
        X,
        y,
        prop,
        indcs_maj,
        indcs_min,
        y_majority,
        y_minority,
        h_mult_maj,
        h_mult_min,
        )
    result = X_new
    result[y_column] = y_new
    return result
end


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

using ValueShapes
using Test

using Distributions

@testset "distributions" begin
    @test @inferred(varshape(Normal())) == ScalarShape{Real}()
    @test @inferred(varshape(MvNormal([2. 1.; 1. 3.]))) == ArrayShape{Real}(2)
    @test @inferred(varshape(MatrixBeta(4, 6, 6))) == ArrayShape{Real}(4, 4)

    @test @inferred(vardof(MvNormal([2. 1.; 1. 3.]))) == 2
end


text:
function plot_basic_state(y, U, E, file)

    p1 = plot(
        y,
        U,
        linewidth = 3,
        label = "U",
        title = "Velocity",
        xlims = (y[end], y[1])
        )

    p2 = plot(
        y,
        E,
        linewidth = 3,
        label = "η",
        xlabel = "y",
        title = "Free-Surface",
        xlims = (y[end], y[1])
        )

    plt1 = plot(p1, p2, layout = (2,1))
    savefig(plt1, file)

end

function plot_growth_rates(ks, σ, Nmodes, file)

    titles = ["1st mode", "2nd mode"]
    plt2 = plot()
    for cnt in 1:Nmodes
        plot!(plt2,
        ks,
        σ[cnt,:],
        linewidth = 3,
        label=titles[cnt],
        title="Growth Rates",
        xlabel="k",
        xlims=(ks[1], ks[end])
        )
    end

    savefig(plt2,file)

end

function plot_1D_streamfunction(k_index, k, y, σmodes, mode_number, file)

      Ny = length(y) - 1
      Lx = 2*pi/k;
      Nx = Ny;
      dx = Lx/Nx;
       x = collect(0:dx:Lx);
    X, Y = meshgrid(x,y)

    # Divide eigenvector matrix into specified spaces: u,v,eta.
    uvec = σmodes[1:Ny+1,        mode_number, k_index];
    vvec = σmodes[Ny+2:2*Ny+2,   mode_number, k_index];
    ηvec = σmodes[2*Ny+1:3*Ny+1, mode_number, k_index];

    p1 = plot(
        y,
        real(uvec),
        linewidth = 3,
        label = "real",
        xlabel = "ϕ",
        title = "u",
        xlims = (y[end], y[1])
        )
    plot!(
        p1,
        y,
        linewidth = 3,
        imag(uvec),
        label = "imag")
    p2 = plot(
        y,
        real(vvec),
        linewidth = 3,
        label = "real",
         title = "v",
        xlims = (y[end], y[1])
        )
    plot!(
        p2,
        y,
        linewidth = 3,
        imag(vvec),
        label = "imag")
    p3 = plot(
        y,
        real(ηvec),
        linewidth = 3,
        label = "real",
         title = "η",
        xlims = (y[end], y[1])
        )
    plot!(
        p3,
        y,
        linewidth = 3,
        imag(vvec),
        label = "imag")
    plt = plot(p1, p2, p3, layout = (3,1))
    savefig(plt, file)
    #display(plt)

end

function plot_2D_streamfunction(k_index, k, y, σmodes, mode_number, file)

    Ny = length(y) - 1
    Lx = 2*pi/k;
    Nx = Ny;
    dx = Lx/Nx;
     x = collect(0:dx:Lx+dx/2);
  X, Y = meshgrid(x,y[end:-1:1])

    uvec = σmodes[1:Ny+1,        mode_number, k_index];
    vvec = σmodes[Ny+2:2*Ny+2,   mode_number, k_index];
    ηvec = σmodes[2*Ny+1:3*Ny+1, mode_number, k_index];

    u =     repeat(real(uvec),1,Nx+1).*cos.(k*X) -    repeat(imag(uvec),1,Nx+1).*sin.(k*X);
    v = -k.*repeat(imag(vvec),1,Nx+1).*cos.(k*X) - k.*repeat(real(vvec),1,Nx+1).*sin.(k*X);
    η =     repeat(real(ηvec),1,Nx+1).*cos.(k*X) -    repeat(imag(ηvec),1,Nx+1).*sin.(k*X);

    kwargs = (
        xlabel = "x",
        ylabel = "y",
          fill = true,
        levels = 20,
     linewidth = 0,
         color = :balance,
      colorbar = true,
          xlim = (    0, Lx),
         ylim = (y[end], y[1])
    )

    u_plt = contour(x, y[end:-1:1], u[:, end:-1:1], title="u"; kwargs...)
    v_plt = contour(x, y[end:-1:1], v[:, end:-1:1], title="v"; kwargs...)
    η_plt = contour(x, y[end:-1:1], η[:, end:-1:1], title="η"; kwargs...)

    plt1 = plot(u_plt, v_plt, η_plt, layout = (1,3))
    savefig(plt1, file)

end

function plot_2D_vorticity(Dy, k_index, k, y, σmodes, mode_number, file)

    Ly = params.Lϕ
    Ny = length(y) - 1
    Lx = 2*pi/k;
    Nx = Ny;
    dx = Lx/Nx;
     x = collect(0:dx:Lx+dx/2);
  X, Y = meshgrid(x,y[end:-1:1])

    uvec = σmodes[1:Ny+1,        mode_number, k_index];
    vvec = σmodes[Ny+2:2*Ny+2,   mode_number, k_index];

    ζvec = im * k * vvec - Dy * uvec
    ζ    = repeat(real(ζvec),1,Nx+1).*cos.(k*X) - repeat(imag(ζvec),1,Nx+1).*sin.(k*X);

    kwargs = (
        xlabel = "x",
        ylabel = "y",
          fill = true,
        levels = 20,
     linewidth = 0,
         color = :balance,
      colorbar = true,
          xlim = (    0, Lx),
         ylim = (y[end], y[1])
    )

    plt = contour(x, y[end:-1:1], ζ[:, end:-1:1], title="ζ"; kwargs...)

    savefig(plt, file)

end


text:
export SetAtom, contains, is_contained_in, is_disjoint_from, intersects
import .LazySets: dim, project

@static if VERSION >= v"1.5"
    import Base: contains
end

using .LazySets: LazySet, ⊆, isdisjoint

"""
    SetAtom{S<:LazySet, T} <: Predicate{Val{1}}

A unary atomic predicate defined with respect to a set.

### Fields

- `X` -- set
- `f` -- function that takes the set `X` as first argument
"""
struct SetAtom{S<:LazySet, T} <: Predicate{Val{1}}
    X::S
    f::T

    function SetAtom(X::S, f::T) where {S<:LazySet, T}
        return new{S, T}(X, f)
    end
end

# function-like evaluation
@inline function (sa::SetAtom)(args...)
    evaluate(sa, args...)
end

function evaluate(sa::SetAtom, args...)
    assert_same_length(sa, args...)
    return sa.f(sa.X, args...)
end

function Base.:(==)(sa1::SetAtom, sa2::SetAtom)
    return sa1.X == sa2.X && sa1.f == sa2.f
end

# ===========================
# Dimension of set predicates
# ===========================

function dim(sa::SetAtom)
    return dim(sa.X)
end

# fallback
function dim(p::Predicate)
    throw(ArgumentError("`dim` cannot be applied to a `$(typeof(p))`; use a " *
                        "`SetAtom` instead"))
end

# ============================
# Projection of set predicates
# ============================

function project(sa::SetAtom, vars::AbstractVector{Int})
    X_proj = project(sa.X, vars)
    return SetAtom(X_proj, sa.f)
end

# fallback
function project(p::Predicate, vars::AbstractVector{Int})
    throw(ArgumentError("`project` cannot be applied to a `$(typeof(p))`; " *
                        "use a `SetAtom` instead"))
end

function project(n::Negation, vars::AbstractVector{Int})
    return Negation(project(n.p, vars))
end

function project(c::Conjunction, vars::AbstractVector{Int})
    conjuncts = similar(c.conjuncts)
    @inbounds for (i, (conjunct, n_args)) in enumerate(c.conjuncts)
        conjuncts[i] = (project(conjunct, vars), n_args)
    end
    return Conjunction(conjuncts)
end

function project(d::Disjunction, vars::AbstractVector{Int})
    disjuncts = similar(d.disjuncts)
    @inbounds for (i, (disjunct, n_args)) in enumerate(d.disjuncts)
        disjuncts[i] = (project(disjunct, vars), n_args)
    end
    return Disjunction(disjuncts)
end

# ==========================================
# Convenience functions to create predicates
# ==========================================

function contains(X::LazySet)
    return SetAtom(X, ⊆)
end

function is_contained_in(X::LazySet)
    return SetAtom(X, (X, Y) -> Y ⊆ X)
end

function is_disjoint_from(X::LazySet)
    return SetAtom(X, isdisjoint)
end

function intersects(X::LazySet)
    return Negation(is_disjoint_from(X))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ bff7f1e0-212b-11eb-0149-ed5bc1197f8b
# I used the Images Package.

using Images

# ╔═╡ d7b16030-2123-11eb-0cea-bf0b04b9f172
md"# Hello!
### This is my third Project Euler question! In this question, I will be trying to find the largest prime factor of the number: 600,851,475,143 ! "

# ╔═╡ c945b7d2-2128-11eb-3bc5-6f1da4fe5115
number = 600851475143

# ╔═╡ c4b7cb40-2128-11eb-2a57-333c8475e6d8
md"## What even is a prime number?"

# ╔═╡ f6687700-2125-11eb-3dea-cd257bcf1809
md"### A prime number is a number with only two factors - one and itself."

# ╔═╡ 13435e50-2129-11eb-3c21-9b3876d9c8f0
md"##### List of Prime numbers: 2, 3, 5, 7, 11 ..."

# ╔═╡ 35f18ad0-2129-11eb-33d3-e599a59910ce
md"### I cannot just **brute-force** and divide 600,851,475,143 by every number between 1 and the number."

# ╔═╡ c1d11420-212a-11eb-26da-ed0f2336bff4
md"## This is why I am using the Sieve of Eratosthenes to calculate a list of prime numbers up to 600,851,475,143 ."

# ╔═╡ 4523b9b0-2133-11eb-11fd-855e41b99d21
md"### This is the _Sieve of Eratosthenes_. It is a sieve that filters out the composite numbers and only primes remain."

# ╔═╡ 382d31e0-2148-11eb-0783-15e4e05d35c4
function Sieve_Eratosthenes(x)
	y = 1:x
	y.^1
end

# ╔═╡ aec527f0-283b-11eb-1d14-13bce147a4c7
Sieve_Eratosthenes(100)

# ╔═╡ 1f49e600-2855-11eb-1cb5-537b9494e376
countries = ["America", "UK", "Hong Kong", "Canada", "Japan", "Taiwan", "India", "Kenya"]

# ╔═╡ c081a3c0-290c-11eb-2182-8de89de4b8f4
countries

# ╔═╡ e8f09782-290c-11eb-0013-d715def4cf7f
countries[6] = "Taiwan"

# ╔═╡ 9dd279fe-2855-11eb-37b7-3b20511c9e03
Tuple = ("America", "UK", "Hong Kong", "Canada", "Japan", "Taiwan", "India", "Kenya")

# ╔═╡ b54aee10-2855-11eb-052c-db94049259c0
Tuple[1]

# ╔═╡ f50fb2b0-2855-11eb-1512-9faddf7ed772
countries[1] = "Mexico"

# ╔═╡ 146da2c0-2856-11eb-0662-29fa86a13c86
countries

# ╔═╡ 65894c00-2855-11eb-08ab-1b755fd41fa3
countries[8]

# ╔═╡ 2d8f58ce-283c-11eb-02ee-ad6b2242de59
[1, 2, 3]

# ╔═╡ 4e7a4750-2132-11eb-3c69-97021bf77041
md"# Appendix"

# ╔═╡ 9cc9a370-212c-11eb-3e55-8b1f4ec8ce2d
# I found a good diagram on the internet.

Image = load(download("https://www.cut-the-knot.org/Curriculum/Arithmetic/Eratosthenes.jpg"))

# ╔═╡ 15fb6022-212e-11eb-08f5-37d4443e81a1
size(Image)

# ╔═╡ 0f318d50-212e-11eb-284d-1522f5c7810c
# I cropped the original image.

Sieve_of_Eratosthenes = Image[35:502, 30:511]

# ╔═╡ f08b7b90-212d-11eb-3d4f-5120c5327fe8
Sieve_of_Eratosthenes

# ╔═╡ Cell order:
# ╟─d7b16030-2123-11eb-0cea-bf0b04b9f172
# ╟─c945b7d2-2128-11eb-3bc5-6f1da4fe5115
# ╟─c4b7cb40-2128-11eb-2a57-333c8475e6d8
# ╟─f6687700-2125-11eb-3dea-cd257bcf1809
# ╟─13435e50-2129-11eb-3c21-9b3876d9c8f0
# ╟─35f18ad0-2129-11eb-33d3-e599a59910ce
# ╟─c1d11420-212a-11eb-26da-ed0f2336bff4
# ╠═f08b7b90-212d-11eb-3d4f-5120c5327fe8
# ╟─4523b9b0-2133-11eb-11fd-855e41b99d21
# ╠═382d31e0-2148-11eb-0783-15e4e05d35c4
# ╠═aec527f0-283b-11eb-1d14-13bce147a4c7
# ╠═1f49e600-2855-11eb-1cb5-537b9494e376
# ╠═c081a3c0-290c-11eb-2182-8de89de4b8f4
# ╠═e8f09782-290c-11eb-0013-d715def4cf7f
# ╠═9dd279fe-2855-11eb-37b7-3b20511c9e03
# ╠═b54aee10-2855-11eb-052c-db94049259c0
# ╠═f50fb2b0-2855-11eb-1512-9faddf7ed772
# ╠═146da2c0-2856-11eb-0662-29fa86a13c86
# ╠═65894c00-2855-11eb-08ab-1b755fd41fa3
# ╠═2d8f58ce-283c-11eb-02ee-ad6b2242de59
# ╟─4e7a4750-2132-11eb-3c69-97021bf77041
# ╠═bff7f1e0-212b-11eb-0149-ed5bc1197f8b
# ╠═9cc9a370-212c-11eb-3e55-8b1f4ec8ce2d
# ╠═15fb6022-212e-11eb-08f5-37d4443e81a1
# ╠═0f318d50-212e-11eb-284d-1522f5c7810c


text:
import .GridWorlds

include("JuliaRL_BasicDQN_EmptyRoom.jl")


text:
using Test
using ITensorTDVP

test_path = joinpath(pkgdir(ITensorTDVP), "test")
test_files = filter(
  file -> startswith(file, "test_") && endswith(file, ".jl"), readdir(test_path)
)
@testset "ITensorTDVP.jl" begin
  @testset "$filename" for filename in test_files
    println("Running $filename")
    include(joinpath(test_path, filename))
  end
end

nothing