text:
using ArgCheck
using Distributions
using Parameters
using HamiltonianABC
using StatsBase
using StatPlots                 # for kernel density
import HamiltonianABC: logdensity, simulate!
using Base.Test

using Plots
plotlyjs()                      # Tamas likes this, optional

"""
    simulate_stochastic(ρ, σ_v, ϵs, νs)

Take in the parameter values (ρ, σ) for the latent volatility process, the errors ϵs used for the measurement equation and the errors νs used for the transition equation.

The discrete-time version of the Ornstein-Ulenbeck Stochastic - volatility model:

    y_t = x_t + ϵ_t where ϵ_t ∼ χ^2(1)
    x_t = ρ * x_(t-1) + σ * ν_t  where ν_t ∼ N(0, 1)

"""
function simulate_stochastic(ρ, σ, ϵs, νs)
    N = length(ϵs)
    @argcheck N == length(νs)
    xs = Vector{Float64}(N)
    for i in 1:N
        xs[i] = (i == 1) ? (νs[1]*σ*(1 - ρ^2)^(-0.5)) : (ρ*xs[i-1] + σ*νs[i])
    end
    xs + log.(ϵs) + 1.27
end

simulate_stochastic(ρ, σ, N) = simulate_stochastic(ρ, σ, rand(Chisq(1), N), randn(N))

struct Toy_Vol_Problem
    "observed data"
    ys
    "prior for ρ (persistence)"
    prior_ρ
    "prior for σ_v (volatility of volatility)"
    prior_σ
    "χ^2 draws for simulation"
    ϵ::Vector{Float64}
    "Normal(0,1) draws for simulation"
    ν::Vector{Float64}
end

"""
    Toy_Vol_Problem(ys, prior_ρ, prior_σ, M)

Convenience constructor for Toy_Vol_Problem.
Take in the observed data, the priors, and number of simulations (M).
"""
function Toy_Vol_Problem(ys, prior_ρ, prior_σ, M)
    Toy_Vol_Problem(ys, prior_ρ, prior_σ, rand(Chisq(1), M), randn(M))
end

"""
    simulate!(pp::Toy_Vol_Problem)

Updates the shocks of the model.
"""
function simulate!(pp::Toy_Vol_Problem)
    @unpack ϵ, ν = pp
    rand!(Chisq(1), ϵ)
    randn!(ν)
end

"""
    β, v = OLS(y, x)

Take in the dependant variable (y) and the regressor (x), give back the estimated coefficients (β) and the variance (v).
"""
function OLS(y, x)
    β = x \ y
    err = (y - x * β)
    v = mean(abs2, err)
    return β, v
end

## In this form, we use an AR(2) process of the first differences with an intercept as the auxiliary model.

"""
    lag(xs, n, K)

Lag-`n` operator on vector `xs` (maximum `K` lags).
"""
lag(xs, n, K) = xs[((K+1):end)-n]

lag_matrix(xs, ns, K = maximum(ns)) = hcat([lag(xs, n, K) for n in ns]...)

@test lag_matrix(1:5, 1:3) == [3 2 1; 4 3 2]

"first auxiliary regression y, X, meant to capture first differences"
function yX1(zs, K)
    Δs = diff(zs)
    lag(Δs, 0, K), hcat(lag_matrix(Δs, 1:K, K), ones(length(Δs)-K), lag(zs, 1, K+1))
end

"second auxiliary regression y, X, meant to capture levels"
function yX2(zs, K)
    lag(zs, 0, K), hcat(ones(length(zs)-K), lag_matrix(zs, 1:K, K))
end

function logdensity(pp::Toy_Vol_Problem, θ)
    @unpack ys, prior_ρ, prior_σ, ν, ϵ = pp
    ρ, σ = θ
    N = length(ϵ)

    if (abs(ρ) > 1 - eps() || σ ≤ 0 )
        return -Inf
    else

        logprior = logpdf(prior_ρ, ρ) + logpdf(prior_σ, σ)

        # Generating xs, which is the latent volatility process

        zs = simulate_stochastic(ρ, σ, ϵ, ν)
        β₁, v₁ = OLS(yX1(zs, 2)...)
        β₂, v₂ = OLS(yX2(zs, 2)...)

        # We work with first differences
        y₁, X₁ = yX1(ys, 2)
        log_likelihood1 = sum(logpdf.(Normal(0, √v₁), y₁ - X₁ * β₁))
        y₂, X₂ = yX2(ys, 2)
        log_likelihood2 = sum(logpdf.(Normal(0, √v₂), y₂ - X₂ * β₂))
        logprior + log_likelihood1 + log_likelihood2
    end
end

# Trial
ρ = 0.8
σ = 1.2
y = simulate_stochastic(ρ, σ, 10000)
pp = Toy_Vol_Problem(y, Uniform(-1, 1), InverseGamma(1, 1), 10000)

