text:
# -*- coding: utf-8 -*-
# Load necessary libraries

# +
using Gen
using Statistics
using ProfileView #Make sure this is loaded before PyPlot to avoid stupid gtk version issues
using Profile
using PyPlot

include("DistributionsBacked.jl")
include("IterDeep.jl")
include("AnimatedPyplot.jl")
include("nuts.jl")

include("unrealizedNormals.jl")

const my_inv_gamma = DistributionsBacked{Float64}((alpha, theta) -> Distributions.InverseGamma(alpha, theta), [true, true], true)
# -
# Space for declaring constants.

disable_sample_logging()
INV_GAMMA_PRIOR_CONSTANT = 0.1;

# Canonical Rats data, from from section 6 of Gelfand et al (1990). See http://www.openbugs.net/Examples/Rats.html

ys_raw = ([151, 145, 147, 155, 135, 159, 141, 159, 177, 134,
    160, 143, 154, 171, 163, 160, 142, 156, 157, 152, 154, 139, 146,
    157, 132, 160, 169, 157, 137, 153, 199, 199, 214, 200, 188, 210,
    189, 201, 236, 182, 208, 188, 200, 221, 216, 207, 187, 203, 212,
    203, 205, 190, 191, 211, 185, 207, 216, 205, 180, 200, 246, 249,
    263, 237, 230, 252, 231, 248, 285, 220, 261, 220, 244, 270, 242,
    248, 234, 243, 259, 246, 253, 225, 229, 250, 237, 257, 261, 248,
    219, 244, 283, 293, 312, 272, 280, 298, 275, 297, 350, 260, 313,
    273, 289, 326, 281, 288, 280, 283, 307, 286, 298, 267, 272, 285,
    286, 303, 295, 289, 258, 286, 320, 354, 328, 297, 323, 331, 305,
    338, 376, 296, 352, 314, 325, 358, 312, 324, 316, 317, 336, 321,
    334, 302, 302, 323, 331, 345, 333, 316, 291, 324])
ys = reshape([Float64(y) for y = ys_raw],30,5)
xs = [8.0, 15.0, 22.0, 29.0, 36.0]
ys = ys ./ 100.
(N,T) = size(ys)

# These data are the weights of 30 rats, measured at 5 common age values in days. The model is a simple hierarchical/random effects one; a random intercept and random slope for each rat, plus an error term (presumably representing not measurement error, but just the random variation of individual growth curves).

for i in 1:30
    plot(xs,ys[i,:])
end
xlabel("Age (days)")
ylabel("Weight (g)")
title("Rat growth")

# +
@gen function oneRat(xs, T, mu_alpha, sigmasq_alpha, mu_beta, sigmasq_beta, sigmasq_y)

    alpha ~ normal(mu_alpha, sqrt(sigmasq_alpha))
    beta ~ normal(mu_beta, sqrt(sigmasq_beta))
    y ~ broadcasted_normal([alpha + beta * (xs[t]) for t = 1:T],
                                        sqrt(sigmasq_y))
end

@gen function ratsModel(xs::Vector,N::Int32,)
    T = length(xs)
    xbar = mean(xs) #could be precomputed, but YKWTS about premature optimization...

    mu_alpha ~ normal(0, 10.0)
    mu_beta ~ normal(0, 1.0)
    sigmasq_y ~ my_inv_gamma(INV_GAMMA_PRIOR_CONSTANT, INV_GAMMA_PRIOR_CONSTANT)
    sigmasq_alpha ~ my_inv_gamma(INV_GAMMA_PRIOR_CONSTANT, INV_GAMMA_PRIOR_CONSTANT)
    sigmasq_beta ~ my_inv_gamma(INV_GAMMA_PRIOR_CONSTANT, INV_GAMMA_PRIOR_CONSTANT)

    alpha = Vector{Float64}(undef,N)
    beta = Vector{Float64}(undef,N)
    y = Vector{Vector{Float64}}(undef,N)
    for n in 1:N
        {:data => n} ~ oneRat(xs .- xbar, T, mu_alpha, sigmasq_alpha, mu_beta, sigmasq_beta, sigmasq_y)
    end
end

model_transformations = choicemap()
model_transformations[:sigmasq_y] = transform_log
model_transformations[:sigmasq_alpha] = transform_log
model_transformations[:sigmasq_beta] = transform_log

;
# -

function make_constraints(ys::Array)
    (N,T) = size(ys)
    constraints = Gen.choicemap()
    for n in 1:N
        constraints[:data => n => :y] = ys[n,:]
    end
    constraints
end
;

function mcmc_inference(xs, ys, num_iters, update, model = ratsModel)
    results = Vector{Gen.DynamicDSLTrace{DynamicDSLFunction{Any}}}(undef, num_iters+1)
    N = size(ys)[1]
    observations = make_constraints(ys)
    (trace, _) = generate(model, (xs, N), observations)
    results[1] = trace
    selection = Gen.complement(select(observations))
    for iter=1:num_iters
        trace = update(trace,N,selection)
        results[iter+1] = trace
    end
    results
end
;

# +

@gen function oneDeferredRat(xs, T, mu_alpha, sigmasq_alpha, mu_beta, sigmasq_beta, sigmasq_y)

    alpha ~ unrealizedIndyNormal(mu_alpha, sqrt(sigmasq_alpha))
    beta ~ unrealizedIndyNormal(mu_beta, sqrt(sigmasq_beta))
    for t = 1:T
        @trace(observedNormal(alpha + beta * (xs[t]), sqrt(sigmasq_y)), "y$(t)")
    end
end

@gen function deferredRatsModel(xs::Vector,N::Int32,)
    T = length(xs)
    xbar = mean(xs) #could be precomputed, but YKWTS about premature optimization...

    mu_alpha ~ normal(0, 10.0)
    mu_beta ~ normal(0, 1.0)
    sigmasq_y ~ my_inv_gamma(INV_GAMMA_PRIOR_CONSTANT, INV_GAMMA_PRIOR_CONSTANT)
    sigmasq_alpha ~ my_inv_gamma(INV_GAMMA_PRIOR_CONSTANT, INV_GAMMA_PRIOR_CONSTANT)
    sigmasq_beta ~ my_inv_gamma(INV_GAMMA_PRIOR_CONSTANT, INV_GAMMA_PRIOR_CONSTANT)

    alpha = Vector{Float64}(undef,N)
    beta = Vector{Float64}(undef,N)
    y = Vector{Vector{Float64}}(undef,N)
    for n in 1:N
        {:data => n} ~ oneDeferredRat(xs .- xbar, T, mu_alpha, sigmasq_alpha, mu_beta, sigmasq_beta, sigmasq_y)
    end
end

model_transformations = choicemap()
model_transformations[:sigmasq_y] = transform_log
model_transformations[:sigmasq_alpha] = transform_log
model_transformations[:sigmasq_beta] = transform_log

;

# +
Profile.clear()

function block_mh_d(tr,N,selection)
    (tr, _) = mh(tr, select(:mu_alpha, :mu_beta))
    (tr, _) = mh(tr, select(:sigmasq_alpha, :sigmasq_beta, :sigmasq_y))

    for n in 1:N
        (tr, _) = mh(tr, select(:data => n => :y1,:data => n => :y2,:data => n => :y3,:data => n => :y4,
                                :data => n => :alpha,:data => n => :beta))
    end
    tr
end
;

VIZ_INTERVAL = 1
VIZ_FRAMES = 20
trs3= Profile.@profile mcmc_inference(xs, ys, VIZ_INTERVAL * VIZ_FRAMES, block_mh_d, deferredRatsModel)
tr3 = trs3[end]
print("mu_alpha: $(tr3[:mu_alpha]), mu_beta: $(tr3[:mu_beta]), sigma_y: $(sqrt(tr3[:sigmasq_y]))\n")
;

# +
c = make_constraints(ys)
alphahats = mean(ys,dims=2)
alphahat = mean(alphahats)
betahats = transpose(ys)\xs
betahat = mean(betahats)
c[:mu_alpha] = alphahat
c[:mu_beta] = betahat

c[:sigmasq_beta] = var(betahats) * 2.1601
c[:sigmasq_alpha] = var(alphahats) * 2.1601
c[:sigmasq_y] = 0.000001

tr, w = generate(deferredRatsModel, (ys, N), c) #w == -N/2

s = get_score(tr)

c[:sigmasq_y] = 1000.

c[:sigmasq_beta] = var(betahats) * 2.1603
c[:sigmasq_alpha] = var(alphahats) * 2.1603
c[:sigmasq_y] = 0.000001

tr2, w2 = generate(deferredRatsModel, (ys, N), c) #w == -N/2

s2 = get_score(tr2)

w, s, w2, s2
# -

#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
# ...
#
#
#
#
#
#
#
#

function block_mh(tr,N,selection)
    (tr, _) = mh(tr, select(:mu_alpha, :mu_beta))
    (tr, _) = mh(tr, select(:sigmasq_alpha, :sigmasq_beta, :sigmasq_y))

    for n in 1:N
        (tr, _) = mh(tr, select(:data => n => :y,
                                :data => n => :alpha,:data => n => :beta))
    end
    tr
end
;

function simple_hmc(tr,N,selection)
    (tr, _) = hmc(tr, selection)
    tr
end
;

VIZ_INTERVAL = 2
VIZ_FRAMES = 10
trs= mcmc_inference(xs, ys, VIZ_INTERVAL * VIZ_FRAMES, block_mh)
tr = trs[end]
get_choices(tr)
;

print("mu_alpha: $(tr[:mu_alpha]), mu_beta: $(tr[:mu_beta]), sigma_y: $(sqrt(tr[:sigmasq_y]))")

# +
fewcolors = ["r", "b", "g"] #In order to avoid busy graphs, I graph just 3 rats per plot
numPlots = 2 #...but tile 2 plots
numRats = length(fewcolors)
CONF95SIGS = 1.96

function make_axes()
    fig = figure("A few rats",figsize=(5,5))
    my_axes = fig.subplots(nrows=numPlots)
    for ax in my_axes
        ax.set_xlim((5,40))
    end
    (fig,my_axes)
end