# visualize posterior
ρ_grid = linspace(0.75, 0.9, 50)
σ_grid = linspace(0.9, 1.4, 30)
logdensity_on_grid = ((ρ, σ)->logdensity(pp, [ρ, σ])).(ρ_grid', σ_grid)
logdensity_on_grid .-= maximum(logdensity_on_grid)
heatmap(ρ_grid, σ_grid, exp.(logdensity_on_grid), xlab = "ρ", ylab = "σ",
        title = "posterior")

chain, a = mcmc(RWMH(diagm([0.02, 0.1])), pp, [ρ, σ], 10000) # Σ was hand-tuned
posterior = vcat(chain[5000:end]'...)

# convergence diagnostics (visual)
plot(posterior, label = ["ρ" "σ"])

# marginal posterior densities vs prior
plt = plot(density(posterior[:, 1]), label = "posterior", title = "ρ")
plot!(plt, a -> pdf(pp.prior_ρ, a), linspace(0.6, 1, 100), label = "prior")
vline!(plt, [ρ], label = "true value")

plt = plot(density(posterior[:, 2]), label = "posterior", title = "σ²")
plot!(plt, a -> pdf(pp.prior_σ, a), linspace(0, 5, 100), label = "prior")
vline!(plt, [σ], label = "true value")

@testset "first difference" begin
    ρ̄, σ̄ = vec(mean(posterior, 1))
    @test ρ̄ ≈ ρ rtol = 0.02
    @test σ̄ ≈ σ rtol = 0.02
end


text:
module AssignTaxonomy

using FASTX, BioSequences, LoopVectorization

include("classifier.jl")
# Write your package code here.

end

ext: jl


text:
#=
test_parser:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-02-13
=#
using Test

@testset "parser" begin
    using HyperGraphTools.Parser
    using HyperGraphTools.SimpleTAGMLTokenizer

    r = my_parse()
    @test r == "Hello World"

    function has_no_errors(x)
        isempty(filter(x -> x.type == SimpleTAGMLTokenizer.ERROR, x))
    end

    # use case 1: 2 layers, full overlap
    tagml1 = "[root>[tag|a>[other|b>one two three<other]<tag]<root]"
    grammar1 = """
    S -> ROOT
    ROOT -> [root> TAG <root]
    TAG -> [tag|a> OTHER <tag]
    OTHER -> [other|b> TEXT <other]
    TEXT -> "one two three"
    """
    t1 = tokenize(tagml1)
    @test has_no_errors(t1)

    # use case 2: 2 layers, partial overlap
    tagml2 = "[root>[tag|a>one [other|b>two<tag] three<other]<root]"
    grammar2 = """
    S -> ROOT
    ROOT -> [root> ROOTBODY <root]
    ROOTBODY -> (OPENTAG | CLOSETAG | TEXT)+
    """
    t2 = tokenize(tagml2)
    @test has_no_errors(t2)

    # use case 3: 2 layers, no overlap
    tagml3 = "[root>[tag|a>one<tag] two [other|b>three<other]<root]"
    t3 = tokenize(tagml3)
    @test has_no_errors(t3)

    # use case 4: 2 layers, shared tags
    tagml4a = "[root+a,+b>[tag|a,b>[other|a,b>one two three<other]<tag]<root]"
    tagml4b = "[root+a,+b>[tag|a,b>[other|a,b>[w|a>one<w] [b|b>two<b] three<other]<tag]<root]"
    t4a = tokenize(tagml4a)
    t4b = tokenize(tagml4b)
    @test has_no_errors(t4a)
    @test has_no_errors(t4b)

end

max_stars_repo_name: rhdekker/HyperGraphTools

max_issues_repo_name: rhdekker/HyperGraphTools

max_forks_repo_name: rhdekker/HyperGraphTools


text:
using CodecZstd
using TranscodingStreams
using Test

@testset "Zstd Codec" begin
    codec = ZstdCompressor()
    @test codec isa ZstdCompressor
    @test occursin(r"^ZstdCompressor\(level=\d+\)$", sprint(show, codec))
    @test CodecZstd.initialize(codec) === nothing
    @test CodecZstd.finalize(codec) === nothing

    codec = ZstdDecompressor()
    @test codec isa ZstdDecompressor
    @test occursin(r"^ZstdDecompressor\(\)$", sprint(show, codec))
    @test CodecZstd.initialize(codec) === nothing
    @test CodecZstd.finalize(codec) === nothing

    data = [0x28, 0xb5, 0x2f, 0xfd, 0x04, 0x50, 0x19, 0x00, 0x00, 0x66, 0x6f, 0x6f, 0x3f, 0xba, 0xc4, 0x59]
    @test read(ZstdDecompressorStream(IOBuffer(data))) == b"foo"
    @test read(ZstdDecompressorStream(IOBuffer(vcat(data, data)))) == b"foofoo"

    @test ZstdCompressorStream <: TranscodingStreams.TranscodingStream
    @test ZstdDecompressorStream <: TranscodingStreams.TranscodingStream

    TranscodingStreams.test_roundtrip_read(ZstdCompressorStream, ZstdDecompressorStream)
    TranscodingStreams.test_roundtrip_write(ZstdCompressorStream, ZstdDecompressorStream)
    TranscodingStreams.test_roundtrip_lines(ZstdCompressorStream, ZstdDecompressorStream)
    TranscodingStreams.test_roundtrip_transcode(ZstdCompressor, ZstdDecompressor)
end


text:
# ------------------------------------------------------------------------------
# Read converted NCA5 files and add 15 mm baseline correction to account for the
# 2000 - 2005 difference between both
# ------------------------------------------------------------------------------
module GriddedProjections
using NetCDF
using NCDatasets

function RunGriddedProjections(settings)
    println("\nGridded scenarios...")
    # Locations and time
    fn = settings["dir_NCA5"]*"NCA5_Low_grid.nc"
    lon = ncread(fn,"lon")
    lat = ncread(fn,"lat")
    years_NCA5 = ncread(fn,"years",start=[1],count=[14])

    # Create output NetCDF file
    fh = Dataset(settings["fn_proj_gri"],"c")
    fh.attrib["title"] = "Gridded projections"
    fh.attrib["description"] = "Gridded projections for the interagency report: Global and Regional Sea Level Rise Scenarios for the United States: Updated Mean Projections and Extreme Water Level Probabilities Along U.S. Coastlines"
    fh.attrib["processes"] = "AIS: Antarctic Ice Sheet, GIS: Greenland, glaciers: Glaciers and Ice Caps, landwaterstorage: Liquid water storage changes on land, oceandynamics: Ocean dynamics and global thermosteric expansion, verticallandmotion: Vertical land motion, total: All processes combined."

    defDim(fh,"years", length(years_NCA5))
    defDim(fh,"percentiles",3)
    defDim(fh,"lon", length(lon))
    defDim(fh,"lat", length(lat))

    defVar(fh,"lon",lon,("lon",),deflatelevel=5, attrib = Dict("description" => "Longitude", "units" => "degrees East"))
    defVar(fh,"lat",lat,("lat",),deflatelevel=5, attrib = Dict("description" => "Latitude", "units" => "degrees North"))
    defVar(fh,"years",years_NCA5,("years",),deflatelevel=5 ,attrib = Dict("units" => "years"))
    defVar(fh,"percentiles",convert.(Int32,[17,50,83]),("percentiles",),deflatelevel=5, attrib = Dict("description" => "The percentile of the projection. The 50th percentile shows the median projection, and the 17th and 83rd show the upper- and lower bound on the 1 sigma level.", "units" => "-"))
    # Write scenarios
    for scn in settings["NCA5_scenarios"]
        fn_in = settings["dir_NCA5"]*"NCA5_"*scn*"_grid.nc"
        for prc in settings["processes"]
            NCA5_prc = ncread(fn_in,prc,start=[1,1,1,1],count=[-1,-1,14,-1]);
            if (prc == "total") | (prc == "verticallandmotion")
                NCA5_prc .+= 15
            end
            defVar(fh,"rsl_"*prc*"_"*scn,NCA5_prc,("lon","lat","years","percentiles",),deflatelevel=5, attrib = Dict("description" => "Projected local sea level for scenario "*scn*" and process "*prc*", relative to year 2000.", "units" => "mm"))
        end
    end
    close(fh)
    println("Gridded scenarios done\n")
end

end


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Heap" begin
    heap = DS.MaxHeap()
    @test heap isa DS.Heap
    for i=1:27
        push!(heap, i)
    end
    DS.build!(heap)
    @test sort!(heap) == collect(27:-1:1)

    heap = DS.MinHeap()
    for i=30:-1:1
        push!(heap, i)
    end
    DS.build!(heap)

    @test length(heap) == 30
    @test DS.level(heap, 5) == 3
    @test DS.root(heap) == 1
    @test DS.left_child(heap, 5) == 10
    @test DS.right_child(heap, 5) == 11

    @test repr(heap) == "Heap(\n" *
        "\n\t\t\tHeapNode(16)\n" *
        "\t\t\t\tHeapNode(24)\n" *
        "\t\tHeapNode(3)\n" *
        "\t\t\t\tHeapNode(30)\n" *
        "\t\t\tHeapNode(17)\n" *
        "\t\t\t\tHeapNode(28)\n" *
        "\tHeapNode(2)\n" *
        "\t\t\t\tHeapNode(18)\n" *
        "\t\t\tHeapNode(5)\n" *
        "\t\t\t\tHeapNode(25)\n" *
        "\t\tHeapNode(4)\n" *
        "\t\t\t\tHeapNode(19)\n" *
        "\t\t\tHeapNode(6)\n" *
        "\t\t\t\tHeapNode(7)\n" *
        "HeapNode(1)\n" *
        "\t\t\t\tHeapNode(29)\n" *
        "\t\t\tHeapNode(20)\n" *
        "\t\t\t\tHeapNode(26)\n" *
        "\t\tHeapNode(9)\n" *
        "\t\t\t\tHeapNode(21)\n" *
        "\t\t\tHeapNode(10)\n" *
        "\t\t\t\tHeapNode(11)\n" *
        "\tHeapNode(8)\n" *
        "\t\t\t\tHeapNode(22)\n" *
        "\t\t\tHeapNode(13)\n" *
        "\t\t\t\tHeapNode(27)\n" *
        "\t\tHeapNode(12)\n" *
        "\t\t\t\tHeapNode(23)\n" *
        "\t\t\tHeapNode(14)\n" *
        "\t\t\t\tHeapNode(15)\n" *
    ")"

    @test sort!(heap) == collect(1:30)
end


text:
module Frequency

import LombScargle
import Trapz

#=
This function calculates a lomb scargle transformation
:param n: is the array that contains the NN-inetrvals
:return: the result of the lomb scargle transformation
=#
function lomb_scargle(n)
    t=cumsum(n).-n[1]
    t=t./1000
    plan=LombScargle.plan(t,n,normalization=:psd,minimum_frequency=0.003,maximum_frequency=0.4)
    return LombScargle.lombscargle(plan)
end # lomb_scargle

#=
This function calculates the power of a frequency band between two given frequencys
:param freq: The frequency of a lomb scargle transformation
:param power: The power of a lomb scargle transformation
:param min: The minimum value of the frequency band to be calculated
:param max: The maximum value of the frequency band to be calculated
:return p: The power of the frequency band
=#
function get_power(freq,power,min,max)
    count=1
    index=[]
    for f in freq
        if f>=min && f<max
            push!(index,count)
        end
        count+=1
    end
    p=Trapz.trapz(freq[index[1]:index[end]],power[index[1]:index[end]])
    return p
end # get_power

end # module


text:
"""
    init_model_heisenberg!(
        J :: Vector{Vector{Float64}},
        l :: Lattice
        ) :: Nothing

Init Heisenberg model on a given lattice by overwriting the respective bonds.
Here, J[n] is the coupling to the n-th nearest neighbor (Euclidean norm).
If there are m symmetry inequivalent n-th nearest neighbors, these are
* uniformly initialized if J[n] is a single value
* initialized in ascending bond distance from the origin, if J[n] is an array of length m
"""
function init_model_heisenberg!(
    J :: Vector{Vector{Float64}},
    l :: Lattice
    ) :: Nothing

    # iterate over sites and add Heisenberg couplings to lattice bonds
    for i in eachindex(l.sites)
       for n in eachindex(J)
            # find n-th nearest neighbors
            nbs = get_nbs(n, l.sites[i], l.sites)

            # uniform initialization for n-th nearest neighbor, if no further couplings provided
            if length(J[n]) == 1
                for j in nbs
                    add_bond!(J[n][1], l.bonds[i, j], 1, 1)
                    add_bond!(J[n][1], l.bonds[i, j], 2, 2)
                    add_bond!(J[n][1], l.bonds[i, j], 3, 3)
                end
            # initialize symmetry non-equivalent bonds interactions in ascending bond-length order
            else
                # get bond distances of neighbors
                dist = Int64[get_metric(l.sites[j], l.sites[i], l.uc) for j in nbs]

                # filter out classes of bond distances
                nbkinds = sort(unique(dist))

                # sanity check
                @assert length(J[n]) == length(nbkinds) "$(l.name) has $(length(nbkinds)) inequivalent $(n)-th nearest neighbors, but $(length(J[n])) couplings were supplied."

                for nk in eachindex(nbkinds)
                    # filter out neighbors with dist == nbkinds[nk]
                    nknbs = nbs[findall(x -> x == nbkinds[nk], dist)]

                    for j in nknbs
                        add_bond!(J[n][nk], l.bonds[i, j], 1, 1)
                        add_bond!(J[n][nk], l.bonds[i, j], 2, 2)
                        add_bond!(J[n][nk], l.bonds[i, j], 3, 3)
                    end
                end
            end
        end
    end

    return nothing
end


text:
#!/usr/bin/env julia
# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy
using RainFARM
using ArgParse
using Compat, Compat.Printf

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--radius", "-r"
            help = "Smoothing radius (in grid units)"
            arg_type = Float64
            default = 0.
        "--lapse", "-l"
            help = "Lapse rate [K/km]"
            arg_type =  Float64
            default = 6.5
        "--varname", "-v"
            help = "Input variable name (in orofile)"
            arg_type = AbstractString
            default = ""
        "--coord", "--box", "-b"
            help = "Coordinates of box in which to downscale [lon1 lon2 lat1 lat2]"
 	    nargs = 4
            arg_type =  Float64
	    default =  [ 7, 7.65, 45.35, 45.75 ]
        "--orocoarse", "-o"
            help = "Coarse orography of original data"
            arg_type = AbstractString
            default = ""
        "orofile"
            help = "The input file to use for orography"
            arg_type = AbstractString
            required = true
        "infile"
            help = "The file to downscale"
            arg_type = AbstractString
            required = true
        "outfile"
            help = "The output file name"
            arg_type = AbstractString
            required = true
    end

    s.description="Downscale temperature with a given lapse rate"
    s.version="0.1"
    s.add_version=true

    return parse_args(s)
end

args = parse_commandline()
lapse=abs(args["lapse"])
radius=args["radius"]
filein=args["infile"]
fileout=args["outfile"]
fileoro=args["orofile"]
fileorocoarse=args["orocoarse"]
varname=args["varname"]
coord=args["coord"]

lon1=coord[1]
lon2=coord[2]
lat1=coord[3]
lat2=coord[4]

println("Downscaling in box: ",coord[1],"/",coord[2],"/",coord[3],"/",coord[4])

(tin0,lonl0,latl0)=read_netcdf2d(filein,varname);
if(length(size(lonl0))>1)
dxl=max(lonl0[2,1]-lonl0[1,1],lonl0[1,2]-lonl0[1,1]);
else
dxl=lonl0[2]-lonl0[1];
end

println("dx=",dxl)
if(radius==0)
    radius=dxl/2
end

#Add buffer
lonl1=lon1-radius
lonl2=lon2+radius
latl1=lat1-radius
latl2=lat2+radius
println("box + buffer=",lonl1,"/",lonl2,"/",latl1,"/",latl2)

rr=round.(Int32,rand()*1000000)

run(`cdo -s sellonlatbox,$lonl1,$lonl2,$latl1,$latl2 $fileoro orocut$rr.nc`)
(oro,lonl2,latl2,oroname)=read_netcdf2d("orocut$rr.nc","");
if(length(size(lonl2))>1)
dxf=max(lonl2[2,1]-lonl2[1,1],lonl2[1,2]-lonl2[1,1]);
else
dxf=lonl2[2]-lonl2[1];
end

println("Remapping input data ...")
run(`cdo -s -b F32 remapnn,orocut$rr.nc $filein input_nn$rr.nc`)
(tin,lonl,latl,varname)=read_netcdf2d("input_nn$rr.nc",varname);
oro=float(oro); # convert to float
(nx,ny,nt)=size(tin)

nf2=div(radius,abs(dxf))
println("Smoothing radius = ",radius," = ",nf2, " pixel")

println("Preparing correction ...")
if(fileorocoarse=="")
   oros=smooth(oro,nf2)
else
   run(`cdo -s -b F32 remapnn,orocut$rr.nc $fileorocoarse orocut_coarse$rr.nc`)
   (oroc,lonl,latl,orocname)=read_netcdf2d("orocut_coarse$rr.nc","");
   oros=smooth(oroc,nf2)
   run(`rm orocut_coarse$rr.nc`)
end
#println("oro=",mean(oro)," oros=",mean(oros))
oro=-(oro-oros)*lapse/1000.

write_netcdf2d("orocorr_temp$rr.nc",oro,lonl,latl,oroname,"orocut$rr.nc")
run(`cdo -s sellonlatbox,$lon1,$lon2,$lat1,$lat2 orocorr_temp$rr.nc orocorr$rr.nc `)

println("Downscaling ...")
for i=1:nt
    println("t=",i)
    tins=smooth(tin[:,:,i],nf2)
    tin[:,:,i]=tins+oro
#println("tin=",mean(tin[:,:,i])," tins=",mean(tins))
end

write_netcdf2d("out_temp$rr.nc",tin,lonl,latl,varname,"input_nn$rr.nc")
run(`cdo -s sellonlatbox,$lon1,$lon2,$lat1,$lat2 out_temp$rr.nc $fileout `)
run(`rm orocut$rr.nc input_nn$rr.nc out_temp$rr.nc orocorr_temp$rr.nc orocorr$rr.nc`)


text:
module SocialNet

export PriorityQueue, peek, decrease_priority!

export Graph, adj_list, add_vertex!, add_edge!, del_edge!,
    ne, nv, shortest_path_length, vertices

export cycle_graph, complete_graph, star_graph, star_cycle_graph, random_graph

export Network, person, persons, add_person!, friend_pair!, defriend_pair!, friends,
    degree_of_separation

include("priority_queue.jl")
include("graph.jl")
include("network.jl")
include("graph_library.jl")

end # module SocialNet


text:
atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch
    end
end

using MyJulia


text:
export Wishart

import Distributions: Wishart
import Base: ndims

vague(::Type{ <: Wishart }, dims::Int) = Wishart(dims, Matrix(Diagonal(huge .* ones(dims))))

Base.ndims(dist::Wishart) = first(size(dist))

prod_analytical_rule(::Type{ <: Wishart }, ::Type{ <: Wishart }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::Wishart, right::Wishart)
    d = dim(left)

    ldf, lS = params(left)
    rdf, rS = params(right)

    V  = (lS * cholinv(lS + rS) * rS) |> Hermitian |> Matrix
    df = ldf + rdf - d - 1

    return Wishart(df, V)
end


text:
import Libdl

const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depsjl_path)
    # User-provided SCIP library
    include(depsjl_path)
else
    # Artifact from BinaryBuilder package
    import SCIP_PaPILO_jll
    if SCIP_PaPILO_jll.is_available()
        using SCIP_PaPILO_jll: libscip
    else
        using SCIP_jll: libscip
    end
end

function __init__()
    major = SCIPmajorVersion()
    minor = SCIPminorVersion()
    patch = SCIPtechVersion()
    current = VersionNumber("$major.$minor.$patch")
    required = VersionNumber("8")
    upperbound = VersionNumber("9")
    if current < required || current >= upperbound
        @error("SCIP is installed at version $current, " *
              "supported are $required up to (excluding) $upperbound.")
    end
end


text:
# general array for testing input to SpecialFunctions.jl
x32 = [-1.0; 0.0; 1e-6; 0.01:0.01:60.0]
x64 = [-1.0; 0.0; 1e-6; 0.01:0.01:600.0]

### Tests for besseli0
i032_SpecialFunctions = SpecialFunctions.besseli.(0, x32)
i064_SpecialFunctions = SpecialFunctions.besseli.(0, x64)

i0_64 = besseli0.(Float64.(x64))
i0_32 = besseli0.(Float32.(x32))

# make sure output types match input types
@test i0_64[1] isa Float64
@test i0_32[1] isa Float32

# test against SpecialFunctions.jl
@test i0_64 ≈ i064_SpecialFunctions
@test i0_32 ≈ i032_SpecialFunctions

### Tests for besseli0x
i0x32_SpecialFunctions = SpecialFunctions.besselix.(0, x32)
i0x64_SpecialFunctions = SpecialFunctions.besselix.(0, x64)

i0x_64 = besseli0x.(Float64.(x64))
i0x_32 = besseli0x.(Float32.(x32))

# make sure output types match input types
@test i0x_64[1] isa Float64
@test i0x_32[1] isa Float32

# test against SpecialFunctions.jl
@test i0x_64 ≈ i0x64_SpecialFunctions
@test i0x_32 ≈ i0x32_SpecialFunctions

### Tests for besseli1
i132_SpecialFunctions = SpecialFunctions.besseli.(1, x32)
i164_SpecialFunctions = SpecialFunctions.besseli.(1, x64)

i1_64 = besseli1.(Float64.(x64))
i1_32 = besseli1.(Float32.(x32))

# make sure output types match input types
@test i1_64[1] isa Float64
@test i1_32[1] isa Float32

# test against SpecialFunctions.jl
@test i1_64 ≈ i164_SpecialFunctions
@test i1_32 ≈ i132_SpecialFunctions

### Tests for besseli1x
i1x32_SpecialFunctions = SpecialFunctions.besselix.(1, x32)
i1x64_SpecialFunctions = SpecialFunctions.besselix.(1, x64)

i1x_64 = besseli1x.(Float64.(x64))
i1x_32 = besseli1x.(Float32.(x32))

# make sure output types match input types
@test i1x_64[1] isa Float64
@test i1x_32[1] isa Float32

# test against SpecialFunctions.jl
@test i1x_64 ≈ i1x64_SpecialFunctions
@test i1x_32 ≈ i1x32_SpecialFunctions

# test for besseli
# test small arguments and order
m = 0:1:200; x = 0.5f0:0.5f0:90.0f0
@test besseli(10, 1.0f0) isa Float32
@test besseli(2, 80.0f0) isa Float32
@test besseli(112, 80.0f0) isa Float32

for m in m, x in x
    @test besseli(m, x) ≈ Float32(SpecialFunctions.besseli(m, x))
end

#Float 64
m = 0:1:200; x = 0.1:0.5:150.0
@test besseli(10, 1.0) isa Float64
@test besseli(2, 80.0) isa Float64
@test besseli(112, 80.0) isa Float64
t = [besseli(m, x) for m in m, x in x]

@test t[10] isa Float64
@test t ≈ [SpecialFunctions.besseli(m, x) for m in m, x in x]

t = [besselix(m, x) for m in m, x in x]
@test t[10] isa Float64
@test t ≈ [SpecialFunctions.besselix(m, x) for m in m, x in x]


text:
function music(cur_player::AbstractPlayer)

  cur_js = """

    if (typeof tetrisMusic !== 'undefined') {

      if ( tetrisMusic.playing() ) {
        tetrisMusic.pause();
        \$(".js-music-button").closest("a").removeClass("cs-active");
      } else {
        tetrisMusic.play();
        \$(".js-music-button").closest("a").addClass("cs-active");
      }

    }

  """

  tetris_js(
    cur_player.game.scope,
    JSString(cur_js)
  )

end


text:
@testset "EuclideanVector"  begin

#---

using ACE
using Random, Printf, Test, LinearAlgebra, ACE.Testing
using ACE: evaluate, evaluate_d, SymmetricBasis, NaiveTotalDegree, PIBasis
using ACE.Random: rand_rot, rand_refl

# construct the 1p-basis
D = NaiveTotalDegree()
maxdeg = 6
ord = 3

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, D = D)