function visualize_rats(trace,my_axes)
    xbar = mean(xs)
    xgrid = collect(min(xs...):.5:max(xs...))
    xlims = [min(xs...), max(xs...)]
    centeredgrid = xgrid .- xbar
    linesigmas = sqrt.(trace[:sigmasq_alpha] .+ trace[:sigmasq_beta] .* centeredgrid .^ 2)
    meanline = trace[:mu_alpha] .+ trace[:mu_beta] .* centeredgrid
    upperline = meanline .+ CONF95SIGS .* linesigmas
    lowerline = meanline .- CONF95SIGS .* linesigmas
    for (j,ax) in enumerate(my_axes)
        ax.clear()
        ax.fill_between(xgrid,lowerline,upperline,alpha=.1,color="black")
        ax.fill_between([6.,7],
                        fill(250. + CONF95SIGS * sqrt(trace[:sigmasq_y]),2),
                        fill(250. - CONF95SIGS * sqrt(trace[:sigmasq_y]),2),
            alpha=1,color="black")

        ax.set_ylim((min(1.00,lowerline[1],lowerline[end]),
                max(4.00,upperline[1],upperline[end])))
        for k in 1:numRats
            i = j*numRats + k
            ax.plot(xs,ys[i,:],color=fewcolors[k], alpha=1.)

            ratline = trace[:data => i => :alpha] .+ trace[:data => i => :beta] .* (xlims .- xbar)
            ax.plot(xlims, ratline,color=fewcolors[k], alpha=.3)
        end
        ax.set_ylabel("Weight (g)")
        if j == 1
            ax.set_title("Rat growth")
        else
            ax.set_xlabel("Age (days)")
        end

    end
end

# +

visualize_rats(tr,make_axes()[2])
# -

display_animation("rats_mh", make_axes,
    (ax)->visualize_rats(trs[1],ax),
    (ax,i)->visualize_rats(trs[1+i*VIZ_INTERVAL],ax),
    VIZ_FRAMES)

VIZ_INTERVAL = 2
VIZ_FRAMES = 10
trs2= mcmc_inference(xs, ys, VIZ_INTERVAL * VIZ_FRAMES, simple_hmc)
tr2 = trs2[end]
get_choices(tr2)
;

print("mu_alpha: $(tr2[:mu_alpha]), mu_beta: $(tr2[:mu_beta]), sigma_y: $(sqrt(tr2[:sigmasq_y]))\n")

if false
    display_animation("rats_hmc", make_axes,
        (ax)->visualize_rats(trs[1],ax),
        (ax,i)->visualize_rats(trs[1+i*VIZ_INTERVAL],ax),
        VIZ_FRAMES)
end

# +
Profile.clear()

VIZ_INTERVAL = 1
VIZ_FRAMES = 20
trs3= Profile.@profile mcmc_inference(xs, ys, VIZ_INTERVAL * VIZ_FRAMES, (tr, N, selection)->my_nuts(tr, selection, model_transformations))
tr3 = trs3[end]
print("mu_alpha: $(tr3[:mu_alpha]), mu_beta: $(tr3[:mu_beta]), sigma_y: $(sqrt(tr3[:sigmasq_y]))\n")
;
# -

ProfileView.view()

open("profile.txt", "w") do s
    Profile.print(IOContext(s, :displaysize => (24, 500)))
end

display_animation("rats_nuts", make_axes,
    (ax)->visualize_rats(trs3[1],ax),
    (ax,i)->visualize_rats(trs3[1+i*VIZ_INTERVAL],ax),
    VIZ_FRAMES,200)

# The below is just a sketch of how I think "semi-automatic conjugacy" might work. It's not working code, yet.

l = [(4,6)]
for (a, (b,c)) in enumerate(l)
    @warn "hi" a b c
end

max_stars_repo_name: jamesonquinn/Gen

max_issues_repo_name: jamesonquinn/Gen

max_forks_repo_name: jamesonquinn/Gen


text:
import ForwardDiff  # Needs to be first due to https://github.com/JuliaMath/Interpolations.jl/issues/207
using AdaptiveDistanceFields
using Base.Test
using StaticArrays: SVector

@testset "coarse" begin
    rtol = 0.05
    atol = 0.1
    s = x -> sqrt(sum((x - SVector(1, 2)).^2))
    adf = AdaptiveDistanceField(s, SVector(-1., -1), SVector(5., 5), rtol, atol)

    for x in linspace(-1, 4)
        for y in linspace(-1, 4)
            @test isapprox(adf(SVector(x, y)), s(SVector(x, y)), atol=atol, rtol=rtol)
        end
    end
end

@testset "fine" begin
    rtol = 0.001
    atol = 0.001
    s = x -> sqrt(sum((x - SVector(-2, -1)).^2))
    adf = AdaptiveDistanceField(s, SVector(-5., -5), SVector(10., 10), rtol, atol)

    for x in linspace(-5, 5)
        for y in linspace(-5, 5)
            @test isapprox(adf(SVector(x, y)), s(SVector(x, y)), atol=0.005, rtol=0.005)
        end
    end
end

@testset "extrapolation" begin
    s = x -> sqrt(sum((x - SVector(-2, -1)).^2))
    adf = AdaptiveDistanceField(s, SVector(-5., -5), SVector(10., 10))

    @test isapprox(ForwardDiff.derivative(x -> adf(SVector(x, -1)), 100), 1, atol=1e-2, rtol=1e-2)
    @test isapprox(ForwardDiff.derivative(x -> adf(SVector(x, -1)), -100), -1, atol=1e-2, rtol=1e-2)
    @test isapprox(ForwardDiff.derivative(y -> adf(SVector(-2, y)), 100), 1, atol=1e-2, rtol=1e-2)
    @test isapprox(ForwardDiff.derivative(y -> adf(SVector(-2, y)), -100), -1, atol=1e-2, rtol=1e-2)
end

ext: jl


text:
using Test
using MIPVerify
using MIPVerify: LInfNormBoundedPerturbationFamily
@isdefined(TestHelpers) || include("../TestHelpers.jl")

@timed_testset "integration.jl" begin
    mnist = read_datasets("MNIST")
    nn_wk17a = get_example_network_params("MNIST.WK17a_linf0.1_authors")
    @testset "batch_find_untargeted_attack" begin
        mktempdir() do dir
            MIPVerify.batch_find_untargeted_attack(
                nn_wk17a,
                mnist.test,
                [1], # robust sample
                TestHelpers.get_main_solver(),
                solve_rerun_option=MIPVerify.never,
                pp=MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                norm_order=Inf,
                rebuild=true,
                tightening_algorithm=lp,
                tightening_solver=TestHelpers.get_tightening_solver(),
                cache_model=false,
                solve_if_predicted_in_targeted=false,
                save_path=dir
            )
        end

        mktempdir() do dir
            MIPVerify.batch_find_untargeted_attack(
                nn_wk17a,
                mnist.test,
                [9, 248], # non-robust and misclassified sample
                TestHelpers.get_main_solver(),
                solve_rerun_option=MIPVerify.never,
                pp=MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                norm_order=Inf,
                rebuild=true,
                tightening_algorithm=interval_arithmetic,
                cache_model=false,
                solve_if_predicted_in_targeted=false,
                save_path=dir
            )
        end
    end

    @testset "batch_find_targeted_attack" begin
        mktempdir() do dir
            MIPVerify.batch_find_targeted_attack(
                nn_wk17a,
                mnist.test,
                [1],
                TestHelpers.get_main_solver(),
                solve_rerun_option=MIPVerify.never,
                pp=MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                norm_order=Inf,
                tightening_algorithm=interval_arithmetic,
                tightening_solver=TestHelpers.get_tightening_solver(),
                cache_model=false,
                solve_if_predicted_in_targeted=false,
                target_labels=[1, 8],
                save_path=dir
            )
        end
    end

    @testset "batch_build_model" begin
        mktempdir() do dir
            MIPVerify.batch_build_model(
                nn_wk17a,
                mnist.test,
                [1],
                TestHelpers.get_tightening_solver(),
                pp=MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                tightening_algorithm=interval_arithmetic
            )
        end
    end
end


text:
using JuliaDot
using Test
using PyCall

@testset "LoadPydotPlus" begin
    JuliaDot.PydotPlus != PyCall.PyPtr_NULL
end

@testset "draw graph" begin
    dot_data = """
digraph G {

  subgraph cluster_0 {
    style=filled;
    color=lightgrey;
    node [style=filled,color=white];
    a0 -> a1 -> a2 -> a3;
    label = "process #1";
  }

  subgraph cluster_1 {
    node [style=filled];
    b0 -> b1 -> b2 -> b3;
    label = "process #2";
    color=blue
  }
  start -> a0;
  start -> b0;
  a1 -> b3;
  b2 -> a3;
  a3 -> a0;
  a3 -> end;
  b3 -> end;

  start [shape=Mdiamond];
  end [shape=Msquare];
}
"""
    println(dot_data)
end

@testset "draw graph" begin
    dot_data = """
digraph G {

  subgraph cluster_0 {
    style=filled;
    color=lightgrey;
    node [style=filled,color=white];
    a0 -> a1 -> a2 -> a3;
    label = "process #1";
  }

  subgraph cluster_1 {
    node [style=filled];
    b0 -> b1 -> b2 -> b3;
    label = "process #2";
    color=blue
  }
  start -> a0;
  start -> b0;
  a1 -> b3;
  b2 -> a3;
  a3 -> a0;
  a3 -> end;
  b3 -> end;

  start [shape=Mdiamond];
  end [shape=Msquare];
}
"""
    println(length(dot(dot_data)))
end


text:
#### Operators for a system with a body and stationary points ####

export bc_constraint_rhs!

function bc_constraint_rhs!(us::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    # Here, need surface_velocity! - _velocity_single_layer!(u,sys,t) - _velocity_freestream!(u,sys,t)
    sys.Vv .= 0.0
    _velocity_single_layer!(sys.Vv,sys,t)
    _velocity_freestream!(sys.Vv,sys,t)
    us .= sys.Ef*sys.Vv
    us .*= -1.0
    surface_velocity!(sys.Δus,sys,t)
    us .+= sys.Δus
    return us
end

# Constraint operators, using stored regularization and interpolation operators
# B₁ᵀ = CᵀEᵀ
function ns_op_constraint_force!(out::Nodes{Dual,NX,NY},τ::VectorData{N},sys::NavierStokes{NX,NY,N}) where {NX,NY,N}
    sys.Vf .= sys.Rf*τ
    out .= 0.0
    curl!(out,sys.Vf)
    return out
end

function _vel_ns_op_constraint_force!(u::Edges{Primal,NX,NY},τ,sys::NavierStokes{NX,NY,0}) where {NX,NY}
    return u
end

function _vel_ns_op_constraint_force!(u::Edges{Primal,NX,NY},τ::VectorData{N},sys::NavierStokes{NX,NY,N}) where {NX,NY,N}
    u .= sys.Rf*τ
    return u
end

# B₂ = -ECL⁻¹
function bc_constraint_op!(out::VectorData{N},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N}) where {NX,NY,N}
    sys.Vv .= 0.0
    ViscousFlow._velocity_vorticity!(sys.Vv,w,sys)
    out .= sys.Ef*sys.Vv
    return out
end


text:
using Test
using CSV
include("../src/TableProcessing.jl")
#include("src/TableProcessing.jl")