# generate a configuration
nX = 10
Xs = rand(EuclideanVectorState, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

#---

@info("SymmetricBasis construction and evaluation: EuclideanVector")

φ = ACE.EuclideanVector(Complex{Float64})
pibasis = PIBasis(B1p, ord, maxdeg; property = φ, isreal=false)
basis = SymmetricBasis(pibasis, φ)

BB = evaluate(basis, cfg)

# a stupid but necessary test
BB1 = basis.A2Bmap * evaluate(basis.pibasis, cfg)
println(@test isapprox(BB, BB1, rtol=1e-10))

##

@info("Test equivariance properties")

tol = 1E-9

@info("check for rotation, permutation and inversion equivariance")
for ntest = 1:20
   Xs = rand(EuclideanVectorState, B1p.bases[1], nX)
   BB = evaluate(basis, ACEConfig(Xs))
   Q = rand([-1,1]) * ACE.Random.rand_rot()
   Xs_rot = Ref(Q) .* shuffle(Xs)
   BB_rot = evaluate(basis, ACEConfig(Xs_rot))
   print_tf(@test all([ norm(Q' * b1 - b2) < tol
                        for (b1, b2) in zip(BB_rot, BB)  ]))
end
println()

##

@info(" ... derivatives")
for ntest = 1:30
   Us = randn(SVector{3, Float64}, length(Xs))
   C = randn(typeof(φ.val), length(basis))
   F = t -> sum( sum(c .* b.val)
                 for (c, b) in zip(C, ACE.evaluate(basis, ACEConfig(Xs + t[1] * Us))) )
   dF = t -> [ sum( sum(c .* db)
                    for (c, db) in zip(C, ACE.evaluate_d(basis, ACEConfig(Xs + t[1] * Us)) * Us) ) ]
   print_tf(@test fdtest(F, dF, [0.0], verbose=false))
end
println()

##

end


text: Dict{String,Any}("negpi" => Dict{String,Any}("value" => "-3.14","type" => "float"),"pospi" => Dict{String,Any}("value" => "3.14","type" => "float"),"pi" => Dict{String,Any}("value" => "3.14","type" => "float"),"zero-intpart" => Dict{String,Any}("value" => "0.123","type" => "float"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
struct CircularReference
   objref::UInt64
end

mutable struct SharedObject
   obj::Any
   refcount::UInt64
end

function SharedObject(obj)
   SharedObject(obj, UInt32(1))
end

# Encapsulate the io stream and a lock to prevent concurrent writing
# by multiple Tasks because this would mess up the messages
# Making the struct mutable makes it easier to share the object with R,
# as it always has a fixed address
mutable struct CommunicatoR{T}
   # for synchronizing the communication
   lock::ReentrantLock

   # the iostream for communicating with R
   io::T

   # used to detect circular references
   objrefs::Set{UInt64}

   # for sharing references between Julia and R
   sharedheap::Dict{UInt64, SharedObject}

   # whether the next translation is forced to be a full translation
   # in contrast to a mere passing of a reference
   full_translation::Bool

   # Contains identifiers of finalized callback functions.
   # These can be communicated to R to be removed there as well.
   finalized_callbacks::Vector{String}

   # contains callback functions and their corresponding identifiers
   # for the communication with R
   registered_callbacks::Dict{Function, String}
end

function CommunicatoR(io)
   CommunicatoR(
      ReentrantLock(),
      io,
      Set{UInt64}(),
      Dict{UInt64, SharedObject}(),
      false,
      Vector{String}(),
      Dict{Function, String}()
   )
end

function write_bin(c::CommunicatoR, x)
   write(c.io, x)
end

function read_bin(c::CommunicatoR, n::Int)
   ret = read(c.io, n)
   while length(ret) < n
      ret = [ret; read(c.io, length(ret) - n)]
   end
   ret
end

function read_byte(c::CommunicatoR)
   read(c.io, UInt8)
end

# GLOBAL variable nclients to count the number of clients that are connected
const nclients = Ref{Int}(0)
const nclients_lock = ReentrantLock()

function inc_nclients()
   lock(nclients_lock)
   nclients.x += 1
   unlock(nclients_lock)
end

# decrement number of clients. Stop server if it is reduced to 0.
function dec_nclients()
   lock(nclients_lock)
   nclients.x -= 1
   if nclients.x <= 0
      exit(0)
   end
   unlock(nclients_lock)
end

"""
    serve(port_hint; keeprunnning = false, portfile = "")
Starts the JuliaConnectoR server, which will accept connections.

It is attempted to start the server at the port given via `port_hint`.
The actual port may be different. If the argument `portfile` is specified,
the real port is written to the file with the given name
as decimal number in form of a string.

If the argument `multiclient` is `false`,
the server accepts only one client and stops after the client disconnects.
If the argument `multiclient` is set to `true`,
the server will accept multiple clients and keeps running
after a client disconnects.
Note that multiple clients share the same global workspace.
"""
function serve(port_hint::Int;
      multiclient::Bool = false,
      portfile::String = "")

   realport, server = listenany(port_hint)
   println("Julia is listening on port $realport.")
   # write real port to portfile
   if (portfile != "")
      open(portfile, "w") do f
         write(f, string(Int(realport)))
         write(f, '\n')
      end
   end

   # Prevents crash when run as script and interrupted
   ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)

   if multiclient
      while true
         sock = accept(server)
         @async serve_repl(sock)
         inc_nclients()
      end
   else
      sock = accept(server)
      close(server) # don't listen for other clients
      serve_repl(sock)
   end
end

function serve_repl(sock)
   communicator = CommunicatoR(sock)
   capture_outputstreams(communicator)

   while isopen(sock)
      call = Call()
      try
         firstbyte = read_bin(communicator, 1)
         if length(firstbyte) == 1 && firstbyte[1] == CALL_INDICATOR
            # Parse incoming function call
            call = read_call(communicator)
         else
            if !(length(firstbyte) == 0 || firstbyte[1] == BYEBYE)
               close(sock)
               error("Unexpected leading character $(firstbyte[1])")
            else
               close(sock)
               dec_nclients()
               return
            end
         end
      catch ex
         @warn "Unexpected parsing error. Stopping listening. " *
               "Original error: $ex"
         dec_nclients()
         return
      end

      result = evaluate_checked!(call, communicator)
      write_message(communicator, result)
   end

   dec_nclients()
end

function callbackfun(callbackid::String, communicator::CommunicatoR)

   callbackfinalizer = CallbackFinalizer(callbackid, communicator)

   ret = (args...; kwargs...) -> begin
      stayalivewithme(callbackfinalizer)

      posargs = isempty(args) ? Vector{Any}() : collect(Any, args)
      callbackargs = ElementList(posargs, Dict{Symbol, Any}(kwargs))
      write_callback_message(communicator, callbackid, callbackargs)

      # read, parse and return answer
      while true
         firstbyte = read_bin(communicator, 1)[1]
         if firstbyte == RESULT_INDICATOR
            answer = read_element(communicator)
            parsingcheck(answer)
            return evaluate!(answer, communicator)
         elseif firstbyte == CALL_INDICATOR
            call = read_call(communicator)
            result = evaluate_checked!(call, communicator)
            write_message(communicator, result)
         elseif firstbyte == FAIL_INDICATOR
            r_errormsg = read_string(communicator)
            error("Error in R callback: " * r_errormsg)
         else
            error("Unexpected leading byte: " * string(firstbyte))
         end
      end
   end

   communicator.registered_callbacks[ret] = callbackid
   ret
end

function callanonymous(communicator::CommunicatoR, functionref::Vector{UInt8},
      args...; kwargs...)

   objref::ImmutableObjectReference =
         communicator.sharedheap[parseheapref(functionref)].obj
   f::Function = objref.obj
   ret = f(args...; kwargs...)
   ret
end

function start_result_message(c::CommunicatoR)
   # Ensure that output is transferred before the result:
   # Let the tasks redirecting the outputstreams take over
   # (see function redirect_outputstream)
   flush(c.io)
   flush(stdout)
   yield()
   flush(stderr)
   yield()

   # after the output is handled,
   # take control of the outputstream to send the result
   lock(c.lock)
end

start_fail_message = start_result_message
start_callback_message = start_result_message

function end_result_message(c::CommunicatoR)
   # end of circular reference detection
   empty!(c.objrefs)

   # release lock of outputstream for next call
   unlock(c.lock)
end

end_fail_message = end_result_message
end_callback_message = end_result_message

function start_output_message(c::CommunicatoR)
   lock(c.lock)
end

function end_output_message(c::CommunicatoR)
   unlock(c.lock)
end

function capture_outputstreams(communicator::CommunicatoR)
   redirect_outputstream(communicator, redirect_stdout()[1], STDOUT_INDICATOR)
   rd, wr = redirect_stderr()

   # Also redirect the logger.
   global_logger(SimpleLogger(wr))
   # For debugging, comment out the line above, use @debug lines and
   # set the environment variable "JULIA_DEBUG" to "all", e. g. in PowerShell
   # $env:JULIA_DEBUG = "all".
   # Start Julia and execute:
   # `import Tables; include("inst/Julia/RConnector.jl"); RConnector.serve(11980)`
   # Set the environment variable JULIACONNECTOR_SERVER=localhost:11980
   # for the R session. Then start R and load the JuliaConnectoR.
   # If you execute commands, debugging info is displayed on the command line.

   redirect_outputstream(communicator, rd, STDERR_INDICATOR)
end

# Starts an asynchronous task that sends the output to R
function redirect_outputstream(communicator::CommunicatoR, reader,
      output_type_indicator::UInt8)

   @async while isopen(communicator.io) && isopen(reader)
      iobuffer = readavailable(reader)
      write_output_message(communicator, iobuffer, output_type_indicator)
   end
end

function write_output_message(communicator::CommunicatoR, msg::Vector{UInt8},
      output_type_indicator::UInt8)

   start_output_message(communicator)
   write_bin(communicator, output_type_indicator)
   write_bin(communicator, Int32(length(msg)))
   write_bin(communicator, msg)
   end_output_message(communicator)
end

const DEFAULT_DISPLAY_LINES = 24
const DEFAULT_DISPLAY_COLUMNS = 80

function showobj(x, width::Int)
   io = IOBuffer()
   displaysize = (DEFAULT_DISPLAY_LINES, width)
   context = IOContext(io, :limit => true, :displaysize => displaysize)
   show(context, "text/plain", x)
   String(take!(io))
end

function showobj(x, width)
   try
      intwidth = convert(Int, width)
   catch ex
      intwidth = DEFAULT_DISPLAY_COLUMNS
   end
   showobj(x, width)
end


text:
macro assert_equal(a, b)
    quote
        aval = $a
        bval = $b
        if aval != bval
            println("Failure: (", aval, " != ", bval, "), expected it to ==")
        else
            print(".")
        end
    end
end

macro assert_not_equal(a, b)
    quote
        aval = $a
        bval = $b
        if aval == bval
            println("Failure: (", aval, " == ", bval, "), expected it to !=")
        else
            print(".")
        end
    end
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8]);F = [6724*x1^2*x3^2+2916*x1^2*x5^2+9072*x1^2*x5*x6+7056*x1^2*x6^2+1600*x1^2*x8^2-1440*x1*x2^2*x8+13448*x1*x2*x3*x4-5724*x1*x2*x4*x5-8904*x1*x2*x4*x6-5120*x1*x2*x5*x8-4800*x1*x2*x6*x8+11316*x1*x3^2*x4+7020*x1*x3^2*x5+10920*x1*x3^2*x6-972*x1*x3*x5*x6-492*x1*x3*x5*x8-12828*x1*x3*x6^2-6400*x1*x8^3+324*x2^4+2304*x2^3*x5+2160*x2^3*x6+9533*x2^2*x4^2+4096*x2^2*x5^2+7680*x2^2*x5*x6+3616*x2^2*x6^2+296*x2^2*x6*x8+4249*x2^2*x8^2-6890*x2*x3^2*x4+11316*x2*x3*x4^2+954*x2*x3*x4*x6-492*x2*x4*x5*x8-11316*x2*x4*x6^2+272*x2*x5*x6*x8+12756*x2*x5*x8^2+9600*x2*x6*x8^2+4225*x3^4-1170*x3^3*x6+4761*x3^2*x4^2+81*x3^2*x6^2-414*x3*x4*x5*x8-9522*x3*x4*x6^2+1165*x5^2*x8^2+414*x5*x6^2*x8+4761*x6^4+6400*x8^4-152*x1*x2*x6-1406*x1*x2*x8+7668*x1*x5*x7+5108*x1*x5*x8+11928*x1*x6*x7-2880*x2^2*x5+256*x2^2*x6+2368*x2^2*x8-632*x2*x4*x6-7526*x2*x4*x7-5846*x2*x4*x8-10240*x2*x5^2-9600*x2*x5*x6+2176*x2*x5*x8+9230*x3^2*x7-1278*x3*x6*x7-5372*x4*x5*x8-12800*x5*x8^2+361*x1^2-1216*x1*x2+6560*x1*x3+3002*x1*x4+1024*x2^2+1504*x2*x4+5520*x3*x4+6241*x4^2+6400*x5^2-240*x5*x8-5520*x6^2+5041*x7^2+1600, -2880*x1*x2*x8+13448*x1*x3*x4-5724*x1*x4*x5-8904*x1*x4*x6-5120*x1*x5*x8-4800*x1*x6*x8+1296*x2^3+6912*x2^2*x5+6480*x2^2*x6+19066*x2*x4^2+8192*x2*x5^2+15360*x2*x5*x6+7232*x2*x6^2+592*x2*x6*x8+8498*x2*x8^2-6890*x3^2*x4+11316*x3*x4^2+954*x3*x4*x6-492*x4*x5*x8-11316*x4*x6^2+272*x5*x6*x8+12756*x5*x8^2+9600*x6*x8^2-152*x1*x6-1406*x1*x8-5760*x2*x5+512*x2*x6+4736*x2*x8-632*x4*x6-7526*x4*x7-5846*x4*x8-10240*x5^2-9600*x5*x6+2176*x5*x8-1216*x1+2048*x2+1504*x4, 13448*x1^2*x3+13448*x1*x2*x4+22632*x1*x3*x4+14040*x1*x3*x5+21840*x1*x3*x6-972*x1*x5*x6-492*x1*x5*x8-12828*x1*x6^2-13780*x2*x3*x4+11316*x2*x4^2+954*x2*x4*x6+16900*x3^3-3510*x3^2*x6+9522*x3*x4^2+162*x3*x6^2-414*x4*x5*x8-9522*x4*x6^2+18460*x3*x7-1278*x6*x7+6560*x1+5520*x4, 13448*x1*x2*x3-5724*x1*x2*x5-8904*x1*x2*x6+11316*x1*x3^2+19066*x2^2*x4-6890*x2*x3^2+22632*x2*x3*x4+954*x2*x3*x6-492*x2*x5*x8-11316*x2*x6^2+9522*x3^2*x4-414*x3*x5*x8-9522*x3*x6^2-632*x2*x6-7526*x2*x7-5846*x2*x8-5372*x5*x8+3002*x1+1504*x2+5520*x3+12482*x4, 5832*x1^2*x5+9072*x1^2*x6-5724*x1*x2*x4-5120*x1*x2*x8+7020*x1*x3^2-972*x1*x3*x6-492*x1*x3*x8+2304*x2^3+8192*x2^2*x5+7680*x2^2*x6-492*x2*x4*x8+272*x2*x6*x8+12756*x2*x8^2-414*x3*x4*x8+2330*x5*x8^2+414*x6^2*x8+7668*x1*x7+5108*x1*x8-2880*x2^2-20480*x2*x5-9600*x2*x6+2176*x2*x8-5372*x4*x8-12800*x8^2+12800*x5-240*x8, 9072*x1^2*x5+14112*x1^2*x6-8904*x1*x2*x4-4800*x1*x2*x8+10920*x1*x3^2-972*x1*x3*x5-25656*x1*x3*x6+2160*x2^3+7680*x2^2*x5+7232*x2^2*x6+296*x2^2*x8+954*x2*x3*x4-22632*x2*x4*x6+272*x2*x5*x8+9600*x2*x8^2-1170*x3^3+162*x3^2*x6-19044*x3*x4*x6+828*x5*x6*x8+19044*x6^3-152*x1*x2+11928*x1*x7+256*x2^2-632*x2*x4-9600*x2*x5-1278*x3*x7-11040*x6, 7668*x1*x5+11928*x1*x6-7526*x2*x4+9230*x3^2-1278*x3*x6+10082*x7, 3200*x1^2*x8-1440*x1*x2^2-5120*x1*x2*x5-4800*x1*x2*x6-492*x1*x3*x5-19200*x1*x8^2+296*x2^2*x6+8498*x2^2*x8-492*x2*x4*x5+272*x2*x5*x6+25512*x2*x5*x8+19200*x2*x6*x8-414*x3*x4*x5+2330*x5^2*x8+414*x5*x6^2+25600*x8^3-1406*x1*x2+5108*x1*x5+2368*x2^2-5846*x2*x4+2176*x2*x5-5372*x4*x5-25600*x5*x8-240*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
using SplineFEEC
using Documenter

makedocs(;
    modules=[SplineFEEC],
    authors="JuliaDEC",
    repo="https://github.com/JuliaDEC/SplineFEEC.jl/blob/{commit}{path}#L{line}",
    sitename="SplineFEEC.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaDEC.github.io/SplineFEEC.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Basis" => "bases.md",
        "Operator" => "operator.md",
        "Examples" => "examples.md",
        "Library" => "library.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaDEC/SplineFEEC.jl",
    devbranch = "documentation",
)


text:
using CausalDiscovery.Model: generateTree, proposeTree, findNodeWithPosition, getExpr, getConditionalLogProb, getPriorProb, getPriorLogProb, isValid

# using CausalDiscovery.Grammer:

#tree1
# TaggedParseTree(
# NonTerminalNode((1, 2), "expr", (1, 2), 0.5, nothing,
#     Any[NonTerminalNode((2, 1), "line", (1, 2, 2, 1), 0.5,parent, 1/2
#         Any[NonTerminalNode((3, 2), "exo_line", (1, 2, 2, 1, 3, 2), 0.5, parent, 1/2
#             Any[TerminalNode((17, 5), "float_var", (1, 2, 2, 1, 3, 2, 17, 5), 0.2, parent, 1/2
#                    :float_var_5), 1/5
#             NonTerminalNode((21, 1), "float_dist", (1, 2, 2, 1, 3, 2, 21, 1), 0.5, parent,
#                 Any[TerminalNode((19, 1), "normal_params", (1, 2, 2, 1, 3, 2, 21, 1, 19, 1), 0.5, parent, 1/2
#                        (0, 1))]) 1/2
#             ])
#         ])
#     ])
@testset "MCMC" begin
    tree1 = generateTree(100)
    proposed = proposeTree(tree1)

    #tree1 should equal proposed expect with float_var_5 replacing float_var_1
    #following test should be added when expressions are able to be evaluated
    #@test(isapprox(eval(getExpr(tree1)), eval(getExpr(proposed)); 0.01)

    #Tests findNodeWithPosition function
    expected_node = tree1.root_node.children[1].children[1]
    found_node = findNodeWithPosition(tree1.root_node, expected_node.node_position)
    @test(expected_node == found_node)

    #Testing random seed version of generateTree
    tree2 = generateTree(100)
    @test(tree1.node_positions == tree2.node_positions)

    expr1 = getExpr(tree1)
    expr2 = getExpr(tree2)
    @test(expr1==expr2)

    #test getPriorProb with exogenous line
    priorProb = getPriorProb(tree1)
    expectedProb = 0.00625
    @test(isapprox(priorProb, expectedProb; .01))

    #test getPriorProb with proposed tree
    priorProb = getPriorProb(proposed)
    @test(isapprox(priorProb, expectedProb; .01))

    #test getPriorLogProb
    priorLogProb = getPriorLogProb(proposed)
    expected = log(expectedProb)
    @test(isapprox(priorLogProb, expected; 0.01))

    #testing conditional
    expectedConditional = log((1/6)*(1/5) + (1/6)*.00625 + (1/30)*(1/8) + (1/30)*(1/16))
    conditional = getConditionalLogProb(tree1, proposed)
    @test(isapprox(conditional, expectedConditional; .01))

    tree3 = generateTree(4)

    #Test getPriorProb with endogenous line
    priorProb = getPriorProb(tree3)
    expectedProb = 1/(2*2*3*4*5*5)
    @test(isapprox(priorProb, expectedProb; 0.01))

    proposed = proposeTree(tree3)
    #Tests for conditional prob with endo line
    expectedConditional = log((1/6)*(1/25)*((1/48)+(1/24)+(1/12)+(1/4)))
    conditional = getConditionalLogProb(tree3, proposed)
    @test(isapprox(conditional, expectedConditional; 0.01))

    new_tree = generateTree(3)
    #bool_var_4 = bool_var_5
    #bool_var_5 ~ Bernoulli(0.1)
    @test(isValid(new_tree, (a = "a",)) == false)

    new_tree = generateTree(5)

    @test(isValid(new_tree, (bool_var_4 = 1,)) == true)
    @test(isValid(new_tree, (bool_var_7 = 1,)) == false)
end


text:
# Force the random number generator to be the same for all runs of the script
using Random
Random.seed!(777)

include("global_curved.jl")

import PGFPlots

let
  τs = 1

  verts = [-1 0 1 -1 0 1;
            0 0 0  1 1 1]
  EToV = [1 2 4 5;
          2 3 5 6]'
  EToF = [1 2 3 4;
          2 5 6 7]'
  FToB = fill(BC_DIRICHLET, (maximum(EToF),))
  FToB[2] = BC_LOCKED_INTERFACE

  (FToE, FToLF, EToO, EToS) = connectivityarrays(EToV, EToF)
  # plot_connectivity(verts, EToV)

  (nelems, nfaces) = (size(EToV, 2), size(FToB, 1))

  # Size of the grid
  SBP_orders = (2, 4, 6)

  num_samp = 1000
  min_eig_noSchur = zeros(length(SBP_orders), num_samp)
  max_eig_noSchur = zeros(length(SBP_orders), num_samp)
  min_eig_Schur_D = zeros(length(SBP_orders), num_samp)
  max_eig_Schur_D = zeros(length(SBP_orders), num_samp)
  min_eig_Schur_M = zeros(length(SBP_orders), num_samp)
  max_eig_Schur_M = zeros(length(SBP_orders), num_samp)
  for k = 1:num_samp
    k % 10 == 1 && println("sample $k of $num_samp")
    for (i, SBPp) in enumerate(SBP_orders)
      Ns = Nr = ntuple(e -> 3*SBPp - 1, nelems)

      # Build the local operators
      lop = ntuple(nelems) do e
        metrics = create_metrics(SBPp, Nr[e], Ns[e])

        # Build coefficients which are positive definite
        λ1 = rand(Nr[e]+1, Ns[e]+1)
        λ2 = rand(Nr[e]+1, Ns[e]+1) / 10000

        q  = π * rand(Nr[e]+1, Ns[e]+1)
        @. metrics.crr = λ1 * cos(q)^2 + λ2 * sin(q)^2
        @. metrics.css = λ1 * sin(q)^2 + λ2 * cos(q)^2
        @. metrics.crs = (λ2 - λ1) * cos(q) * sin(q)

        # Build local operators
        locoperator(SBPp, Nr[e], Ns[e], metrics; τscale = τs)
      end

      #
      # Assemble the global volume operators
      #
      (M, FbarT, D, vstarts, FToλstarts) =
      LocalGlobalOperators(lop, Nr, Ns, FToB, FToE, FToLF, EToO, EToS,
                           (x) -> cholesky(Symmetric(x)))
      locfactors = M.F

      # Get a unique array indexes for the face to jumps map
      FToδstarts = bcstarts(FToB, FToE, FToLF, BC_JUMP_INTERFACE, Nr, Ns)

      # Compute the number of volume, trace (λ), and jump (δ) points
      VNp = vstarts[nelems+1]-1
      λNp = FToλstarts[nfaces+1]-1
      δNp = FToδstarts[nfaces+1]-1

      Np = (Nr.+1) .* (Ns.+1)
      A11 = [lop[1].M̃ spzeros(Np[1], Np[2])
             spzeros(Np[2], Np[1]) lop[2].M̃]
      A = [A11 FbarT';
           FbarT Diagonal(D)]

      # @show (VNp, λNp, δNp)
      # Build the (sparse) λ matrix using the schur complement and factor
      B = assembleλmatrix(FToλstarts, vstarts, EToF, FToB, locfactors, D, FbarT)
      @assert B ≈ Diagonal(D) - FbarT * (A11 \ Matrix(FbarT'))
      C = A11 - FbarT' * Diagonal(1 ./ D) * FbarT

      evA = eigen(Matrix(A)).values
      evB = eigen(Matrix(B)).values
      evC = eigen(Matrix(C)).values

      min_eig_noSchur[i, k], max_eig_noSchur[i, k] = extrema(evA)
      min_eig_Schur_M[i, k], max_eig_Schur_M[i, k] = extrema(evB)
      min_eig_Schur_D[i, k], max_eig_Schur_D[i, k] = extrema(evC)
    end
  end

  # Plot No Schur eigenvalues
  Lx = (1, num_samp)
  Ly = extrema(min_eig_noSchur)
  Ly = (floor(10000 * Ly[1]) / 10000, ceil(10000 * Ly[2]) / 10000)
  plt = Plot(BrailleCanvas(40, 20,
                           origin_x = Lx[1], origin_y = Ly[1],
                           width = Lx[2] - Lx[1], height = Ly[2] - Ly[1]))

  annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
  annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
  annotate!(plt, :bl, string(Lx[1]), color = :light_black)
  annotate!(plt, :br, string(Lx[2]), color = :light_black)
  for (i, SBPp) = enumerate(SBP_orders)
    lineplot!(plt, min_eig_noSchur[i, :])
  end
  display(plt)

  pgf_axis = PGFPlots.Axis(style="width=5cm, height=5cm", xlabel="realization",
                           ylabel=PGFPlots.L"$\min \lambda$", xmin = 1,
                           xmax = num_samp)
  for (i, SBPp) = enumerate(SBP_orders)
    push!(pgf_axis, PGFPlots.Linear(1:num_samp, min_eig_noSchur[i, :],
                                    style="no marks, thick",
                                    legendentry = "order = $SBPp"))
  end
  PGFPlots.save("eigenvalues_noSchur.tikz", pgf_axis)

  # Plot Schur of D block eigenvalues
  Lx = (1, num_samp)
  Ly = extrema(min_eig_Schur_D)
  Ly = (floor(10000 * Ly[1]) / 10000, ceil(10000 * Ly[2]) / 10000)
  plt = Plot(BrailleCanvas(40, 20,
                           origin_x = Lx[1], origin_y = Ly[1],
                           width = Lx[2] - Lx[1], height = Ly[2] - Ly[1]))

  annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
  annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
  annotate!(plt, :bl, string(Lx[1]), color = :light_black)
  annotate!(plt, :br, string(Lx[2]), color = :light_black)
  for (i, SBPp) = enumerate(SBP_orders)
    lineplot!(plt, min_eig_Schur_D[i, :])
  end
  display(plt)

  pgf_axis = PGFPlots.Axis(style="width=5cm, height=5cm", xlabel="realization",
                           ylabel=PGFPlots.L"$\min \lambda$", xmin = 1,
                           xmax = num_samp)
  for (i, SBPp) = enumerate(SBP_orders)
    push!(pgf_axis, PGFPlots.Linear(1:num_samp, min_eig_Schur_D[i, :],
                                    style="no marks, thick",
                                    legendentry = "order = $SBPp"))
  end
  PGFPlots.save("eigenvalues_Schur_D.tikz", pgf_axis)

  # Plot Schur of M block eigenvalues
  Lx = (1, num_samp)
  Ly = extrema(min_eig_Schur_M)
  Ly = (floor(10000 * Ly[1]) / 10000, ceil(10000 * Ly[2]) / 10000)
  plt = Plot(BrailleCanvas(40, 20,
                           origin_x = Lx[1], origin_y = Ly[1],
                           width = Lx[2] - Lx[1], height = Ly[2] - Ly[1]))

  annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
  annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
  annotate!(plt, :bl, string(Lx[1]), color = :light_black)
  annotate!(plt, :br, string(Lx[2]), color = :light_black)
  for (i, SBPp) = enumerate(SBP_orders)
    lineplot!(plt, min_eig_Schur_M[i, :])
  end
  display(plt)

  pgf_axis = PGFPlots.Axis(style="width=5cm, height=5cm", xlabel="realization",
                           ylabel=PGFPlots.L"$\min \lambda$", xmin = 1,
                           xmax = num_samp)
  for (i, SBPp) = enumerate(SBP_orders)
    push!(pgf_axis, PGFPlots.Linear(1:num_samp, min_eig_Schur_M[i, :],
                                    style="no marks, thick",
                                    legendentry = "order = $SBPp"))
  end
  PGFPlots.save("eigenvalues_Schur_M.tikz", pgf_axis)
end


text:
export get_metadata

function get_metadata(search_path::String; include_parents=true)
    search_path = abspath(search_path)
    while search_path != ""
        search_path == projectdir() && return nothing
        id = find_file_in_index(search_path)
        if id !== nothing
            return Metadata(search_path)
        end
        include_parents || return nothing
        search_path, _ = splitdir(search_path)
    end
    return nothing
end

function get_metadata(f::Function)
    ms = Metadata[]
    for file in filter(x->endswith(x,".bson"),readdir(metadatadir()))
        m = load_metadata(joinpath(metadatadir(),file), ignore_exceptions=true)
        m === nothing && continue
        f(m) && push!(ms, m)
    end
    return ms
end

function get_metadata()
    return get_metadata() do m
        true
    end
end

function get_metadata(field::String, value)
    return get_metadata() do m
        field in keys(m) && m[field] == value
    end
end

max_stars_repo_name: sebastianpech/DrWatsonSim

max_issues_repo_name: sebastianpech/DrWatsonSim

max_forks_repo_name: sebastianpech/DrWatsonSim


text:
macro c(ret_type, func, arg_types, lib)
    local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
    quote
        $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, lib)) ),
                                            $ret_type, $arg_types, $(args_in...) )
    end