#Check that reading in test data works
df = CSV.read("test/data/test_counts_matrix.txt",
header = true, delim = '\t')
@test ReadCountsMatrix("test/data/test_counts_matrix.txt") == df

#Check that an error is thrown when no filename is passed
@test_throws MethodError ReadCountsMatrix()

#Same tests for gene-isoform table
genes = CSV.read("test/data/all_isoforms.txt", header = false, delim = ' ')
@test ReadGeneIsoformRelationships("test/data/all_isoforms.txt") == genes
@test_throws MethodError ReadGeneIsoformRelationships()

#Check function returns true when genes match isoforms
df = DataFrame(Column1 = ["I1", "I2", "I3", "I4"], Cell1 = [0, 0, 0, 0])
genes = DataFrame(Column1 = ["G1", "G1", "G2", "G3"],
Column2 = ["I2", "I3", "I1", "I4"])
@test CheckGenesMatchIsoforms(df, genes) == true

#Check function returns false when genes don't match isoforms
@test CheckGenesMatchIsoforms(df, genes[1:2, :]) == false

#Check we get expected output when adding gene information
df_test = DataFrame(Genes = ["G1", "G1", "G2", "G3"],
Isoforms = ["I2", "I3", "I1", "I4"], Cell1 = [0, 0, 0, 0])
@test AddGeneInfo(df, genes) == df_test

#Check that we correctly convert from dataframe to dict
df_test = DataFrame(Genes = ["G1", "G1", "G2", "G3"],
Isoforms = ["I2", "I3", "I1", "I4"], Cell1 = [5, 5, 0, 5], Cell2 = [5,5,5,5])

dict_test_two = Dict("G1" => Dict("I2" => [5,5], "I3" => [5,5]))
dict_test_one = Dict("G3" => Dict("I4" => [5,5]))

#Check that we keep expected genes when N = 1 to N = 4
@test KeepGenesWithNIsoforms(df_test, 1) == dict_test_one
@test KeepGenesWithNIsoforms(df_test, 2) == dict_test_two

#Check that an error is thrown when N<1
@test_throws ArgumentError KeepGenesWithNIsoforms(df_test, 0)
@test_throws ArgumentError KeepGenesWithNIsoforms(df_test, -10)

#Check that it works when two of three isoforms are expressed
df_test = DataFrame(Genes = ["G1", "G1", "G1", "G3"],
Isoforms = ["I2", "I3", "I1", "I4"], Cell1 = [5, 5, 0, 5], Cell2 = [5,5,5,5])
@test KeepGenesWithNIsoforms(df_test, 2) == dict_test_two

df_test = DataFrame(Genes = ["G1", "G1", "G1"],
Isoforms = ["I2", "I3", "I1"], Cell1 = [5, 5, 0], Cell2 = [5,5,5])
dict_test_two = Dict("G1" => Dict("I2" => [5,5], "I3" => [5,5]))
@test ConvertExpressedIsosToDict(df_test) == dict_test_two

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
function darcy_velocity(h, Ks, mins, maxs, ns)#note: this is currently not very Zygote-compatible
	allpoints = map(getpoints, mins, maxs, ns)
	h_itp_unscaled = Interpolations.interpolate(h, Interpolations.BSpline(Interpolations.Quadratic(Interpolations.Line(Interpolations.OnCell()))))
	h_itp = Interpolations.scale(h_itp_unscaled, reverse(allpoints)...)
	K_itp_unscaled = Interpolations.interpolate(Ks, Interpolations.BSpline(Interpolations.Quadratic(Interpolations.Line(Interpolations.OnCell()))))
	K_itp = Interpolations.scale(K_itp_unscaled, reverse(allpoints)...)
	return (x...)->reverse(K_itp(x...) * Interpolations.gradient(h_itp, x...))
end

function getpoints(min, max, n)
	return range(min, max; length=n)
end

function load_uge(filename)
	fuge = open(filename)
	line = readline(fuge)
	numcells= parse(Int, split(line)[2])
	volumes = Array{Float64}(undef, numcells)
	coords = Array{Float64}(undef, 3, numcells)
	for i in 1:numcells
		line = readline(fuge)
		lineparts = split(line)
		coords[1, i] = parse(Float64, lineparts[2])
		coords[2, i] = parse(Float64, lineparts[3])
		coords[3, i] = parse(Float64, lineparts[4])
		volumes[i] = parse(Float64, lineparts[5])
	end
	line = readline(fuge)
	numconnections = parse(Int, split(line)[2])
	neighbors = Array{Pair{Int, Int}}(undef, numconnections)
	areas = Array{Float64}(undef, numconnections)
	for i in 1:numconnections
		line = readline(fuge)
		lineparts = split(line)
		neighbors[i] = parse(Int, lineparts[1])=>parse(Int, lineparts[2])
		areas[i] = parse(Float64, lineparts[6])
	end
	close(fuge)
	lengths = map(n->sqrt(sum((coords[j, n[1]] - coords[j, n[2]]) .^ 2 for j = 1:3)), neighbors)
	return coords, volumes, neighbors, areas, lengths
end

function regulargrid1d(mins, maxs, ns, dy, dz)
	coords = Array{Float64}(undef, 1, prod(ns))
	xs = getpoints(mins[1], maxs[1], ns[1])
	dx = xs[2] - xs[1]
	j = 1
	neighbors = Array{Pair{Int, Int}}(undef, ns[1] - 1)
	areasoverlengths = Array{Float64}(undef, ns[1] - 1)
	volumes = fill(dx * dy * dz, ns[1])
	areasoverlengths = fill(dy * dz / dx, ns[1])
	for i1 = 1:ns[1]
			coords[1, i1] = xs[i1]
			if i1 < ns[1]
				neighbors[j] = i1=>i1 + 1
				j += 1
			end
	end
	return coords, neighbors, areasoverlengths, volumes
end

function regulargrid2d(mins, maxs, ns, dz)
	linearindex = (i1, i2)->i2 + ns[2] * (i1 - 1)
	coords = Array{Float64}(undef, 2, prod(ns))
	xs = getpoints(mins[1], maxs[1], ns[1])
	ys = getpoints(mins[2], maxs[2], ns[2])
	dx = xs[2] - xs[1]
	dy = ys[2] - ys[1]
	j = 1
	neighbors = Array{Pair{Int, Int}}(undef, 2 * prod(ns) - ns[1] - ns[2])
	areasoverlengths = Array{Float64}(undef, 2 * prod(ns) - ns[1] - ns[2])
	volumes = Array{Float64}(undef, 0)
	for i1 = 1:ns[1]
		for i2 = 1:ns[2]
			push!(volumes, dx * dy * dz)
			coords[1, linearindex(i1, i2)] = xs[i1]
			coords[2, linearindex(i1, i2)] = ys[i2]
			if i1 < ns[1]
				neighbors[j] = linearindex(i1, i2)=>linearindex(i1 + 1, i2)
				areasoverlengths[j] = dy * dz / dx
				j += 1
			end
			if i2 < ns[2]
				neighbors[j] = linearindex(i1, i2)=>linearindex(i1, i2 + 1)
				areasoverlengths[j] = dx * dz / dy
				j += 1
			end
		end
	end
	return coords, neighbors, areasoverlengths, volumes
end

function regulargrid3d(mins, maxs, ns)
	linearindex = (i1, i2, i3)->i3 + ns[3] * (i2 - 1) + ns[3] * ns[2] * (i1 - 1)
	coords = Array{Float64}(undef, 3, prod(ns))
	xs = getpoints(mins[1], maxs[1], ns[1])
	ys = getpoints(mins[2], maxs[2], ns[2])
	zs = getpoints(mins[3], maxs[3], ns[3])
	dx = xs[2] - xs[1]
	dy = ys[2] - ys[1]
	dz = zs[2] - zs[1]
	j = 1
	neighbors = Array{Pair{Int, Int}}(undef, 3 * prod(ns) - ns[1] * ns[2] - ns[1] * ns[3] - ns[2] * ns[3])
	areasoverlengths = Array{Float64}(undef, 3 * prod(ns) - ns[1] * ns[2] - ns[1] * ns[3] - ns[2] * ns[3])
	volumes = fill(dx * dy * dz, size(coords, 2))
	for i1 = 1:ns[1]
		for i2 = 1:ns[2]
			for i3 = 1:ns[3]
				coords[1, linearindex(i1, i2, i3)] = xs[i1]
				coords[2, linearindex(i1, i2, i3)] = ys[i2]
				coords[3, linearindex(i1, i2, i3)] = zs[i3]
				if i1 < ns[1]
					neighbors[j] = linearindex(i1, i2, i3)=>linearindex(i1 + 1, i2, i3)
					areasoverlengths[j] = dy * dz / dx
					j += 1
				end
				if i2 < ns[2]
					neighbors[j] = linearindex(i1, i2, i3)=>linearindex(i1, i2 + 1, i3)
					areasoverlengths[j] = dx * dz / dy
					j += 1
				end
				if i3 < ns[3]
					neighbors[j] = linearindex(i1, i2, i3)=>linearindex(i1, i2, i3 + 1)
					areasoverlengths[j] = dx * dy / dz
					j += 1
				end
			end
		end
	end
	return coords, neighbors, areasoverlengths, volumes
end


text:
struct Fermion{P,I<:Sector} <: Sector
    sector::I
    function Fermion{P,I}(sector::I) where {P, I<:Sector}
        @assert BraidingStyle(I) isa Bosonic
        return new{P,I}(sector)
    end
end
Fermion{P}(sector::I) where {P, I<:Sector} = Fermion{P,I}(sector)
Fermion{P,I}(sector) where {P, I<:Sector} = Fermion{P,I}(convert(I, sector))
Base.convert(::Type{Fermion{P,I}}, a::Fermion{P,I}) where {P, I<:Sector} = a
Base.convert(::Type{Fermion{P,I}}, a) where {P, I<:Sector} = Fermion{P,I}(convert(I, a))

fermionparity(f::Fermion{P}) where P = P(f.sector)

Base.IteratorSize(::Type{SectorValues{Fermion{P,I}}}) where {P, I<:Sector} =
    Base.IteratorSize(SectorValues{I})
Base.length(::SectorValues{Fermion{P,I}}) where {P, I<:Sector} = length(values(I))
function Base.iterate(::SectorValues{Fermion{P, I}}) where {P, I<:Sector}
    next = iterate(values(I))
    @assert next !== nothing
    value, state = next
    return Fermion{P}(value), state
end
function Base.iterate(::SectorValues{Fermion{P, I}}, state) where {P, I<:Sector}
    next = iterate(values(I), state)
    if next === nothing
        return nothing
    else
        value, state = next
        return Fermion{P}(value), state
    end
end
Base.getindex(::SectorValues{Fermion{P, I}}, i) where {P, I<:Sector} =
    Fermion{P}(values(I)[i])