end

const unix = 1
const linux = 1
const EOF = (-1)
# Skipping MacroDefinition: NULL((void*)0)
const BUFSIZ =
const SEEK_SET = 0
const SEEK_CUR = 1
const SEEK_END = 2
const P_tmpdir = "/tmp"
const L_tmpnam = 20
const TMP_MAX = 238328
const FILENAME_MAX = 4096
const L_ctermid = 9
const FOPEN_MAX = 16
const stdin =
const stdout =
const stderr =
# Skipping MacroDefinition: getc(_fp)_IO_getc(_fp)
# Skipping MacroDefinition: putc(_ch,_fp)_IO_putc(_ch,_fp)
# Skipping MacroDefinition: InvalidOid((Oid)0)
const OID_MAX =
const PG_DIAG_SEVERITY = 'S'
const PG_DIAG_SQLSTATE = 'C'
const PG_DIAG_MESSAGE_PRIMARY = 'M'
const PG_DIAG_MESSAGE_DETAIL = 'D'
const PG_DIAG_MESSAGE_HINT = 'H'
const PG_DIAG_STATEMENT_POSITION = 'P'
const PG_DIAG_INTERNAL_POSITION = 'p'
const PG_DIAG_INTERNAL_QUERY = 'q'
const PG_DIAG_CONTEXT = 'W'
const PG_DIAG_SOURCE_FILE = 'F'
const PG_DIAG_SOURCE_LINE = 'L'
const PG_DIAG_SOURCE_FUNCTION = 'R'
const PG_COPYRES_ATTRS = 0x01
const PG_COPYRES_TUPLES = 0x02
const PG_COPYRES_EVENTS = 0x04
const PG_COPYRES_NOTICEHOOKS = 0x08
# manually added pg formats
const PGF_TEXT = 0
const PGF_BINARY = 1
# Skipping MacroDefinition: PQsetdb(M_PGHOST,M_PGPORT,M_PGOPT,M_PGTTY,M_DBNAME)PQsetdbLogin(M_PGHOST,M_PGPORT,M_PGOPT,M_PGTTY,M_DBNAME,NULL,NULL)
# Skipping MacroDefinition: PQfreeNotify(ptr)PQfreemem(ptr)
const PQnoPasswordSupplied = "fe_sendauth: no password supplied\n"
typealias _IO_lock_t Union{}
typealias va_list Cint
# typealias off_t __off_t
# typealias ssize_t __ssize_t
# typealias fpos_t _G_fpos_t
typealias Oid UInt32
# begin enum ConnStatusType
typealias ConnStatusType UInt32
const CONNECTION_OK = 0
const CONNECTION_BAD = 1
const CONNECTION_STARTED = 2
const CONNECTION_MADE = 3
const CONNECTION_AWAITING_RESPONSE = 4
const CONNECTION_AUTH_OK = 5
const CONNECTION_SETENV = 6
const CONNECTION_SSL_STARTUP = 7
const CONNECTION_NEEDED = 8
# end enum ConnStatusType
# begin enum PostgresPollingStatusType
typealias PostgresPollingStatusType UInt32
const PGRES_POLLING_FAILED = 0
const PGRES_POLLING_READING = 1
const PGRES_POLLING_WRITING = 2
const PGRES_POLLING_OK = 3
const PGRES_POLLING_ACTIVE = 4
# end enum PostgresPollingStatusType
# begin enum ExecStatusType
typealias ExecStatusType UInt32
const PGRES_EMPTY_QUERY = 0
const PGRES_COMMAND_OK = 1
const PGRES_TUPLES_OK = 2
const PGRES_COPY_OUT = 3
const PGRES_COPY_IN = 4
const PGRES_BAD_RESPONSE = 5
const PGRES_NONFATAL_ERROR = 6
const PGRES_FATAL_ERROR = 7
const PGRES_COPY_BOTH = 8
const PGRES_SINGLE_TUPLE = 9
# end enum ExecStatusType
# begin enum PGTransactionStatusType
typealias PGTransactionStatusType UInt32
const PQTRANS_IDLE = 0
const PQTRANS_ACTIVE = 1
const PQTRANS_INTRANS = 2
const PQTRANS_INERROR = 3
const PQTRANS_UNKNOWN = 4
# end enum PGTransactionStatusType
# begin enum PGVerbosity
typealias PGVerbosity UInt32
const PQERRORS_TERSE = 0
const PQERRORS_DEFAULT = 1
const PQERRORS_VERBOSE = 2
# end enum PGVerbosity
# begin enum PGPing
typealias PGPing UInt32
const PQPING_OK = 0
const PQPING_REJECT = 1
const PQPING_NO_RESPONSE = 2
const PQPING_NO_ATTEMPT = 3
# end enum PGPing
typealias PQnoticeReceiver Ptr{Void}
typealias PQnoticeProcessor Ptr{Void}
typealias pqbool UInt8
typealias pgthreadlock_t Ptr{Void}
typealias PGconn Void
typealias PGresult Void

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
function _ensure_not_running_on_functor(f, name)
    # if x itself is a Type, then it is a constructor, thus not a functor.
    # This also catchs UnionAll constructors which have a `:var` and `:body` fields
    f isa Type && return

    if fieldcount(typeof(f)) > 0
        throw(ArgumentError(
            "$name cannot be used on closures/functors (such as $f)"
        ))
    end