findindex(::SectorValues{Fermion{P, I}}, f::Fermion{P, I}) where {P, I<:Sector} =
    findindex(values(I), f.sector)

Base.one(::Type{Fermion{P, I}}) where {P, I<:Sector} = Fermion{P}(one(I))
Base.conj(f::Fermion{P}) where {P} = Fermion{P}(conj(f.sector))

dim(f::Fermion) = dim(f.sector)

FusionStyle(::Type{<:Fermion{<:Any,I}}) where {I<:Sector} = FusionStyle(I)
BraidingStyle(::Type{<:Fermion}) = Fermionic()
Base.isreal(::Type{Fermion{<:Any,I}}) where {I<:Sector} = isreal(I)

⊗(a::F, b::F) where {F<:Fermion} = SectorSet{F}(a.sector ⊗ b.sector)

Nsymbol(a::F, b::F, c::F) where {F<:Fermion} = Nsymbol(a.sector, b.sector, c.sector)

Fsymbol(a::F, b::F, c::F, d::F, e::F, f::F) where {F<:Fermion} =
    Fsymbol(a.sector, b.sector, c.sector, d.sector, e.sector, f.sector)

function Rsymbol(a::F, b::F, c::F) where {F<:Fermion}
    if fermionparity(a) && fermionparity(b)
        return -Rsymbol(a.sector, b.sector, c.sector)
    else
        return +Rsymbol(a.sector, b.sector, c.sector)
    end
end

twist(a::Fermion) = ifelse(fermionparity(a), -1, +1)*twist(a.sector)

type_repr(::Type{Fermion{P,I}}) where {P, I<:Sector} = "Fermion{$P, " * type_repr(I) * "}"

function Base.show(io::IO, a::Fermion{P, I}) where {P, I<:Sector}
    if get(io, :typeinfo, nothing) !== Fermion{P, I}
        print(io, type_repr(typeof(a)), "(")
    end
    print(IOContext(io, :typeinfo => I), a.sector)
    if get(io, :typeinfo, nothing) !== Fermion{P, I}
        print(io, ")")
    end
end

Base.hash(f::Fermion, h::UInt) = hash(f.sector, h)
Base.isless(a::F, b::F) where {F<:Fermion} = isless(a.sector, b.sector)

_fermionparity(a::Z2Irrep) = isodd(a.n)
_fermionnumber(a::U1Irrep) = isodd(convert(Int, a.charge))
_fermionspin(a::SU2Irrep) = isodd(twice(a.j))

const FermionParity = Fermion{_fermionparity, Z2Irrep}
const FermionNumber = Fermion{_fermionnumber, U1Irrep}
const FermionSpin = Fermion{_fermionspin, SU2Irrep}
const fℤ₂ = FermionParity
const fU₁ = FermionNumber
const fSU₂ = FermionSpin

type_repr(::Type{FermionParity}) = "FermionParity"
type_repr(::Type{FermionNumber}) = "FermionNumber"
type_repr(::Type{FermionSpin}) = "FermionSpin"


text:
## removed in CUDA 11