end

"""
    _wrap_function(f, xs, ignores)

Return a new version of `f`, `fnew`, that ignores some of the arguments `xs`.

# Arguments
- `f`: The function to be wrapped.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === nothing`.
"""
function _wrap_function(f, xs, ignores)
    function fnew(sigargs...)
        callargs = Any[]
        j = 1

        for (i, (x, ignore)) in enumerate(zip(xs, ignores))
            if ignore
                push!(callargs, x)
            else
                push!(callargs, sigargs[j])
                j += 1
            end
        end
        @assert j == length(sigargs) + 1
        @assert length(callargs) == length(xs)
        return f(callargs...)
    end
    return fnew
end

"""
    _make_j′vp_call(fdm, f, ȳ, xs, ignores) -> Tuple

Call `FiniteDifferences.j′vp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `ȳ`: The adjoint w.r.t. output of `f`.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === nothing`.

# Returns
- `∂xs::Tuple`: Derivatives estimated by finite differencing.
"""
function _make_j′vp_call(fdm, f, ȳ, xs, ignores)
    f2 = _wrap_function(f, xs, ignores)

    ignores = collect(ignores)
    args = Any[nothing for _ in 1:length(xs)]
    all(ignores) && return (args...,)
    sigargs = xs[.!ignores]
    arginds = (1:length(xs))[.!ignores]
    fd = j′vp(fdm, f2, ȳ, sigargs...)
    @assert length(fd) == length(arginds)

    for (dx, ind) in zip(fd, arginds)
        args[ind] = _maybe_fix_to_composite(dx)
    end
    return (args...,)
end

"""
    _make_jvp_call(fdm, f, xs, ẋs, ignores)

Call `FiniteDifferences.jvp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ẋs`: The directional derivatives of `xs` w.r.t. some real number `t`.
- `ignores`: Collection of `Bool`s, the same length as `xs` and `ẋs`.
   If `ignores[i] === true`, then `ẋs[i]` is ignored for derivative estimation.

# Returns
- `Ω̇`: Derivative of output w.r.t. `t` estimated by finite differencing.
"""
function _make_jvp_call(fdm, f, xs, ẋs, ignores)
    f2 = _wrap_function(f, xs, ignores)

    ignores = collect(ignores)
    all(ignores) && return ntuple(_->nothing, length(xs))
    sigargs = zip(xs[.!ignores], ẋs[.!ignores])
    return _maybe_fix_to_composite(jvp(fdm, f2, sigargs...))
end

# TODO: remove after https://github.com/JuliaDiff/FiniteDifferences.jl/issues/97
# For functions which return a tuple, FD returns a tuple to represent the differential. Tuple
# is not a natural differential, because it doesn't overload +, so make it a Composite.
_maybe_fix_to_composite(x::Tuple) = Composite{typeof(x)}(x...)
_maybe_fix_to_composite(x::NamedTuple) = Composite{typeof(x)}(;x...)
_maybe_fix_to_composite(x) = x

"""
    test_scalar(f, z; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), check_inferred=true, kwargs...)

Given a function `f` with scalar input and scalar output, perform finite differencing checks,
at input point `z` to confirm that there are correct `frule` and `rrule`s provided.

# Arguments
- `f`: Function for which the `frule` and `rrule` should be tested.
- `z`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).

`fkwargs` are passed to `f` as keyword arguments.
If `check_inferred=true`, then the type-stability of the `frule` and `rrule` are checked.
All remaining keyword arguments are passed to `isapprox`.
"""
function test_scalar(f, z; rtol=1e-9, atol=1e-9, fdm=_fdm, fkwargs=NamedTuple(), check_inferred=true, kwargs...)
    # To simplify some of the calls we make later lets group the kwargs for reuse
    rule_test_kwargs = (; rtol=rtol, atol=atol, fdm=fdm, fkwargs=fkwargs, check_inferred=check_inferred, kwargs...)
    isapprox_kwargs = (; rtol=rtol, atol=atol, kwargs...)

    _ensure_not_running_on_functor(f, "test_scalar")
    # z = x + im * y
    # Ω = u(x, y) + im * v(x, y)
    Ω = f(z; fkwargs...)

    # test jacobian using forward mode
    Δx = one(z)
    @testset "$f at $z, with tangent $Δx" begin
        # check ∂u_∂x and (if Ω is complex) ∂v_∂x via forward mode
        frule_test(f, (z, Δx); rule_test_kwargs...)
        if z isa Complex
            # check that same tangent is produced for tangent 1.0 and 1.0 + 0.0im
            _, real_tangent = frule((Zero(), real(Δx)), f, z; fkwargs...)
            _, embedded_tangent = frule((Zero(), Δx), f, z; fkwargs...)
            check_equal(real_tangent, embedded_tangent; isapprox_kwargs...)
        end
    end
    if z isa Complex
        Δy = one(z) * im
        @testset "$f at $z, with tangent $Δy" begin
            # check ∂u_∂y and (if Ω is complex) ∂v_∂y via forward mode
            frule_test(f, (z, Δy); rule_test_kwargs...)
        end
    end

    # test jacobian transpose using reverse mode
    Δu = one(Ω)
    @testset "$f at $z, with cotangent $Δu" begin
        # check ∂u_∂x and (if z is complex) ∂u_∂y via reverse mode
        rrule_test(f, Δu, (z, Δx); rule_test_kwargs...)
        if Ω isa Complex
            # check that same cotangent is produced for cotangent 1.0 and 1.0 + 0.0im
            _, back = rrule(f, z)
            _, real_cotangent = back(real(Δu))
            _, embedded_cotangent = back(Δu)
            check_equal(real_cotangent, embedded_cotangent; isapprox_kwargs...)
        end
    end
    if Ω isa Complex
        Δv = one(Ω) * im
        @testset "$f at $z, with cotangent $Δv" begin
            # check ∂v_∂x and (if z is complex) ∂v_∂y via reverse mode
            rrule_test(f, Δv, (z, Δx); rule_test_kwargs...)
        end
    end
end

"""
    _test_inferred(f, args...; kwargs...)

Simple wrapper for `@inferred f(args...: kwargs...)`, avoiding the type-instability in not
knowing how many `kwargs` there are.
"""
function _test_inferred(f, args...; kwargs...)
    if isempty(kwargs)
        @inferred f(args...)
    else
        @inferred f(args...; kwargs...)
    end
end

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), check_inferred=true, kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

`fkwargs` are passed to `f` as keyword arguments.
If `check_inferred=true`, then the type-stability of the `frule` is checked.
All remaining keyword arguments are passed to `isapprox`.
"""
function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, fkwargs=NamedTuple(), check_inferred=true, kwargs...)
    # To simplify some of the calls we make later lets group the kwargs for reuse
    isapprox_kwargs = (; rtol=rtol, atol=atol, kwargs...)

    _ensure_not_running_on_functor(f, "frule_test")

    xs = first.(xẋs)
    ẋs = last.(xẋs)
    check_inferred && _test_inferred(frule, (NO_FIELDS, deepcopy(ẋs)...), f, deepcopy(xs)...; deepcopy(fkwargs)...)
    Ω_ad, dΩ_ad = frule((NO_FIELDS, deepcopy(ẋs)...), f, deepcopy(xs)...; deepcopy(fkwargs)...)
    Ω = f(deepcopy(xs)...; deepcopy(fkwargs)...)
    check_equal(Ω_ad, Ω; isapprox_kwargs...)

    ẋs_is_ignored = ẋs .== nothing
    # Correctness testing via finite differencing.
    dΩ_fd = _make_jvp_call(fdm, (xs...) -> f(deepcopy(xs)...; deepcopy(fkwargs)...), xs, ẋs, ẋs_is_ignored)
    check_equal(dΩ_ad, dΩ_fd; isapprox_kwargs...)

    # No tangent is passed in to test accumlation, so generate one
    # See: https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/66
    acc = rand_tangent(Ω)
    _check_add!!_behavour(acc, dΩ_ad; rtol=rtol, atol=atol, kwargs...)
end

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), check_inferred=true, kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
  Should be same structure as `f(x)` (so if multiple returns should be a tuple)
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

`fkwargs` are passed to `f` as keyword arguments.
If `check_inferred=true`, then the type-stability of the `rrule` and the pullback it
returns are checked.
All remaining keyword arguments are passed to `isapprox`.
"""
function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, check_inferred=true, fkwargs=NamedTuple(), kwargs...)
    # To simplify some of the calls we make later lets group the kwargs for reuse
    isapprox_kwargs = (; rtol=rtol, atol=atol, kwargs...)

    _ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    xs = first.(xx̄s)
    accumulated_x̄ = last.(xx̄s)
    check_inferred && _test_inferred(rrule, f, xs...; fkwargs...)
    y_ad, pullback = rrule(f, xs...; fkwargs...)
    y = f(xs...; fkwargs...)
    check_equal(y_ad, y; isapprox_kwargs...)  # make sure primal is correct

    check_inferred && _test_inferred(pullback, ȳ)
    ∂s = pullback(ȳ)
    ∂self = ∂s[1]
    x̄s_ad = ∂s[2:end]
    @test ∂self === NO_FIELDS  # No internal fields

    # Correctness testing via finite differencing.
    x̄s_is_dne = accumulated_x̄ .== nothing
    x̄s_fd = _make_j′vp_call(fdm, (xs...) -> f(xs...; fkwargs...), ȳ, xs, x̄s_is_dne)
    for (accumulated_x̄, x̄_ad, x̄_fd) in zip(accumulated_x̄, x̄s_ad, x̄s_fd)
        if accumulated_x̄ === nothing  # then we marked this argument as not differentiable
            @assert x̄_fd === nothing  # this is how `_make_j′vp_call` works
            @test x̄_ad isa DoesNotExist  # we said it wasn't differentiable.
        else
            x̄_ad isa AbstractThunk && check_inferred && _test_inferred(unthunk, x̄_ad)

            # The main test of the actual deriviative being correct:
            check_equal(x̄_ad, x̄_fd; isapprox_kwargs...)
            _check_add!!_behavour(accumulated_x̄, x̄_ad; isapprox_kwargs...)
        end
    end

    check_thunking_is_appropriate(x̄s_ad)
end

function check_thunking_is_appropriate(x̄s)
    @testset "Don't thunk only non_zero argument" begin
        num_zeros = count(x->x isa AbstractZero, x̄s)
        num_thunks = count(x->x isa Thunk, x̄s)
        if num_zeros + num_thunks == length(x̄s)
            @test num_thunks !== 1
        end
    end
end


text:
function part1(input)
    rules, medicine = parse_input(input)
    return length(transform(rules, medicine))
end

function part2(input)
    rules, medicine = parse_input(input)
    rules = reverse.(rules)
    immediate_rules = filter(p -> occursin("Ca", first(p)), rules)
    filter!(p -> !occursin("Ca", first(p)), rules)
    append!(immediate_rules, filter(p -> occursin("TiTi", first(p)), rules))
    filter!(p -> !occursin("TiTi", first(p)), rules)
    RnAr = filter(r -> endswith(r, "Ar"), first.(rules))
    RnAr = unique([r[first(findfirst("Rn", r)):end] for r in RnAr])
    i = 1
    molecules = Dict(medicine => 0)
    queue = [medicine]
    depth = 0
    while !isempty(queue)
        molecule = popfirst!(queue)
        if molecules[molecule] > depth
            depth = molecules[molecule]
        end
        new_molecules = transform_one(immediate_rules, molecule)
        if isempty(new_molecules)
            new_molecules = transform_reduced(rules, RnAr, molecule)
        end
        for new_molecule in new_molecules
            new_molecule == "e" && return molecules[molecule] + 1
            if !haskey(molecules, new_molecule)
                molecules[new_molecule] = molecules[molecule] + 1
                push!(queue, new_molecule)
            end
        end
    end
end

function parse_input(input)
    rules = Tuple{String, String}[]
    for line in eachline(input)
        if occursin(" => ", line)
            push!(rules, Tuple(string.(split(line, " => "))))
        elseif !isempty(line)
            return rules, line
        end
    end
end

function transform(rules, molecule)
    products = Set{String}()
    for (a, b) in rules
        n = length(a) - 1
        for i = 1:(length(molecule) - n)
            if molecule[i:(i + n)] == a
                push!(products, join([molecule[1:(i - 1)], b, molecule[(i + n + 1):end]]))
            end
        end
    end
    return products
end

function transform_one(rules, molecule)
    for (a, b) in rules
        n = length(a) - 1
        for i = 1:(length(molecule) - n)
            if molecule[i:(i + n)] == a
                return [join([molecule[1:(i - 1)], b, molecule[(i + n + 1):end]])]
            end
        end
    end
    return String[]
end

function transform_reduced(rules, RnAr, molecule)
    products = String[]
    if !occursin("Ar", molecule)
        start = 1
        stop = length(molecule)
    else
        stop = last(findfirst("Ar", molecule))
        Rn = findall("Rn", molecule)
        starts = filter(<(stop), first.(Rn))
        start = starts[end]
        if molecule[start:stop] ∈ RnAr
            if length(starts) < 2
                start = 1
            else
                start = starts[end - 1]
            end
        end
    end
    for (a, b) in rules
        n = length(a) - 1
        for i = start:(stop - n)
            if molecule[i:(i + n)] == a
                push!(products, join([molecule[1:(i - 1)], b, molecule[(i + n + 1):end]]))
            end
        end
    end
    return products
end


text:
using DashBootstrapComponents

left_jumbotron = dbc_col(
    html_div(
        [
            html_h2("Change the background", className = "display-3"),
            html_hr(className = "my-2"),
            html_p(
                "Swap the background-color utility and add a `.text-*` color " *
                "utility to mix up the look.",
            ),
            dbc_button("Example Button", color = "light", outline = true),
        ],
        className = "h-100 p-5 text-white bg-dark rounded-3",
    ),
    md = 6,
)

right_jumbotron = dbc_col(
    html_div(
        [
            html_h2("Add borders", className = "display-3"),
            html_hr(className = "my-2"),
            html_p(
                "Or, keep it light and add a border for some added definition " *
                "to the boundaries of your content.",
            ),
            dbc_button("Example Button", color = "secondary", outline = true),
        ],
        className = "h-100 p-5 bg-light border rounded-3",
    ),
    md = 6,
)

jumbotron = dbc_row([left_jumbotron, right_jumbotron], className = "align-items-md-stretch")

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Revise
using BNPRegressionGGA2021
using CSV
using DataStructures
using DataFrames
using Distributions
using Gadfly
using Random
using Statistics
const BNP = BNPRegressionGGA2021

# Define the DGP
function simulate_sample(N0, N1)
    dy1 = LogNormal(0.3, 0.7);
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (0.3, 1.2)], [0.6, 0.4]);
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    X1 = kron([ones(2, 5) [1, 0]], ones(N1))
    y0 = [X0[i, end] == 1 ? rand(dy2) : rand(dy1) for i in 1:N0]
    c0 = rand(Exponential(10), N0)
    y1 = LinRange(0, 6, N1) |> collect |> x -> repeat(x, size(X1, 1) ÷ N1)
    z0 = deepcopy(y0)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

N0, N1, Niter = 500, 2, 1;
df = zeros(Bool, Niter, 6);
Random.seed!(1);
dy1, dy2, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
mean(y0 .> c0)
grid = LinRange(0, 6, 50) |> collect;
plot(
    layer(x = grid, y = 1 .- cdf.(dy1, grid), Geom.line()),
    layer(x = grid, y = 1 .- cdf.(dy2, grid), Geom.line())
)
plot(
    layer(x = grid, y = pdf.(dy1, grid) ./ (1 .- cdf.(dy1, grid)), Geom.line()),
    layer(x = grid, y = pdf.(dy2, grid) ./ (1 .- cdf.(dy2, grid)), Geom.line())
)
m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
_, _, _, chaing = BNP.sample!(m; mcmcsize = 10000, burnin = 5000);
mean(chaing)

# Get the MAP of γ for 100 simulated samples
begin
    N0, N1, Niter = 500, 2, 100;
    df = zeros(Bool, Niter, 6);
    Random.seed!(1);
    for iter in 1:Niter
        _, _, _, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
        m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
        _, _, _, chaing = BNP.sample!(m; mcmcsize = 10000, burnin = 5000);
        df[iter, :] .= begin
            chaing |>
            DataStructures.counter |>
            collect |>
            x -> sort(x, by = x -> x[2], rev = true) |>
            x -> getindex(x, 1) |>
            x -> getindex(x, 1)
        end
        println(iter)
        println(df[iter, :])
    end
    CSV.write("data/simulation-erlang-nph-gamma-bnp.csv", DataFrame(df, :auto))
end

# Save the 100 simulated samples
N0, N1, Niter = 500, 2, 100;
Random.seed!(1);
df = map(1:Niter) do iter
    _, _, _, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
    df = DataFrame(X0, :auto)
    df[!, :iter] .= iter
    df[!, :y0] = y0
    df[!, :c0] = c0
    df
end;
df = reduce(vcat, df);
CSV.write("data/simulation-erlang-nph-data.csv", df);


text: using Base.Test


text:
abstract type AbstractSemiconductor{T} <: AbstractObject{T} end