@checked function cusparseScsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_context()
    ccall((:cusparseScsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{Cfloat}, CuPtr{Cint},
                    CuPtr{Cint}, CuPtr{Cfloat}, CuPtr{Cint}, CuPtr{Cint}, cusparseAction_t,
                    cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseDcsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_context()
    ccall((:cusparseDcsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{Cdouble}, CuPtr{Cint},
                    CuPtr{Cint}, CuPtr{Cdouble}, CuPtr{Cint}, CuPtr{Cint},
                    cusparseAction_t, cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseCcsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_context()
    ccall((:cusparseCcsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{cuComplex}, CuPtr{Cint},
                    CuPtr{Cint}, CuPtr{cuComplex}, CuPtr{Cint}, CuPtr{Cint},
                    cusparseAction_t, cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseZcsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_context()
    ccall((:cusparseZcsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{cuDoubleComplex},
                    CuPtr{Cint}, CuPtr{Cint}, CuPtr{cuDoubleComplex}, CuPtr{Cint},
                    CuPtr{Cint}, cusparseAction_t, cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseScsrmm2(handle, transA, transB, m, n, k, nnz, alpha, descrA,
                                  csrSortedValA, csrSortedRowPtrA, csrSortedColIndA, B,
                                  ldb, beta, C, ldc)
    initialize_context()
    ccall((:cusparseScsrmm2, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, cusparseOperation_t, cusparseOperation_t, Cint,
                    Cint, Cint, Cint, Ref{Cfloat}, cusparseMatDescr_t, CuPtr{Cfloat},
                    CuPtr{Cint}, CuPtr{Cint}, CuPtr{Cfloat}, Cint, Ref{Cfloat},
                    CuPtr{Cfloat}, Cint),
                   handle, transA, transB, m, n, k, nnz, alpha, descrA, csrSortedValA,
                   csrSortedRowPtrA, csrSortedColIndA, B, ldb, beta, C, ldc)
end

@checked function cusparseDcsrmm2(handle, transA, transB, m, n, k, nnz, alpha, descrA,
                                  csrSortedValA, csrSortedRowPtrA, csrSortedColIndA, B,
                                  ldb, beta, C, ldc)
    initialize_context()
    ccall((:cusparseDcsrmm2, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, cusparseOperation_t, cusparseOperation_t, Cint,
                    Cint, Cint, Cint, Ref{Cdouble}, cusparseMatDescr_t, CuPtr{Cdouble},
                    CuPtr{Cint}, CuPtr{Cint}, CuPtr{Cdouble}, Cint, Ref{Cdouble},
                    CuPtr{Cdouble}, Cint),
                   handle, transA, transB, m, n, k, nnz, alpha, descrA, csrSortedValA,
                   csrSortedRowPtrA, csrSortedColIndA, B, ldb, beta, C, ldc)
end

@checked function cusparseCcsrmm2(handle, transA, transB, m, n, k, nnz, alpha, descrA,
                                  csrSortedValA, csrSortedRowPtrA, csrSortedColIndA, B,
                                  ldb, beta, C, ldc)
    initialize_context()
    ccall((:cusparseCcsrmm2, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, cusparseOperation_t, cusparseOperation_t, Cint,
                    Cint, Cint, Cint, Ref{cuComplex}, cusparseMatDescr_t, CuPtr{cuComplex},
                    CuPtr{Cint}, CuPtr{Cint}, CuPtr{cuComplex}, Cint, Ref{cuComplex},
                    CuPtr{cuComplex}, Cint),
                   handle, transA, transB, m, n, k, nnz, alpha, descrA, csrSortedValA,
                   csrSortedRowPtrA, csrSortedColIndA, B, ldb, beta, C, ldc)
end

@checked function cusparseZcsrmm2(handle, transA, transB, m, n, k, nnz, alpha, descrA,
                                  csrSortedValA, csrSortedRowPtrA, csrSortedColIndA, B,
                                  ldb, beta, C, ldc)
    initialize_context()
    ccall((:cusparseZcsrmm2, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, cusparseOperation_t, cusparseOperation_t, Cint,
                    Cint, Cint, Cint, Ref{cuDoubleComplex}, cusparseMatDescr_t,
                    CuPtr{cuDoubleComplex}, CuPtr{Cint}, CuPtr{Cint},
                    CuPtr{cuDoubleComplex}, Cint, Ref{cuDoubleComplex},
                    CuPtr{cuDoubleComplex}, Cint),
                   handle, transA, transB, m, n, k, nnz, alpha, descrA, csrSortedValA,
                   csrSortedRowPtrA, csrSortedColIndA, B, ldb, beta, C, ldc)
end


text:
df        = CSV.File(joinpath(path, "csv", "Penicillin.csv"); types = [String, Float64, String, String]) |> DataFrame

@testset "  Model 7: Penicillin.csv parameters                       " begin
    lmm = Metida.LMM(@formula(diameter ~ 1), df;
    random = [Metida.VarEffect(Metida.@covstr(1|plate), Metida.SI), Metida.VarEffect(Metida.@covstr(1|sample), Metida.SI)]
    )
    Metida.fit!(lmm)
    @test Metida.coef(lmm)[1]     ≈ 22.9722 atol = 1E-4
    @test Metida.stderror(lmm)[1] ≈ 0.808573 atol = 1E-4
    @test Metida.theta(lmm)[1]^2 ≈ 0.716908  atol = 1E-4
    @test Metida.theta(lmm)[2]^2 ≈ 3.73092 atol = 1E-4
    @test Metida.theta(lmm)[3]^2 ≈ 0.302415 atol = 1E-4
end


text:
abstract type TemplateStyle end

"""
    struct ImmutableStruct <: TemplateStyle end

Plugin-assembled types marked as `ImmutableStruct` will be generated as a `struct`.
"""
struct ImmutableStruct <: TemplateStyle end

"""
    struct MutableStruct <: TemplateStyle end

Plugin-assembled types marked as `MutableStruct` will be generated as a `mutable struct`.
"""
struct MutableStruct <: TemplateStyle end

"""
    TemplateStyle(::Type) = MutableStruct()

Trait to select the template used for plugin-assembled types

Use [`MutableStruct`](@ref) (default), [`ImmutableStruct`](@ref), or subtype
it when you want to create your own template.

#Examples

Assembling immutable structs:

```julia
abstract type DebugInfo end
Plugins.TemplateStyle(::Type{DebugInfo}) = Plugins.ImmutableStruct()
```

Defining your own template (see also  [`typedef`](@ref) ):

```julia
struct CustomTemplate <: Plugins.TemplateStyle end
Plugins.TemplateStyle(::Type{State}) = CustomTemplate()
```
"""
TemplateStyle(::Type) = MutableStruct()

struct FieldSpec
    name        :: Symbol
    type        :: Type
    constructor :: Union{Function, DataType}
end
"""
    FieldSpec(name, type::Type, constructor::Union{Function, DataType} = type)

Field specification for plugin-assembled types.

Note that every field of an assembled type will be constructed with the same arguments.
 Possibly The constructor will be called when the system
"""
FieldSpec(name, type::Type, constructor::Union{Function, DataType} = type) = FieldSpec(Symbol(name), type, constructor)

@inline structfield(spec::FieldSpec) = :($(spec.name)::$(Meta.parse(string(spec.type))))

struct TypeSpec
    name::Symbol
    parent_type::Type
    fields::Vector{FieldSpec}
    params::Vector{Symbol}
    mod::Module
end

function structfields(spec::TypeSpec)
    return Expr(:block, map(structfield, spec.fields)...)
end

fieldcalls(spec) = map(field -> :($(field.constructor)(args...; kwargs...)), spec.fields)

function default_constructor(spec)
    retval = :(TYPE_NAME(args...; kwargs...) = $(Expr(:call, :TYPE_NAME, fieldcalls(spec)...)))
    return retval
end

"""
    typedef(templatestyle, spec::TypeSpec)::Expr

Return an expression defining a type.

Implement it for your own template styles. More info in the
[tests](https://github.com/tisztamo/Plugins.jl/blob/master/test/customfields.jl).
"""
function typedef end

typedef(::MutableStruct, spec::TypeSpec) = quote
    mutable struct TYPE_NAME{$(spec.params...)} <: $(spec.parent_type)
        $(structfields(spec))
    end;
    $(default_constructor(spec))
    TYPE_NAME
end

typedef(::ImmutableStruct, spec::TypeSpec) = begin
    mutabledef = typedef(MutableStruct(), spec)
    mutabledef.args[2].args[1] = false
    return mutabledef
end

"""
    customtype(stack::PluginStack, typename::Symbol, abstract_type::Type, target_module::Module = Main; unique_name = true)

Assemble a type with fields provided by the plugins in `stack`.

`abstract_type` will be the supertype of the assembled type.

If `unique_name` == `true`, then `typename` will be suffixed with a structure-dependent id.
The  id is generated as a hash of the evaluated expression (with the :TYPE_NAME placeholder
used instead of the name), meaning that for the same id will be generated for a given type
when the same plugins with the same source code are loaded.

# Examples

Assembling a type `AppStateImpl <: AppState` and parametrizing the app with
it.
```julia
abstract type AppState end

mutable struct CustomFieldsApp{TCustomState}
    state::TCustomState
    function CustomFieldsApp(plugins, hookfns, stateargs...)
        stack = PluginStack(plugins, hookfns)
        state_type = customtype(stack, :AppStateImpl, AppState)
        return new{state_type}(Base.invokelatest(state_type, stateargs...))
    end
end
```
!!! note "The need for `invokelatest`"
    We need to use `invokelatest` to instantiate a newly generated type. To use
    the generated type normally, first you have to allow control flow to go to the top-level
    scope after the type was generated. See also the [docs](https://docs.julialang.org/en/v1/manual/methods/#Redefining-Methods)
!!! warning "Antipattern"
    Assembling state types is an antipattern, because plugins can have their own state.
    (This may provide better performance in a few cases though)
    Assembled types can make your code less readable, use them sparingly!
"""
function customtype(
    stack         :: PluginStack,
    typename      :: Symbol,
    abstract_type :: Type = Any,
    params        :: Vector{Symbol} = Symbol[],
    target_module :: Module = Main;
    unique_name   =  true
    )
    hookres = customfields(stack, abstract_type)
    if !hookres.allok
        throw(ErrorException("Cannot define custom type, a plugin throwed an error: $hookres"))
    end
    fields = filter(f -> !isnothing(f), hookres.results)
    spec = TypeSpec(typename, abstract_type, fields, params, target_module)
    def = inject_name(typedef(TemplateStyle(abstract_type), spec), spec, unique_name)
    try
        retval = Base.eval(target_module, def)
    catch e
        @info "Exception while assembling custom type: $def"
        rethrow(e)
    end
end

function inject_name(typedef::Expr, spec::TypeSpec, unique_name)
    suffix = unique_name ? hashstring(typedef) : ""
    name = Symbol(string(spec.name) * "_" * suffix)
    recursive_replace!(typedef, :TYPE_NAME => name)
end

hashstring(typedef::Expr) = string(hash(typedef) % 1679616; base=36) # 4 alphanumeric chars

function recursive_replace!(expr::Expr, source_target_pair)
    replace!(expr.args, source_target_pair)
    foreach(ex -> recursive_replace!(ex, source_target_pair), expr.args)
    return expr
end
recursive_replace!(s, _) = s


text:
"A triangle defined by three points in a plane."
abstract type Triangle{T} <: Domain{SVector{2,T}} end
# This is just a conceptual type that acts as a data container. No functionality is implemented
# to manipulate triangles.

"The triangle below the diagonal of a square `[a,b] × [a,b]`."
struct LowerRightTriangle{T} <: Triangle{T}
   a  :: T
   b  :: T
end

DomainSets.indomain(x, d::LowerRightTriangle) =
   (x[1] >= d.a) && (x[1] <= d.b) && (x[2] >= d.a) && (x[2] <= d.b) && (x[1] <= x[2])

"The triangle above the diagonal of a square `[a,b] × [a,b]`."
struct UpperRightTriangle{T} <: Triangle{T}
   a  :: T
   b  :: T
end

DomainSets.indomain(x, d::UpperRightTriangle) =
   (x[1] >= d.a) && (x[1] <= d.b) && (x[2] >= d.a) && (x[2] <= d.b) && (x[1] >= x[2])

const StaticTypes = DomainSets.StaticTypes


text:
# Copyright (c) 2019 Guido Kraemer
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

using BitConverter, Test

files = readdir(".")

# use only files that are named test_*.jl and start with the newest

filter!(x -> occursin(r"^test_.*\.jl$", x), files)
sort!(files, by = x -> stat(x).mtime, rev = true)

for f in files
    include(f)
end


text:
@testset "363.max-sum-of-rectangle-no-larger-than-k.jl" begin
    @test max_sum_submatrix(hcat([[1, 0, 1], [0, -2, 3]]...), 2) == 2
    @test max_sum_submatrix(hcat([[2, 2, -1]]...), 3) == 3
end


text:
module RocketLastOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: last()" begin

    println("Testing: operator last()")

    run_proxyshowcheck("Last", last())

    run_testset([

        (
            source = from(1:42) |> last(),
            values = @ts([ 42, c ])
        ),
        (
            source = timer(50, 10) |> take(10) |> last(),
            values = @ts(150 ~ [ 9, c ])
        ),
        (
            source = completed() |> last(),
            values = @ts(e(LastNotFoundException()))
        ),
        (
            source      = completed(Int) |> last(default = "String"),
            values      = @ts([ "String", c ]),
            source_type = Union{Int, String}
        ),
        (
            source      = faulted(Int, "e") |> last(),
            values      = @ts(e("e")),
            source_type = Union{Int}
        ),
        (
            source      = faulted(Int, "e") |> last(default = "String"),
            values      = @ts(e("e")),
            source_type = Union{Int, String}
        ),
        (
            source = never() |> last(),
            values = @ts()
        )
    ])

end

end


text:
module CUBLAS

import CUDAdrv: CUDAdrv, CuContext, CuStream_t

using ..CuArrays
using ..CuArrays: libcublas, active_context

using LinearAlgebra

include("libcublas_types.jl")
include("error.jl")
include("libcublas.jl")

const _handles = Dict{CuContext,cublasHandle_t}()
const _handle = Ref{cublasHandle_t}(C_NULL)

function handle()
    if _handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cublasCreate_v2()
            atexit(()->CUDAdrv.isvalid(context) && cublasDestroy_v2(handle))
            handle
        end
    end

    return _handle[]
end

include("util.jl")
include("wrappers.jl")
include("highlevel.jl")

end


text:
using SparseVectors
using Base.Test

import SparseVectors: exact_equal

### Data

rnd_x0 = sprand(50, 0.6)
rnd_x0f = full(rnd_x0)

rnd_x1 = sprand(50, 0.7) * 4.0
rnd_x1f = full(rnd_x1)

spv_x1 = SparseVector(8, [2, 5, 6], [1.25, -0.75, 3.5])
_x2 = SparseVector(8, [1, 2, 6, 7], [3.25, 4.0, -5.5, -6.0])
spv_x2 = view(_x2)

### Arithmetic operations

let x = spv_x1, x2 = spv_x2
    # negate
    @test exact_equal(-x, SparseVector(8, [2, 5, 6], [-1.25, 0.75, -3.5]))

    # abs and abs2
    @test exact_equal(abs(x), SparseVector(8, [2, 5, 6], abs([1.25, -0.75, 3.5])))
    @test exact_equal(abs2(x), SparseVector(8, [2, 5, 6], abs2([1.25, -0.75, 3.5])))

    # plus and minus
    xa = SparseVector(8, [1,2,5,6,7], [3.25,5.25,-0.75,-2.0,-6.0])

    @test exact_equal(x + x, x * 2)
    @test exact_equal(x + x2, xa)
    @test exact_equal(x2 + x, xa)

    xb = SparseVector(8, [1,2,5,6,7], [-3.25,-2.75,-0.75,9.0,6.0])
    @test exact_equal(x - x, SparseVector(8, Int[], Float64[]))
    @test exact_equal(x - x2, xb)
    @test exact_equal(x2 - x, -xb)

    @test full(x) + x2 == full(xa)
    @test full(x) - x2 == full(xb)
    @test x + full(x2) == full(xa)
    @test x - full(x2) == full(xb)

    # multiplies
    xm = SparseVector(8, [2, 6], [5.0, -19.25])
    @test exact_equal(x .* x, abs2(x))
    @test exact_equal(x .* x2, xm)
    @test exact_equal(x2 .* x, xm)

    @test full(x) .* x2 == full(xm)
    @test x .* full(x2) == full(xm)

    # max & min
    @test exact_equal(max(x, x), x)
    @test exact_equal(min(x, x), x)
    @test exact_equal(max(x, x2),
        SparseVector(8, Int[1, 2, 6], Float64[3.25, 4.0, 3.5]))
    @test exact_equal(min(x, x2),
        SparseVector(8, Int[2, 5, 6, 7], Float64[1.25, -0.75, -5.5, -6.0]))
end

### Complex

let x = spv_x1, x2 = spv_x2
    # complex
    @test exact_equal(complex(x, x),
        SparseVector(8, [2,5,6], [1.25+1.25im, -0.75-0.75im, 3.5+3.5im]))
    @test exact_equal(complex(x, x2),
        SparseVector(8, [1,2,5,6,7], [3.25im, 1.25+4.0im, -0.75+0.im, 3.5-5.5im, -6.0im]))
    @test exact_equal(complex(x2, x),
        SparseVector(8, [1,2,5,6,7], [3.25+0.im, 4.0+1.25im, -0.75im, -5.5+3.5im, -6.0+0.im]))

    # real & imag

    @test is(real(x), x)
    @test exact_equal(imag(x), sparsevector(Float64, length(x)))

    xcp = complex(x, x2)
    @test exact_equal(real(xcp), x)
    @test exact_equal(imag(xcp), x2)
end

### Zero-preserving math functions: sparse -> sparse

function check_nz2z_z2z{T}(f::Function, x::SparseVector{T}, xf::Vector{T})
    R = typeof(f(zero(T)))
    r = f(x)
    isa(r, AbstractSparseVector) || error("$f(x) is not a sparse vector.")
    eltype(r) == R || error("$f(x) results in eltype = $(eltype(r)), expect $R")
    all(r.nzval .!= 0) || error("$f(x) contains zeros in nzval.")
    full(r) == f(xf) || error("Incorrect results found in $f(x).")
end

for f in [floor, ceil, trunc, round]
    check_nz2z_z2z(f, rnd_x1, rnd_x1f)
end

for f in [log1p, expm1,
          sin, tan, sinpi, sind, tand,
          asin, atan, asind, atand,
          sinh, tanh, asinh, atanh]
    check_nz2z_z2z(f, rnd_x0, rnd_x0f)
end

### Non-zero-preserving math functions: sparse -> dense

function check_z2nz{T}(f::Function, x::SparseVector{T}, xf::Vector{T})
    R = typeof(f(zero(T)))
    r = f(x)
    isa(r, Vector) || error("$f(x) is not a dense vector.")
    eltype(r) == R || error("$f(x) results in eltype = $(eltype(r)), expect $R")
    r == f(xf) || error("Incorrect results found in $f(x).")
end

for f in [exp, exp2, exp10, log, log2, log10,
          cos, csc, cot, sec, cospi,
          cosd, cscd, cotd, secd,
          acos, acot, acosd, acotd,
          cosh, csch, coth, sech, acsch, asech]
    check_z2nz(f, rnd_x0, rnd_x0f)
end

### Reduction

# sum, sumabs, sumabs2, vecnorm

let x = spv_x1
    @test sum(x) == 4.0
    @test sumabs(x) == 5.5
    @test sumabs2(x) == 14.375

    @test vecnorm(x) == sqrt(14.375)
    @test vecnorm(x, 1) == 5.5
    @test vecnorm(x, 2) == sqrt(14.375)
    @test vecnorm(x, Inf) == 3.5
end

# maximum, minimum, maxabs, minabs

let x = spv_x1
    @test maximum(x) == 3.5
    @test minimum(x) == -0.75
    @test maxabs(x) == 3.5
    @test minabs(x) == 0.0
end

let x = abs(spv_x1)
    @test maximum(x) == 3.5
    @test minimum(x) == 0.0
end

let x = -abs(spv_x1)
    @test maximum(x) == 0.0
    @test minimum(x) == -3.5
end

let x = SparseVector(3, [1, 2, 3], [-4.5, 2.5, 3.5])
    @test maximum(x) == 3.5
    @test minimum(x) == -4.5
    @test maxabs(x) == 4.5
    @test minabs(x) == 2.5
end

let x = sparsevector(Float64, 8)
    @test maximum(x) == 0.0
    @test minimum(x) == 0.0
    @test maxabs(x) == 0.0
    @test minabs(x) == 0.0
end

# Issue https://github.com/JuliaLang/julia/issues/14046
let s14046 = sprand(5, 1.0)
    z = sparsevector(Float64, 5)
    @test z + s14046 == s14046
    @test 2*z == z + z == z
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function busywait(t, n_tries)
    iter = 0
    while iter < n_tries && Profile.len_data() == 0
        iter += 1
        tend = time() + t
        while time() < tend end
    end
end

Profile.clear()
@profile busywait(1, 20)
let iobuf = IOBuffer()
    Profile.print(iobuf, format=:tree, C=true)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:tree, maxdepth=2)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:flat, C=true)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf)
    @test !isempty(String(take!(iobuf)))
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:flat, sortedby=:count)
    @test !isempty(String(take!(iobuf)))
    Profile.clear()
    @test isempty(Profile.fetch())
end


text:
@safetestset "Assets functionality" begin

  @safetestset "Assets paths" begin
    using Genie, Genie.Assets

    @test include_asset(:css, "foo")  == "/assets/css/foo.css"
    @test include_asset(:js, "foo")   == "/assets/js/foo.js"

    @test css_asset("foo") == css("foo") == "/assets/css/foo.css"
    @test js_asset("foo") == js("foo") == "/assets/js/foo.js"
  end;

  @safetestset "Expose settings" begin
    using Genie, Genie.Assets

    @test js_settings() == "window.Genie = {};\nGenie.Settings = {\"webchannels_autosubscribe\":true,\"server_host\":\"127.0.0.1\",\"webchannels_eval_command\":\">eval:\",\"webthreads_js_file\":\"webthreads.js\",\"webchannels_unsubscribe_channel\":\"unsubscribe\",\"webthreads_default_route\":\"__\",\"webchannels_subscribe_channel\":\"subscribe\",\"server_port\":$(Genie.config.server_port),\"base_path\":\"$(Genie.config.base_path)\",\"webthreads_pull_route\":\"pull\",\"webchannels_default_route\":\"__\",\"webchannels_timeout\":1000,\"webthreads_push_route\":\"push\",\"websockets_port\":$(Genie.config.websockets_port)};\n"
  end

  @safetestset "Embedded assets" begin
    using Genie, Genie.Assets

    @test Assets.channels()[1:18] == "window.Genie = {};"
    @test channels_script()[1:27] == "<script>\nwindow.Genie = {};"

    @test channels_support() == "<script src=\"/genie.jl/master/assets/js/__/channels.js\">\n\n</script>\n"
    @test Genie.Router.routes()[1].path == "/genie.jl/master/assets/js/__/channels.js"
    @test Genie.Router.channels()[1].path == "/__/unsubscribe"
    @test Genie.Router.channels()[2].path == "/__/subscribe"

    @test favicon_support() == "<link rel=\"icon\" type=\"image/x-icon\" href=\"/favicon.ico\" />"
  end

end;


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

JLLWrappers.@generate_main_file_header("OpenSlide")
JLLWrappers.@generate_main_file("OpenSlide", UUID("becf559c-afb9-529d-9a2d-92566e0304eb"))
end  # module OpenSlide_jll


text:
using Base.Test
using Einsum

## Test that vars in Main aren't overwritten by einsum
i = -1
y = randn(10)
@einsum x[i] := y[i]
@test i == -1

## Test with preallocated array ##

A = zeros(5,6,7);
X = randn(5,2);
Y = randn(6,2);
Z = randn(7,2);

@einsum A[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]

for i = 1:5
    for j = 1:6
        for k = 1:7
            s = 0
            for r = 1:2
                s += X[i,r]*Y[j,r]*Z[k,r]
            end
            @test A[i,j,k] == s
        end
    end
end

## Test without preallocated array ##

@einsum A2[i,j,k] := X[i,r]*Y[j,r]*Z[k,r]
@test isapprox(A,A2)

# Interesting test case, can throw an error that
# local vars are declared twice. Solution was to wrap
# everything in a let statement.
if true
    @einsum A[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]
else
    @einsum A[i,j,k] = X[i,r]*Y[j,r]*Z[k,r]
end

# At one point this threw an error because the lhs
# had no indices/arguments
x = randn(10)
y = randn(10)
@einsum k := x[i]*y[i]
@test isapprox(k,dot(x,y))

# Elementwise multiplication (this should create nested loops with no
# no summation, due to lack of repeated variables.)
x = randn(10)
y = randn(10)
@einsum k[i] := x[i]*y[i]
@test isapprox(k,x.*y)

# Transpose a block matrix
z = Any[ rand(2,2) for i=1:2, j=1:2]
@einsum t[i,j] := transpose(z[j,i])
@test isapprox(z[1,1], t[1,1]')
@test isapprox(z[2,2], t[2,2]')
@test isapprox(z[1,2], t[2,1]')
@test isapprox(z[2,1], t[1,2]')

# Example from numpy
A = reshape(collect(1:25),5,5)
@einsum B[i] := A[i,i]
@test all(B .== [1,7,13,19,25])


text:
module SoapySDR

using soapysdr_jll
const lib = soapysdr_jll.libsoapysdr

using Intervals
using Unitful
using Unitful.DefaultSymbols
const dB = u"dB"
const GC = Base.GC

export @u_str

include("lowlevel/Constants.jl") # Done
include("lowlevel/Errors.jl")    # Done
include("lowlevel/Formats.jl")   # Done
include("lowlevel/Types.jl")     # Done
include("lowlevel/Device.jl")
include("lowlevel/Modules.jl")
include("lowlevel/Logger.jl")
include("typemap.jl")
include("typewrappers.jl")
include("highlevel.jl")
include("loghandler.jl")

const SDRStream = Stream
export SDRStream

#
# FIXME: Log handler breaks USRP Support
#
#function __init__()
#    julia_log_handler = @cfunction(logger_soapy2jl, Cvoid, (Cint, Cstring))
#    SoapySDR_registerLogHandler(julia_log_handler)
#end

end


text:
###############################################################################
#
# Canopy optical variables
#
###############################################################################
"""
    mutable struct CanopyOpticals{FT}

A struct for canopy optical properties

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct CanopyOpticals{FT}
    # local storage of dimension information
    "Number of azimuth angles"
    nAzi  ::Int = 36
    "Number of inclination agles"
    nIncl ::Int = 9
    "Number of canopy layers"
    nLayer::Int = 5
    "Number of wave lengths"
    nWL   ::Int = 10

    # Single value
    "Solar -> Diffuse backscatter weight"
    sdb::FT = FT(0.0)
    "Solar -> Diffuse forward scatter weight"
    sdf::FT = FT(0.0)
    "Diffuse -> Directional backscatter weight"
    dob::FT = FT(0.0)
    "Diffuse -> Directional forward scatter weight"
    dof::FT = FT(0.0)
    "Diffuse -> Diffuse backscatter weight"
    ddb::FT = FT(0.0)
    "Diffuse -> Diffuse forward scatter weight"
    ddf::FT = FT(0.0)
    "Solar beam extinction coefficient weight"
    ks ::FT = FT(0.0)
    "Outgoing beam extinction coefficient weight"
    ko ::FT = FT(0.0)
    # TODO What is this?
    "?"
    bf ::FT = FT(0.0)
    "Weight of specular2directional backscatter coefficient"
    sob::FT = FT(0.0)
    "Weight of specular2directional forward coefficient"
    sof::FT = FT(0.0)

    # dimension of nLevel
    "Probability of directly viewing a leaf in solar direction"
    Ps ::Array{FT,1} = zeros(FT, nLayer+1)
    "Probability of directly viewing a leaf in viewing direction"
    Po ::Array{FT,1} = zeros(FT, nLayer+1)
    "Bi-directional probability of directly viewing a leaf (solar->canopy->viewing)"
    Pso::Array{FT,1} = zeros(FT, nLayer+1)

    # dimension of nIncl * nAzi
    "conversion factor fs to compute irradiance on inclined leaf"
    fs     ::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "abs(fs)"
    absfs  ::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "abs(fs*fo)"
    absfsfo::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "fs*fo"
    fsfo   ::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "conversion factor fo for angle towards observer (not sun like fs)"
    fo     ::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "abs(fo)"
    absfo  ::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "Cosine of leaf azimuths"
    cosΘ_l ::Array{FT,2} = zeros(FT, (nIncl, nAzi))
    "cos of leaf azimuth sqared"
    cos2Θ_l::Array{FT,2} = zeros(FT, (nIncl, nAzi))

    # The following also depend on leaf reflectance and transmission.
    # Might go into a separate strcuture so that we can have it separately for thermal, SW and SIF?
    # dimension of nWL * nLayer
    "diffuse     backscatter scattering coefficient for diffuse  incidence"
    sigb::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "diffuse     forward     scattering coefficient for diffuse  incidence"
    sigf::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "diffuse     backscatter scattering coefficient for specular incidence"
    sb  ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "diffuse     forward     scattering coefficient for specular incidence"
    sf  ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "directional backscatter scattering coefficient for diffuse  incidence"
    vb  ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "directional forward     scattering coefficient for diffuse  incidence"
    vf  ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "bidirectional scattering coefficent (directional-directional)"
    w   ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "attenuation"
    a   ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "Effective layer transmittance (direct->diffuse)"
    Xsd ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "Effective layer transmittance (diffuse->diffuse)"
    Xdd ::Array{FT,2} = zeros(FT, (nWL,nLayer))

    # dimension of nWL * nLevel
    "Effective layer reflectance (direct->diffuse)"
    R_sd::Array{FT,2} = zeros(FT, (nWL, nLayer+1))
    "Effective layer reflectance (diffuse->diffuse)"
    R_dd::Array{FT,2} = zeros(FT, (nWL, nLayer+1))
    "Solar direct radiation per layer)"
    Es_::Array{FT,2}  = zeros(FT, (nWL, nLayer+1))
end

max_stars_repo_name: climate-machine/Land


text:
# src path
# srcまでの絶対パスを記入
src_path="C:\\Users\\xxxxx\\Desktop\\git\\general_2d_diff\\src\\"

# main (変更しないこと)
src_read="read_grid.jl"
include(src_path*src_read)
src_read="read_para.jl"
include(src_path*src_read)
src_read="output.jl"
include(src_path*src_read)
src_read="boundary.jl"
include(src_path*src_read)
src_read="misc.jl"
include(src_path*src_read)
src_read="value_setup.jl"
include(src_path*src_read)
src_read="rhs.jl"
include(src_path*src_read)
src_read="cal_time_step.jl"
include(src_path*src_read)
src_main="main.jl"
include(src_path*src_main)


text:
struct Dataset{T<:Real}

    X::AbstractMatrix{T}
    y::AbstractVector{T}
    n::Int
    nfeatures::Int
    weighted::Bool
    weights::Union{AbstractVector{T}, Nothing}
    varMap::Array{String, 1}

end

"""
    Dataset(X::AbstractMatrix{T}, y::AbstractVector{T};
            weights::Union{AbstractVector{T}, Nothing}=nothing,
            varMap::Union{Array{String, 1}, Nothing}=nothing)

Construct a dataset to pass between internal functions.
"""
function Dataset(
        X::AbstractMatrix{T},
        y::AbstractVector{T};
        weights::Union{AbstractVector{T}, Nothing}=nothing,
        varMap::Union{Array{String, 1}, Nothing}=nothing
       ) where {T<:Real}

    n = size(X, 2)
    nfeatures = size(X, 1)
    weighted = weights !== nothing
    if varMap == nothing
        varMap = ["x$(i)" for i=1:nfeatures]
    end

    return Dataset{T}(X, y, n, nfeatures, weighted, weights, varMap)

end


text:
module MetaGraphsNext

using JLD2
using LightGraphs

import Base:
    Tuple,
    ==,
    copy,
    delete!,
    eltype,
    getindex,
    haskey,
    issubset,
    length,
    setindex!,
    show,
    size,
    zero

import LightGraphs:
    AbstractGraph,
    AbstractGraphFormat,
    add_edge!,
    add_vertex!,
    edges,
    edgetype,
    has_edge,
    has_vertex,
    induced_subgraph,
    inneighbors,
    is_directed,
    loadgraph,
    ne,
    nv,
    outneighbors,
    rem_edge!,
    rem_vertex!,
    reverse,
    savegraph,
    vertices,
    weights

import LightGraphs.SimpleGraphs: SimpleGraph, SimpleDiGraph, fadj, badj

export MetaGraph, weighttype, defaultweight, weightfunction
export MGFormat, DOTFormat
export reverse
export label_for, code_for

include("metagraph.jl")
include("metaundigraph.jl")
include("metadigraph.jl")
include("lightgraphs.jl")
include("weights.jl")
include("dict_utils.jl")
include("overrides.jl")
include("persistence.jl")

end # module


text:
using CUDAnative, CUDAdrv
import LLVM

using Test

@testset "CUDAnative" begin

include("util.jl")

include("base.jl")
include("pointer.jl")
include("codegen.jl")

if CUDAnative.configured
    @test length(devices()) > 0
    if length(devices()) > 0
        # the API shouldn't have been initialized
        @test CuCurrentContext() == nothing

        device_callbacked = nothing
        device_callback = (dev, ctx) -> begin
            device_callbacked = dev
        end
        push!(CUDAnative.device!_listeners, device_callback)

        # now cause initialization
        Mem.alloc(1)
        @test CuCurrentContext() != nothing
        @test device(CuCurrentContext()) == CuDevice(0)
        @test device_callbacked == CuDevice(0)

        device!(CuDevice(0))
        device!(CuDevice(0)) do
            nothing
        end

        # test the device selection functionality
        if length(devices()) > 1
            device!(1) do
                @test device(CuCurrentContext()) == CuDevice(1)
            end
            @test device(CuCurrentContext()) == CuDevice(0)

            device!(1)
            @test device(CuCurrentContext()) == CuDevice(1)
        end

        # pick most recent device (based on compute capability)
        global dev = first(sort(collect(devices()); by=capability))
        @info("Testing using device $(name(dev))")
        device!(dev)

        if capability(dev) < v"2.0"
            @warn("native execution not supported on SM < 2.0")
        else
            include("device/codegen.jl")
            include("device/execution.jl")
            include("device/pointer.jl")
            include("device/array.jl")
            include("device/intrinsics.jl")

            include("examples.jl")
        end
    end
else
    @warn("CUDAnative.jl has not been configured; skipping on-device tests.")
end

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

"""
    is_unix([os])

Predicate for testing if the OS provides a Unix-like interface.
See documentation in [Handling Operating System Variation](@ref).
"""
function is_unix(os::Symbol)
    if is_windows(os)
        return false
    elseif is_linux(os) || is_bsd(os)
        return true
    else
        throw(ArgumentError("unknown operating system \"$os\""))
    end
end

"""
    is_linux([os])

Predicate for testing if the OS is a derivative of Linux.
See documentation in [Handling Operating System Variation](@ref).
"""
is_linux(os::Symbol) = (os == :Linux)

"""
    is_bsd([os])

Predicate for testing if the OS is a derivative of BSD.
See documentation in [Handling Operating System Variation](@ref).
"""
is_bsd(os::Symbol) = (os == :FreeBSD || os == :OpenBSD || os == :NetBSD || os == :DragonFly || os == :Darwin || os == :Apple)

"""
    is_windows([os])

Predicate for testing if the OS is a derivative of Microsoft Windows NT.
See documentation in [Handling Operating System Variation](@ref).
"""
is_windows(os::Symbol) = (os == :Windows || os == :NT)

"""
    is_apple([os])

Predicate for testing if the OS is a derivative of Apple Macintosh OS X or Darwin.
See documentation in [Handling Operating System Variation](@ref).
"""
is_apple(os::Symbol) = (os == :Apple || os == :Darwin)

"""
    @static

Partially evaluates an expression at parse time.

For example, `@static is_windows() ? foo : bar` will evaluate `is_windows()` and insert either `foo` or `bar` into the expression.
This is useful in cases where a construct would be invalid on other platforms,
such as a `ccall` to a non-existent function.
`@static if is_apple() foo end` and `@static foo <&&,||> bar` are also valid syntax.
"""
macro static(ex)
    if isa(ex, Expr)
        if ex.head === :if || ex.head === :&& || ex.head === :||
            cond = eval(current_module(), ex.args[1])
            if xor(cond, ex.head === :||)
                return esc(ex.args[2])
            elseif length(ex.args) == 3
                return esc(ex.args[3])
            elseif ex.head === :if
                return nothing
            else
                return cond
            end
        end
    end
    throw(ArgumentError("invalid @static macro"))
end

let KERNEL = ccall(:jl_get_UNAME, Any, ())
    # evaluate the zero-argument form of each of these functions
    # as a function returning a static constant based on the build-time
    # operating-system kernel
    for f in (:is_unix, :is_linux, :is_bsd, :is_apple, :is_windows)
        @eval $f() = $(getfield(current_module(),f)(KERNEL))
    end
end


text:
#=
    Generating Hypergeometric distribution

    https://en.wikipedia.org/wiki/Hypergeometric_distribution
    """
    [T]he probability of k successes (random draws for which the object
    drawn has a specified feature) in n draws, without replacement, from
    a finite population of size N that contains exactly K objects with
    that feature, wherein each draw is either a success or a failure.
    In contrast, the binomial distribution describes the probability of
    k successes in n draws with replacement.
    """

    The example from
    https://en.wikipedia.org/wiki/Hypergeometric_distribution#Working_example
    """
    Now, assume (for example) that there are 5 green and 45 red marbles in
    the urn. Standing next to the urn, you close your eyes and draw 10 marbles
    without replacement. What is the probability that exactly 4 of the 10
    are green?

    The probability of drawing exactly 4 marbels is: ~0.00396458305801506542
    The probability of drawing exactly 5 marbels is: ~0.00011893749174045196
    """

    Since the probabilities are so small, we use big() in the calculations.
    Note that gpmf and gmean are the exact values from the formulas,
    while g and gCount are simulated from the function hypergeometric(.).

    * k = 4
        Summary Statistics
        parameters       mean        std   naive_se       mcse          ess       rhat   ess_per_sec
            Symbol   BigFloat   BigFloat   BigFloat   BigFloat     BigFloat   BigFloat      BigFloat

                g     0.0036     0.0601     0.0003     0.0003   49813.4280     1.0000    97291.8515
             gpmf     0.0040     0.0000     0.0000     0.0000          NaN        NaN           NaN
           gCount     0.9935     0.8567     0.0038     0.0038   50059.3725     1.0000    97772.2120
            gmean     1.0000     0.0000     0.0000     0.0000     100.7029     1.0000      196.6854
             gvar     0.0000     0.0000     0.0000     0.0000          NaN        NaN           NaN

        Distributions of variable gCount (num:0)
        1.00000 =>   21503  (0.430060)
        0.00000 =>   15703  (0.314060)
        2.00000 =>   10430  (0.208600)
        3.00000 =>    2155  (0.043100)
        4.00000 =>     200  (0.004000)
        5.00000 =>       9  (0.000180)

        mean(g): BigFloat[0.003619999999999999999999999999999999999999999999999999999999999999999999999999983]
        mean(gpmf): BigFloat[0.003964583058015065664125842204157379455864429473876953125]
        mean(gmean): BigFloat[1.0]
        mean(gvar): BigFloat[0.0]

    * k = 5

        parameters       mean        std   naive_se       mcse          ess       rhat   ess_per_sec
            Symbol   BigFloat   BigFloat   BigFloat   BigFloat     BigFloat   BigFloat      BigFloat

                g     0.0001     0.0110     0.0000     0.0000   50014.6732     1.0000   106414.1984
             gpmf     0.0001     0.0000     0.0000     0.0000     100.7029     1.0000      214.2615
           gCount     0.9991     0.8551     0.0038     0.0037   50011.1699     1.0000   106406.7445
            gmean     1.0000     0.0000     0.0000     0.0000     100.7029     1.0000      214.2615
             gvar     0.0000     0.0000     0.0000     0.0000          NaN        NaN           NaN

        Distributions of variable gCount (num:0)
        1.00000 =>   21567  (0.431340)
        0.00000 =>   15521  (0.310420)
        2.00000 =>   10566  (0.211320)
        3.00000 =>    2134  (0.042680)
        4.00000 =>     208  (0.004160)
        5.00000 =>       4  (0.000080)

        mean(g): BigFloat[0.0001200000000000000000000000000000000000000000000000000000000000000000000000000003]
        mean(gpmf): BigFloat[0.0001189374917404519634185622312116947796312160789966583251953125]
        mean(gmean): BigFloat[1.0]
        mean(gvar): BigFloat[0.0]

   Cf https://github.com/distributions-io/hypergeometric-random/blob/master/lib/number.js

   See also hypergeometric_dist.wppl
=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

#=
   Hypergeometric:
   What is the probability that we draw exactly k "success" objects
   of the n drawn objects of total N objects where there are in total
   K "success" objects

   k: number of successes we want to check
   N: total number of objects
   K: total number of success objects
   n: number of draws

=#
function hypergeometric1(k,N,K,n,count)
    if n == 0 || K <= 0
        return count
    else
        #   we have K successes left and N objects left
        p = K/N #  probability of drawing a success object
        ff = rand(Bernoulli(p))
        if ff == 1.0
            #  We drew a success:
            #  - decrement the total objects (N)
            #  - decrement the number of "success" objects (K)
            #  - decrement the number of drawn objects (n)
            #  - increment the number of successful draws (count)
            return hypergeometric1(k,N-1,K-1,n-1,count+1)
        else
            #  We drew a failure:
            #  - decrement the total objects (N)
            #  - decrement the number of drawn objects (n)
            return hypergeometric1(k,N-1,K,n-1,count)
        end
    end

end

function hypergeometric(k,N,K,n)
    res = hypergeometric1(k,N,K,n,0)
    return res == k
end

#
# Return the number of found successes.
#
function hypergeometricCount(k,N,K,n)
    res = hypergeometric1(k,N,K,n,0)
    return res
end

# PDF, From https://en.wikipedia.org/wiki/Hypergeometric_distribution
# Gives the exact probability mass function
function hypergeometric_pmf(N,K,n,k)
    return  binomial(K,k)*binomial(N-K,n-k) / binomial(N,n)
end

# Exact probability
# From https://en.wikipedia.org/wiki/Hypergeometric_distribution
function hypergeometric_mean(N,K,n)
    n*K/N
end

# From https://en.wikipedia.org/wiki/Hypergeometric_distribution
function hypergeometric_variance(N,K,n)
    # n{K \over N}{(N-K) \over N}{N-n \over N-1}
    n*binomial(K,N)*binomial(N-K,N)*binomial(N-n,N-1)
end

@model function hypergeometric_dist_test()
    #  total: 5 green and 45 red marbles
    #  drawn: 4 green marbles, 6 red marbles
    K = 5 #  total green marbles: 4 drawn + 1 not drawn
    N = 50 #  total marbles: 5 green + 45 red marbles

    k = 4 #  drawn green_marbles
    # k = 5 #  drawn green_marbles
    n = 10 #  total drawn green + red marbles

    # Note: we use big float since the probabilities
    # are so small.
    g ~ Dirac(big(hypergeometric(k,N,K,n)))
    gpmf ~ Dirac(big(hypergeometric_pmf(N,K,n,k))) # Exact
    gCount ~ Dirac(hypergeometricCount(k,N,K,n)) #  Count version
    gmean ~ Dirac(big(hypergeometric_mean(N,K,n)))
    gvar  ~ Dirac(big(hypergeometric_variance(N,K,n)))
end

model = hypergeometric_dist_test()

# We can use Prior since there's not any observations.
chns = sample(model, Prior(), 50_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5),  10_000)
# chns = sample(model, PG(5),  MCMCThreads(), 10_000, 4)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

show_var_dist_pct(chns,:gCount)
println()
println("mean(g): $(mean(chns[[:g]].value.data,dims=1))")
println("mean(gpmf): $(mean(chns[[:gpmf]].value.data,dims=1))")
println("mean(gmean): $(mean(chns[[:gmean]].value.data,dims=1))")
println("mean(gvar): $(mean(chns[[:gvar]].value.data,dims=1))")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
    fS(eqn::Euler{d},UL,UR) where {d}

Entropy conservative fluxes for the compressible Euler equations, where UL,UR are
tuples of left and right states in terms of the conservative variables Q = (ρ,ρuvw...,E)
"""
@inline function fS(eqn::Euler{d},UL,UR) where {d}
    QL = cons_to_prim_beta(eqn,UL)
    QR = cons_to_prim_beta(eqn,UR)
    return fS_prim(eqn,QL,QR)
end

"""
    fS_prim(eqn::Euler{d},UL,UR)

Entropy conservative fluxes for the compressible Euler equations, where
UL,UR are tuples of left and right solutions in terms of β-primitive variables
Q = (rho,uvw...,β) (where uvw = velocity components and β is inverse temperature).
"""
@inline function fS_prim(eqn::Euler{d},UL,UR) where {d}
    logL = map(x->log.(x),(first(UL),last(UL)))
    logR = map(x->log.(x),(first(UR),last(UR)))

    return fS_prim_log(eqn,(UL...,logL...),(UR...,logR...))
end

# manually unroll for d = 1,2,3 to avoid type instability
@inline function fS_prim(eqn::Euler{1},UL,UR)
    UlogL = (UL[1],UL[2],UL[3],log.(UL[1]),log.(UL[3]))
    UlogR = (UR[1],UR[2],UR[3],log.(UR[1]),log.(UR[3]))
    return fS_prim_log(eqn,UlogL,UlogR)
end

@inline function fS_prim(eqn::Euler{2},UL,UR)
    UlogL = (UL[1],UL[2],UL[3],UL[4],log.(UL[1]),log.(UL[4]))
    UlogR = (UR[1],UR[2],UR[3],UR[4],log.(UR[1]),log.(UR[4]))
    return fS_prim_log(eqn,UlogL,UlogR)
end

@inline function fS_prim(eqn::Euler{3},UL,UR)
    UlogL = (UL[1],UL[2],UL[3],UL[4],UL[5],log.(UL[1]),log.(UL[5]))
    UlogR = (UR[1],UR[2],UR[3],UR[4],UR[5],log.(UR[1]),log.(UR[5]))
    return fS_prim_log(eqn,UlogL,UlogR)
end

"""
    fS_prim_log(eqn::Euler{d},UL,UR)

Entropy conservative fluxes for the compressible Euler equations, where
UL,UR are tuples of left and right solutions (in β-primitive variables) and logs
of rho, β, e.g., U = (rho,uvw...,β,log.(rho),log.(β)).
"""
@inline function fS_prim_log(eqn::Euler{1},UL,UR)

    @unpack γ = eqn

    rhoL,uL,betaL,rhologL,betalogL = UL
    rhoR,uR,betaR,rhologR,betalogR = UR

    rholog  = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))

    unorm = (@. uL*uR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. f4aux*uavg)

    return SVector(FxS1,FxS2,FxS3)
end

@inline function fS_prim_log(eqn::Euler{2},UL,UR)

    @unpack γ = eqn

    rhoL,uL,vL,betaL,rhologL,betalogL = UL
    rhoR,uR,vR,betaR,rhologR,betalogR = UR

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg)
    FxS4 = (@. f4aux*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3)
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. f4aux*vavg)

    return SVector(FxS1,FxS2,FxS3,FxS4),SVector(FyS1,FyS2,FyS3,FyS4)
end

@inline function fS_prim_log(eqn::Euler{3},UL,UR)

    @unpack γ = eqn

    rhoL,uL,vL,wL,betaL,rhologL,betalogL = UL
    rhoR,uR,vR,wR,betaR,rhologR,betalogR = UR

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))
    wavg   = (@. .5*(wL+wR))

    unorm = (@. uL*uR + vL*vR + wL*wR)
    pa    = (@. rhoavg/(betaL+betaR))
    E_plus_p  = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg) # rho * u * v
    FxS4 = (@. FxS1*wavg) # rho * u * w
    FxS5 = (@. E_plus_p*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3) # rho * u * v
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. FyS1*wavg) # rho * v * w
    FyS5 = (@. E_plus_p*vavg)

    FzS1 = (@. rholog*wavg)
    FzS2 = (@. FxS4) # rho * u * w
    FzS3 = (@. FyS4) # rho * v * w
    FzS4 = (@. FzS1*wavg + pa) # rho * w^2 + p
    FzS5 = (@. E_plus_p*wavg)

    Fx = SVector(FxS1,FxS2,FxS3,FxS4,FxS5)
    Fy = SVector(FyS1,FyS2,FyS3,FyS4,FyS5)
    Fz = SVector(FzS1,FzS2,FzS3,FzS4,FzS5)
    return Fx,Fy,Fz
end

### didn't add coordinates b/c timing doesn't seem to make a diff.
# julia> Qlog = (1.0, 0.1, 0.2, 2.5250000000000004, 0.0, 0.9262410627273233)
# julia> @btime fS_prim_log($Euler{2}(),$Qlog,$Qlog);
#   2.963 μs (0 allocations: 0 bytes)
#
# julia> @btime fS_prim_log($Euler{2}(),$Qlog,$Qlog,$1);
#   3.044 μs (0 allocations: 0 bytes)


text:
#include("helperfunctions_rectgrid.jl")

using RecipesBase

@recipe function f(rp::RectangularPartition)
    label --> ""
    for rect in rp.rectangles
        @series begin
            connectrect(rect)
        end
    end
end

@recipe function f(r::Rectangle)
    label --> ""
    lc --> :black
    connectrect(r)
end