"""
    struct Semiconductor{T,G,MT,CDM,IDM} <: AbstractSemiconductor{T}

Semiconductor bulk of a [`SolidStateDetector`](@ref).

This is the volume in which electrons and holes will drift during the signal development.

## Parametric types
* `T`: Precision type.
* `G`: Type of `geometry`.
* `MT`: Type of `material`.
* `CDM`: Type of `charge_drift_model`.
* `IDM`: Type of `impurity_density_model`.

## Fields
* `temperature::T`: Temperature (in K) of the semiconductor.
* `material::MT`: Material of the semiconductor.
* `impurity_density_model::IDM`: Impurity density model for the points inside the semiconductor.
* `charge_drift_model::CDM`: Model that describes the drift of electrons and holes inside the semiconductor.
* `geometry::G`: Geometry of the semiconductor, see [Constructive Solid Geometry (CSG)](@ref).

## Definition in Configuration File

A `Semiconductor` is defined through the `semiconductor` field of a detector.
It needs `material` and `geometry`, and can optionally be given a `temperature`, `impurity_density` and `charge_drift_model`.

An example definition of a semiconductor looks like this:
```yaml
semiconductor:
  material: HPGe
  temperature: 78
  impurity_density: # ...
  charge_drift_model: # ...
  geometry: # ...
```
"""
struct Semiconductor{T,G,MT,CDM,IDM} <: AbstractSemiconductor{T}
    temperature::T
    material::MT
    impurity_density_model::IDM
    charge_drift_model::CDM
    geometry::G
end

function Semiconductor{T}(dict::Dict, input_units::NamedTuple, outer_transformations) where T <: SSDFloat
    impurity_density_model = if haskey(dict, "impurity_density")
        ImpurityDensity(T, dict["impurity_density"], input_units)
    elseif haskey(dict, "charge_density_model")
        @warn "Config file deprication: The field \"charge_density_model\" under semiconductor is deprecated.
            It should be changed to \"impurity_density\". In later version this will result in an error.
            For now, it will be treated as an impurity density."
        ImpurityDensity(T, dict["charge_density_model"], input_units)
    else
        ConstantImpurityDensity{T}(0)
    end
    charge_drift_model = if haskey(dict, "charge_drift_model") && haskey(dict["charge_drift_model"], "model")
        cdm = getfield(Main, Symbol(dict["charge_drift_model"]["model"])){T}
        cdm <: AbstractChargeDriftModel{T} ? cdm() : ElectricFieldChargeDriftModel{T}()
    else
        ElectricFieldChargeDriftModel{T}()
    end
    material = material_properties[materials[dict["material"]]]
    temperature = if haskey(dict, "temperature")
        T(dict["temperature"])
    elseif material.name == "High Purity Germanium"
        T(78)
    else
        T(293)
    end
    inner_transformations = parse_CSG_transformation(T, dict, input_units)
    transformations = combine_transformations(inner_transformations, outer_transformations)
    geometry = Geometry(T, dict["geometry"], input_units, transformations)
    return Semiconductor(temperature, material, impurity_density_model, charge_drift_model, geometry)
end

function println(io::IO, d::Semiconductor{T}) where {T <: SSDFloat}
    println("\t---General Properties---")
    println("\t-Detector Material: \t $(d.material.name)")
    println("\t-Charge Drift Model:\t $(typeof(d.charge_drift_model).name.name)")
end

print(io::IO, d::Semiconductor{T}) where {T} = print(io, "Semiconductor{$T} - $(d.material.name)")

show(io::IO, d::Semiconductor) = print(io, d)
show(io::IO,::MIME"text/plain", d::Semiconductor) = show(io, d)

function Semiconductor(sc::Semiconductor{T,G,MT,CDM,IDM}, impurity_density::AbstractImpurityDensity{T}) where {T,G,MT,CDM,IDM}
    Semiconductor(sc.temperature, sc.material, impurity_density, sc.charge_drift_model, sc.geometry)
end
function Semiconductor(sc::Semiconductor{T,G,MT,CDM,IDM}, chargedriftmodel::AbstractChargeDriftModel{T}) where {T,G,MT,CDM,IDM}
    Semiconductor(sc.temperature, sc.material, sc.impurity_density_model, chargedriftmodel, sc.geometry)
end

"""
    scaling_factor_for_permittivity_in_undepleted_region(sc::Semiconductor{T})::T where {T}

This function is called in the calculations of weighting potentials of undepleted detectors.
The electric permittivity, ``ϵ_{r}``, is scaled with this function in areas where the detector is undepleted.
A value between `[0, +Inf]` should be returned. However, `Inf` should not be returned but instead
a very high value should be returned in order to mimic perfect conductivity if that is desired.

## Arguments
* `sc::Semiconductor{T}`: Semiconductor for which the dielectric permittivity should be scaled up.

!!! danger "Experimental feature!"
    This feature is under research. The goal is to study the properties / signal response of undepleted detector.
    This function is indented to be overwritten by the user to study the response.
"""
function scaling_factor_for_permittivity_in_undepleted_region(sc::Semiconductor{T})::T where {T}
    100000
end


text:
include("../emi.jl")
using EMI
using EMI.Draw, EMI.Gmsh

"""Make nxm grid of shapes which are separated by given spacing"""
function fill_canvas(shape::Shape, counts::Tuple{Int, Int}, spacing::Tuple{Real, Real})
    # Compute the transform
    bb = BoundingBox(shape)
    dx = (first(bb.ur) - first(bb.ll)) + first(spacing)
    dy = (last(bb.ur) - last(bb.ll)) + last(spacing)
    shift_x = Point(dx, 0)
    shift_y = Point(0, dy)

    n, m = counts
    model = shape
    shapes = [shape]
    for j in 1:m
        # Fill horizontal o -> o -> o
        for i in 1:n-1
            push!(shapes, model + i*shift_x)
        end
        # Fill vertical
        model = model + shift_y
        j < m && push!(shapes, model)
    end
    Canvas(shapes)
end

# shape = ClosedPolygon([Point(0, 0), Point(1, 0), Point(1, 1)])

const R, dx, dy = 1, 0.2, 0.3
shape = Loop([Line(Point(0, R+dy+dx), Point(0, R+dy)),
              Line(Point(0, R+dy), Point(dy, R+dy)),
              CircleArc(Point(dy, R+dy), Point(R+dy, R+dy), Point(R+dy, dy)),
              Line(Point(R+dy, dy), Point(R+dy, 0)),
              Line(Point(R+dy, 0), Point(R+dy+dx, 0)),
              Line(Point(R+dy+dx, 0), Point(R+dy+dx, dy)),
              CircleArc(Point(R+dy+dx, dy), Point(R+dy+dx, R+dy), Point(2R+dy+dx, R+dy)),
              Line(Point(2R+dy+dx, R+dy), Point(2R+2*dy+dx, R+dy)),
              Line(Point(2R+2*dy+dx, R+dy), Point(2R+2*dy+dx, R+dy+dx)),
              Line(Point(2R+2*dy+dx, R+dy+dx), Point(2R+dy+dx, R+dy+dx)),
              CircleArc(Point(2R+dy+dx, R+dy+dx), Point(R+dy+dx, R+dy+dx), Point(R+dy+dx, 2R+dy+dx)),
              Line(Point(R+dy+dx, 2R+dy+dx), Point(R+dy+dx, 2R+2*dy+dx)),
              Line(Point(R+dy+dx, 2R+2*dy+dx), Point(R+dy, 2R+2*dy+dx)),
              Line(Point(R+dy, 2R+2*dy+dx), Point(R+dy, 2R+dy+dx)),
              CircleArc(Point(R+dy, 2R+dy+dx), Point(R+dy, R+dy+dx), Point(dy, R+dy+dx)),
              Line(Point(dy, R+dy+dx), Point(0, R+dy+dx))])

canvas = fill_canvas(shape, (3, 5), (0.2, 0.1))

set_bbox!(canvas, 0.2, 0.3)
gmsh_script(canvas, 0.2)

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
# t.jl
# Various forms of t-tests
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export OneSampleTTest, TwoSampleTTest, EqualVarianceTTest,
    UnequalVarianceTTest

abstract type TTest <: HypothesisTest end
abstract type TwoSampleTTest <: TTest end

pvalue(x::TTest; tail=:both) = pvalue(TDist(x.df), x.t; tail=tail)

default_tail(test::TTest) = :both

# confidence interval by inversion
function StatsBase.confint(x::TTest; level::Float64=0.95, tail=:both)
    check_level(level)

    if tail == :left
        (-Inf, StatsBase.confint(x, level=1-(1-level)*2)[2])
    elseif tail == :right
        (StatsBase.confint(x, level=1-(1-level)*2)[1], Inf)
    elseif tail == :both
        q = quantile(TDist(x.df), 1-(1-level)/2)
        (x.xbar-q*x.stderr, x.xbar+q*x.stderr)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

## ONE SAMPLE T-TEST

struct OneSampleTTest <: TTest
    n::Int       # number of observations
    xbar::Real   # estimated mean
    df::Int      # degrees of freedom
    stderr::Real # empirical standard error
    t::Real      # t-statistic
    μ0::Real     # mean under h_0
end

testname(::OneSampleTTest) = "One sample t-test"
population_param_of_interest(x::OneSampleTTest) = ("Mean", x.μ0, x.xbar) # parameter of interest: name, value under h0, point estimate

function show_params(io::IO, x::OneSampleTTest, ident="")
    println(io, ident, "number of observations:   $(x.n)")
    println(io, ident, "t-statistic:              $(x.t)")
    println(io, ident, "degrees of freedom:       $(x.df)")
    println(io, ident, "empirical standard error: $(x.stderr)")
end

"""
    OneSampleTTest(xbar::Real, stddev::Real, n::Int, μ0::Real = 0)

Perform a one sample t-test of the null hypothesis that `n` values with mean `xbar` and
sample standard deviation `stddev`  come from a distribution with mean `μ0` against the
alternative hypothesis that the distribution does not have mean `μ0`.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function OneSampleTTest(xbar::Real, stddev::Real, n::Int, μ0::Real=0)
    stderr = stddev/sqrt(n)
    t = (xbar-μ0)/stderr
    df = n-1
    OneSampleTTest(n, xbar, df, stderr, t, μ0)
end

"""
    OneSampleTTest(v::AbstractVector{T<:Real}, μ0::Real = 0)

Perform a one sample t-test of the null hypothesis that the data in vector `v` comes from
a distribution with mean `μ0` against the alternative hypothesis that the distribution
does not have mean `μ0`.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
OneSampleTTest(v::AbstractVector{T}, μ0::Real=0) where {T<:Real} = OneSampleTTest(mean(v), std(v), length(v), μ0)

"""
    OneSampleTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real}, μ0::Real = 0)

Perform a paired sample t-test of the null hypothesis that the differences between pairs of
values in vectors `x` and `y` come from a distribution with mean `μ0` against the
alternative hypothesis that the distribution does not have mean `μ0`.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function OneSampleTTest(x::AbstractVector{T}, y::AbstractVector{S}, μ0::Real=0) where {T<:Real, S<:Real}
    check_same_length(x, y)

    OneSampleTTest(x - y, μ0)
end

## TWO SAMPLE T-TEST (EQUAL VARIANCE)

struct EqualVarianceTTest <: TwoSampleTTest
    n_x::Int     # number of observations
    n_y::Int     # number of observations
    xbar::Real   # estimated mean difference
    df::Int      # degrees of freedom
    stderr::Real # empirical standard error
    t::Real      # t-statistic
    μ0::Real     # mean difference under h_0
end

function show_params(io::IO, x::TwoSampleTTest, ident="")
    println(io, ident, "number of observations:   [$(x.n_x),$(x.n_y)]")
    println(io, ident, "t-statistic:              $(x.t)")
    println(io, ident, "degrees of freedom:       $(x.df)")
    println(io, ident, "empirical standard error: $(x.stderr)")
end

testname(::EqualVarianceTTest) = "Two sample t-test (equal variance)"
population_param_of_interest(x::TwoSampleTTest) = ("Mean difference", x.μ0, x.xbar) # parameter of interest: name, value under h0, point estimate

"""
    EqualVarianceTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real})

Perform a two-sample t-test of the null hypothesis that `x` and `y` come from distributions
with equal means and variances against the alternative hypothesis that the distributions
have different means but equal variances.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function EqualVarianceTTest(x::AbstractVector{T}, y::AbstractVector{S}, μ0::Real=0) where {T<:Real,S<:Real}
    nx, ny = length(x), length(y)
    xbar = mean(x) - mean(y)
    stddev = sqrt(((nx - 1) * var(x) + (ny - 1) * var(y)) / (nx + ny - 2))
    stderr = stddev * sqrt(1/nx + 1/ny)
    t = (xbar - μ0) / stderr
    df = nx + ny - 2
    EqualVarianceTTest(nx, ny, xbar, df, stderr, t, μ0)
end

## TWO SAMPLE T-TEST (UNEQUAL VARIANCE)

struct UnequalVarianceTTest <: TwoSampleTTest
    n_x::Int     # number of observations
    n_y::Int     # number of observations
    xbar::Real   # estimated mean
    df::Real     # degrees of freedom
    stderr::Real # empirical standard error
    t::Real      # t-statistic
    μ0::Real     # mean under h_0
end

testname(::UnequalVarianceTTest) = "Two sample t-test (unequal variance)"

"""
    UnequalVarianceTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real})

Perform an unequal variance two-sample t-test of the null hypothesis that `x` and `y` come
from distributions with equal means against the alternative hypothesis that the
distributions have different means.

This test is sometimes known as Welch's t-test. It differs from the equal variance t-test in
that it computes the number of degrees of freedom of the test using the Welch-Satterthwaite
equation:
```math
    ν_{χ'} ≈ \\frac{\\left(\\sum_{i=1}^n k_i s_i^2\\right)^2}{\\sum_{i=1}^n
        \\frac{(k_i s_i^2)^2}{ν_i}}
```

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function UnequalVarianceTTest(x::AbstractVector{T}, y::AbstractVector{S}, μ0::Real=0) where {T<:Real,S<:Real}
    nx, ny = length(x), length(y)
    xbar = mean(x)-mean(y)
    varx, vary = var(x), var(y)
    stderr = sqrt(varx/nx + vary/ny)
    t = (xbar-μ0)/stderr
    df = (varx / nx + vary / ny)^2 / ((varx / nx)^2 / (nx - 1) + (vary / ny)^2 / (ny - 1))
    UnequalVarianceTTest(nx, ny, xbar, df, stderr, t, μ0)
end


text:
using Plots

# enable plotting by default
if !@isdefined do_visu; do_visu = true end

@views function diffusion_1D_nprocs(; do_visu=true)
    # Physics
    lx  = 10.0
    λ   = 1.0
    nt  = 200
    # Numerics
    np  = 4             # number of procs
    nx  = 32            # local number of gridpoints
    # Derived numerics
    nxg = (nx-2)*np+2   # global number of grid points
    dxg = lx/nxg        # dx for global grid
    dt  = dxg^2/λ/2.1
    # Array allocation
    x   = zeros(nx,np)  # local coord array
    H   = zeros(nx,np)  # local H array
    xt  = zeros(nxg)    # global coord array
    Ht  = zeros(nxg)    # global initial H array
    Hg  = zeros(nxg)    # global H array
    # Initial condition
    for ip = 1:np
        i1 = 1 + (ip-1)*(nx-2)
        for ix = 1:nx
            x[ix,ip] = ( (ip-1)*(nx-2) + (ix-0.5) )*dxg - 0.5*lx
            H[ix,ip] = exp(-x[ix,ip]^2)
        end
        xt[i1:i1+nx-2] .= x[1:end-1,ip]; if (ip==np) xt[i1+nx-1] = x[end,ip] end
        Ht[i1:i1+nx-2] .= H[1:end-1,ip]; if (ip==np) Ht[i1+nx-1] = H[end,ip] end
    end
    # Time loop
    for it = 1:nt
        for ip = 1:np # compute physics locally
            H[2:end-1,ip] .= H[2:end-1,ip] + dt*λ*diff(diff(H[:,ip])/dxg)/dxg
        end
        for ip = 1:np-1 # update boundaries
            H[end,ip  ] = H[    2,ip+1]
            H[  1,ip+1] = H[end-1,ip  ]
        end
        for ip = 1:np # global picture
            i1 = 1 + (ip-1)*(nx-2)
            Hg[i1:i1+nx-2] .= H[1:end-1,ip]
        end
        # Visualise
        if do_visu
            fontsize = 12
            plot(xt, Ht, legend=false, linewidth=1, markershape=:circle, markersize=3, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"), titlefontsize=fontsize, titlefont="Courier")
            display(plot!(xt, Hg, legend=false, linewidth=3, framestyle=:box, xlabel="Lx", ylabel="H", title="diffusion (it=$(it))"))
        end
    end
end

diffusion_1D_nprocs(; do_visu=do_visu)


text:
# ---
# title: 59. Spiral Matrix II
# id: problem59
# author: Qling
# date: 2021-03-16
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/spiral-matrix-ii/description/>
# hidden: true
# ---
#
# Given a positive integer `n`, generate an `n x n` `matrix` filled with
# elements from `1` to `n2` in spiral order.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/spiraln.jpg)
#
#
#
#     Input: n = 3
#     Output: [[1,2,3],[8,9,4],[7,6,5]]
#
#
# **Example 2:**
#
#
#
#     Input: n = 1
#     Output: [[1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 20`
#
#
## @lc code=start
using LeetCode

function generate_matrix(n::Int)::Vector{Vector{Int}}
    mat = [fill(0, n) for _ = 1:n]
    i, j, di, dj = 1, 1, 0, 1

    for k = 1:n*n
        mat[i][j] = k
        if mat[mod1(i + di, n)][mod1(j + dj, n)] != 0
            di, dj = dj, -di
        end

        i += di
        j += dj
    end

    return mat
end
## @lc code=end


text:
# Beta distribution

export Beta

@parameterized Beta(α,β)

@kwstruct Beta(α, β)

@kwalias Beta [
    a     => α
    alpha => α
    b     => β
    beta  => β
]

TV.as(::Beta) = as𝕀

function logdensity(d::Beta{(:α, :β), Tuple{A,B}}, x::X) where {A,B,X}
    if static_hasmethod(xlogy, Tuple{A,X}) && static_hasmethod(xlog1py, Tuple{B,X})
        return xlogy(d.α - 1, x) + xlog1py(d.β - 1, -x)
    else
        return (d.α - 1) * log(x) + (d.β - 1) * log1p(-x)
    end
end

function basemeasure(d::Beta{(:α,:β)})
    inbounds(x) = 0 < x < 1
    constℓ = 0.0
    varℓ() = - logbeta(d.α, d.β)
    base = Lebesgue(ℝ)
    FactoredBase(inbounds, constℓ, varℓ, base)
end

Base.rand(rng::AbstractRNG, T::Type, μ::Beta) = rand(rng, Dists.Beta(μ.α, μ.β))

distproxy(d::Beta{(:α, :β)}) = Dists.Beta(d.α, d.β)

asparams(::Type{<:Beta}, ::Val{:α}) = asℝ₊
asparams(::Type{<:Beta}, ::Val{:β}) = asℝ₊


text:
@testset "SeqSim" begin
  Random.seed!(1234)
  sdata = georef((z=[1.,0.,1.],), [(25.,25.), (50.,75.), (75.,50.)])
  sgrid = CartesianGrid(100,100)

  prob1 = SimulationProblem(sgrid, :z => Float64, 3)
  prob2 = SimulationProblem(sdata, sgrid, :z, 3)

  solver = SeqSim(:z => (estimator=DummyEstimator(),
                         neighborhood=NormBall(10.),
                         minneighbors=1, maxneighbors=10,
                         marginal=Normal(), path=RandomPath(),
                         mapping=NearestMapping()))

  Random.seed!(1234)
  usol = solve(prob1, solver)
  csol = solve(prob2, solver)

  # reals = csol[:z]
  # inds = LinearIndices(size(sgrid))
  # @test all(reals[i][inds[25,25]] == 1. for i in 1:3)
  # @test all(reals[i][inds[50,75]] == 0. for i in 1:3)
  # @test all(reals[i][inds[75,50]] == 1. for i in 1:3)
end


text:
module Metis64

    if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
        include("../deps/deps.jl")      # define libmetis and check it can be dlopen'd
    else
        error("Metis package not properly installed. Please run Pkg.build(\"Metis\")")
    end

    using Graphs,Compat                 # for AdjacencyList types
    using LightGraphs                   # metisform
    export
        partGraphKway,
        partGraphRecursive

    include("metis_h.jl")               # define constants
    include("util.jl")                  # metisform and testgraph functions

# """
# ```
# partGraphKway(G, nparts,
#               [, adjwgt, vwgt, vsize,
#                tpwgts, ubvec, options]) -> objval, part
# ```
#
# Partition a graph with multilevel k-way partitioning. See Metis
# documentation for more detail.
#
# Inputs:
#
# * `G` : Undirected graph. Can be an adjacency matrix (sparse or dense), a Graphs package adjacency list, or a LightGraphs package graph.
# * `nparts` : The number of parts to partition the graph.
# * `adjwgt` : Boolean indicating if edges are weighted. If true, `G` must be an adjacency matrix.
# * `vwgt` : The weights of the vertices.
# * `vsize` : The size of the vertices for computing the total communication volume.
# * `tpwgts` : Array of size `nparts` x `ncon` that specifies the desired weight for each partition and constraint.
# * `ubvec` : Array of size `ncon` that specifies the allowed load imbalance tolerance for each constraint.
# * `options` : Array of options.
#
# Outputs:
#
# * `objval` : The edge-cut or the total communication volume of the partitioning solution.
# * `part` : The partition vector of the graph.
#
# """
function partGraphKway(G, nparts, ::Type{T};
                           options::Array=-ones(METIS_NOPTIONS),
                           adjwgt::Bool=false,
                           vwgt::Array=[], vsize::Array=[],
                           tpwgts::Array=[], ubvec::Array=[],
                           ) where {T<:Integer}

    if !(T in [Int32 Int64])
        error("nparts must be either Int32 or Int64, but is " * string(T))
    end

    # Get adjacency structure of graph
    if adjwgt
        typeof(G)<:SparseMatrixCSC || error("weighted graphs must be represented in CSC format")
        n, xadj, adjncy, _adjwgt = metisform_weighted(G, T)
    else
         n, xadj, adjncy = metisform(G, T)
        _adjwgt = C_NULL
    end

    # Check parameters are valid
    nparts > 1 || error("nparts must be greater than one")
    length(vwgt)==0 || size(vwgt)==(n,) || error("vwgt must have n entries")
    length(vsize)==0 || size(vsize)==(n,) || error("vsize must have n entries")
    if length(tpwgts)>0
        ncon = convert(T, size(tpwgts)[1])
        size(tpwgts)==(ncon,nparts) || error("tpwgts must be an ncon x nparts array")
    else
        ncon = one(T)
    end
    length(ubvec)==0 || size(ubvec)==(ncon,) || error("ubvec must have ncon entries")
    length(options)==0 || size(options)==(METIS_NOPTIONS,) || error("options must have METIS_NOPTIONS entries")

    # Initialize partition parameters
    _vwgt = (length(vwgt)>0) ? round.(T, vwgt) : C_NULL
    _vsize = (length(vsize)>0) ? round.(T, vsize) : C_NULL

    _tpwgts = (length(tpwgts)>0) ? convert(Array{T==Int32 ? Cfloat : Cdouble}, tpwgts) : C_NULL
    _ubvec = (length(ubvec)>0) ? convert(Array{T==Int32 ? Cfloat : Cdouble}, ubvec) : C_NULL

    _options = (length(options)>0) ? convert(Array{T}, options) : C_NULL

    # Allocate memory for outputs
    part = Array{T}(n)
    objval = zeros(T, 1)

    #println("here, we go")
    # Call Metis partitioner
    if T==Int32
        err = ccall((:METIS_PartGraphKway, libmetis), Cint, # Last one is return value
                (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                 Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cfloat}, Ptr{Cfloat},
                 Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                &n, &ncon, xadj, adjncy,
               _vwgt, _vsize, _adjwgt, &convert(T,nparts),
               _tpwgts, _ubvec, _options, objval, part)
    else
        err = ccall((:METIS_PartGraphKway, libmetis64), Clonglong, # Last one is return value
                (Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong},
                 Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble},
                 Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}),
                &n, &ncon, xadj, adjncy,
                _vwgt, _vsize, _adjwgt, &convert(T,nparts),
                _tpwgts, _ubvec, _options, objval, part)
    end

    part = part .+ 1

    # Return results
    return err, objval[1], part

end

function partGraphRecursive(G, nparts, ::Type{T};
                               options::Array=-ones(METIS_NOPTIONS),
                               adjwgt::Bool=false,
                               vwgt::Array=[], vsize::Array=[],
                               tpwgts::Array=[], ubvec::Array=[],
                               ) where {T<:Integer}
    if !(T in [Int32 Int64])
        error("nparts must be either Int32 or Int64")
    end

    # Get adjacency structure of graph
    if adjwgt
        typeof(G)<:SparseMatrixCSC || error("weighted graphs must be represented in CSC format")
        n, xadj, adjncy, _adjwgt = metisform_weighted(G, T)
    else
        n, xadj, adjncy = metisform(G, T)
        _adjwgt = C_NULL
    end

    # Check parameters are valid
    nparts > 1 || error("nparts must be greater than one")
    length(vwgt)==0 || size(vwgt)==(n,) || error("vwgt must have n entries")
    length(vsize)==0 || size(vsize)==(n,) || error("vsize must have n entries")
    if length(tpwgts)>0
        ncon = convert(T, size(tpwgts)[1])
        size(tpwgts)==(ncon,nparts) || error("tpwgts must be an ncon x nparts array")
    else
        ncon = one(T)
    end
    length(ubvec)==0 || size(ubvec)==(ncon,) || error("ubvec must have ncon entries")
    length(options)==0 || size(options)==(METIS_NOPTIONS,) || error("options must have METIS_NOPTIONS entries")

    # Initialize partition parameters
    _vwgt = (length(vwgt)>0) ? round.(T, vwgt) : C_NULL
    _vsize = (length(vsize)>0) ? round.(T, vsize) : C_NULL

    _tpwgts = (length(tpwgts)>0) ? convert(Array{T==Int32 ? Cfloat : Cdouble}, tpwgts) : C_NULL
    _ubvec = (length(ubvec)>0) ? convert(Array{T==Int32 ? Cfloat : Cdouble}, ubvec) : C_NULL

    _options = (length(options)>0) ? convert(Array{T}, options) : C_NULL

    part = Array{T}(n)
    objval = zeros(T, 1)
    if T==Int32
        err = ccall((:METIS_PartGraphRecursive, libmetis), Cint, # Last one is return value
                (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                 Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cfloat}, Ptr{Cfloat},
                 Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                &n, &ncon, xadj, adjncy,
               _vwgt, _vsize, _adjwgt, &convert(T,nparts),
               _tpwgts, _ubvec, _options, objval, part)
    else
        err = ccall((:METIS_PartGraphRecursive, libmetis64), Clonglong, # Last one is return value
                (Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong},
                 Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble},
                 Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Clonglong}),
                &n, &ncon, xadj, adjncy,
                _vwgt, _vsize, _adjwgt, &convert(T, nparts),
                _tpwgts, _ubvec, _options, objval, part)
    end

    part = part .+ 1

    return err, objval[1], part
end

end # module


text:
using Documenter, TaylorSeries

makedocs(
    modules = Module[TaylorSeries],
    doctest = false,
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-cinder", "python-markdown-math"),
    repo   = "github.com/JuliaDiff/TaylorSeries.jl.git",
    osname = "linux"
)